Error "Arrays have incompatible sizes for this operation" when try to integrate

Hello. This is my entire code, but the relevant part is at the end:

syms kb hbar T c l p k d n m xi q d;

T=300;

hbar=1.054571817*10^(-34);

kb=1.380649*10^(-23);

c=300000000;

d1=1e-7;

N=10;

M=300;

a=load('SiO2Tourne.dat');

w1=a(:,1);

ep1=a(:,2);

eprot1=@(omega) interp1(w1,ep1,omega); %% DIELECTRIC FUNCTION %%

xi_k=@(k) 2.*pi.*kb.*T.*k./hbar;

z=d1./N;

q0=@(xi,q) sqrt(xi.^2./c.^2+q.^2);

q1=@(xi,q) sqrt(xi.^2./c.^2.*eprot1(xi)+q.^2);

q0te=@(xi,q) q0(xi,q);

q1te=@(xi,q) q1(xi,q);

q0tm=@(xi,q) q0(xi,q);

q1tm=@(xi,q) q1(xi,q)./eprot1(xi);

ste0=@(xi,q) 1./[q1te(xi,q)+q0te(xi,q)].*...

[1.*(q0te(xi,q)-q1te(xi,q)) 2.*q1te(xi,q); 2.*q0te(xi,q) q1te(xi,q)-q0te(xi,q)];

ste=@(n,xi,q) 1./[2*q1te(xi,q)].*...

[0 2*exp(-q1(xi,q).*z).*q1te(xi,q); 2*exp(-q1(xi,q).*z).*q1te(xi,q) 0];

steN=@(xi,q) 1./[q0te(xi,q)+q1te(xi,q)].*...

[exp(-2*q1(xi,q).*z).*(q1te(xi,q)-q0te(xi,q)) 2*exp(-q1(xi,q).*z).*q0te(xi,q); 2*exp(-q1(xi,q).*z).*q1te(xi,q) q0te(xi,q)-q1te(xi,q)];

stm0=@(xi,q) 1./[q0tm(xi,q)+q1tm(xi,q)].*...

[q0tm(xi,q)-q1tm(xi,q) 2*q1tm(xi,q); 2*q0tm(xi,q) -q0tm(xi,q)+q1tm(xi,q)];

stm=@(n,xi,q) 1./[2*q1tm(xi,q)].*...

[0 2*exp(-q1(xi,q).*z).*q1tm(xi,q); 2*exp(-q1(xi,q).*z).*q1tm(xi,q) 0];

stmN=@(xi,q) 1./[q0tm(xi,q)+q1tm(xi,q)].*...

[exp(-2*q1(xi,q).*z).*(q1tm(xi,q)-q0tm(xi,q)) 2*exp(-q1(xi,q).*z).*q0tm(xi,q); 2*exp(-q1(xi,q).*z).*q1tm(xi,q) q0tm(xi,q)-q1tm(xi,q)];

f11=@(x,y) x(1,1)+x(1,2)./[1-y.*x(2,2)].*y.*x(2,1);

for m=0:M-1

temp_te=@(q) 1./(q0te(xi_k(m),q)+q1te(xi_k(m),q)).*(exp(-2*q1(xi_k(m),q).*z).*(q1te(xi_k(m),q)-q0te(xi_k(m),q)));

temp_tm=@(q) 1./(q0tm(xi_k(m),q)+q1tm(xi_k(m),q)).*(exp(-2*q1(xi_k(m),q).*z).*(q1tm(xi_k(m),q)-q0tm(xi_k(m),q)));

for n=N-1:-1:1

S11te=@(q) f11(ste(n,xi_k(m),q),temp_te(q));

temp_te=@(q) S11te(q);

S11tm=@(q) f11(stm(n,xi_k(m),q),temp_tm(q));

temp_tm=@(q) S11tm(q);

end

S11te=@(q) f11(ste0(xi_k(m),q),temp_te(q));

temp_te=@(q) S11te(q);

S11tm=@(q) f11(stm0(xi_k(m),q),temp_tm(q));

temp_tm=@(q) S11tm(q);

rte{m+1}=temp_te;

rtm{m+1}=temp_tm;

end

g=@(k,q,d) q.*q0(xi_k(k),q).*((exp(2.*q0(xi_k(k),q).*d)./(rte{k+1}(q).^2)-1).^(-1)+...

(exp(2.*q0(xi_k(k),q).*d)./(rtm{k+1}(q).^2)-1).^(-1));

add=@(k,d) quadgk(@(q) g(k,q,d),0,Inf);

Now, when I try to evaluate the function "add" for some specific values I get the following error message:

Arrays have incompatible sizes for this operation.

Error in

untitled>@(xi,q)1./[q1te(xi,q)+q0te(xi,q)].*[1.*(q0te(xi,q)-q1te(xi,q)),2.*q1te(xi,q);2.*q0te(xi,q),q1te(xi,q)-q0te(xi,q)]

(line 38)

ste0=@(xi,q) 1./[q1te(xi,q)+q0te(xi,q)].*...

Error in untitled>@(q)f11(ste0(xi_k(m),q),temp_te(q)) (line 78)

S11te=@(q) f11(ste0(xi_k(m),q),temp_te(q));

Error in untitled>@(q)S11te(q) (line 79)

temp_te=@(q) S11te(q);

Error in

untitled>@(k,q,d)q.*q0(xi_k(k),q).*((exp(2.*q0(xi_k(k),q).*d)./(rte{k+1}(q).^2)-1).^(-1)+(exp(2.*q0(xi_k(k),q).*d)./(rtm{k+1}(q).^2)-1).^(-1))

(line 95)

g=@(k,q,d) q.*q0(xi_k(k),q).*((exp(2.*q0(xi_k(k),q).*d)./(rte{k+1}(q).^2)-1).^(-1)+...

Error in @(q)g(k,q,d)

Error in quadgk/evalFun (line 330)

fx = FUN(x);

Error in quadgk/f2 (line 361)

[y,too_close] = evalFun(t2t);

Error in quadgk/vadapt (line 249)

[fx,too_close] = f(x);

Error in quadgk (line 196)

[q,errbnd] = vadapt(@f2,interval);

Error in @(k,d)quadgk(@(q)g(k,q,d),0,Inf)

Related documentation

The same thing happens using "integral" instead of "quadgk". I partly solved the issue using the "integral" function and putting the optional inputs "'arrayvalue',true". In such a way the code works but it is too slow for the next computations, so I would like to solve the problem more efficiently.

Answer 1

[Walter Roberson]

quadgk and integral pass in vectors for the first parameter.

@(xi,q)1./[q1te(xi,q)+q0te(xi,q)].*[1.*(q0te(xi,q)-q1te(xi,q)),2.*q1te(xi,q);2.*q0te(xi,q),q1te(xi,q)-q0te(xi,q)]

Because q is a vector, q1te and q0te are vectors so the left [] calculates a vector. The right [] thinks that it is building a 2x2 array but because q is a vector it is building something that is larger than 2x2... something that is not compatible for .* purposes with the left []

You need to use arrayvalued, or you need to reorganize your logic. For example perhaps you could build the right [] in the third dimension

[Silvio]

Sorry I don't get it. q is a sym, and q0, q1 etc are handle functions

[Walter Roberson]

q0te=@(xi,q) q0(xi,q);

Except for checking the number of arguments, that code is equivalent to

q0te = @(varargin) q0(varargin{1}, varargin{2} );

Within an anonymous function, any variable named in the @() list is replaced by a reference to the corresponding positional parameter. q is not a sym in that context: it is whatever is passed in in that position.

When you follow the chain of functions passing things around, the ultimate source of the corresponding value is

add=@(k,d) quadgk(@(q) g(k,q,d),0,Inf);

and quadgk is going to be passing in a numeric vector, not a sym.

If q were a sym in that code then you would be calculating a symbolic expression, but quadgk requires that the calculated value be numeric not symbolic.

[Walter Roberson]

syms kb hbar T c l p k d n m xi q d;
T=300;
hbar=1.054571817*10^(-34);

That code does not mean that T should be assigned symbolic 300. That code is equivalent to

T = sym('T') ;
T = 300;

which starts by assigning MATLAB-level T a reference to symbolic T that lives in the symbolic engine. Then it replaces MATLAB-level T with numeric 300, cutting the link to the T that lives in the symbolic engine.

[Walter Roberson]

Remember that MATLAB itself does not have type declarations. When you assign to an unindexed variable, the variable becomes whatever the type of the right hand side is.