Error in Double Integration problem

%anisotropic spatial power spectrum
syms q theta

qx=q*cos(theta);qy=q*sin(theta);
q=sqrt(qx^2+qy^2);

mu=2;gamma=90;C0=0.72;C1=2.35;epsilon=10^-1;chiT=10^-7;

S = [32.4 33 35 34.3 32];
T=[23 25 24 22 21];
z=[11 12 13 14 15 ];

t=T;SP=S;
p=z.*9.8.*.1045;


alpha =  0.02;
betac =0.222;
H = diff(T)./diff(S);
H = [H(1) H];
w = (alpha.*H)./(betac);

clear dr
for i=1:length(w)
   
    if abs(w(i))>=1
        dr(i) = abs(w(i)) + (abs(w(i)).^0.5)*(abs(w(i))-1).^0.5;
    elseif  abs(w(i))< 0.5
        dr(i)=0.15.*abs(w(i));
    else  
        dr(i)=1.85.*abs(w(i))-0.85;
    end
    
end

%%
M1=1.541 + (1.998*10^-2).*T - (9.52*10^-5).*T.^2;
M2=7.974 - 7.561*10^-2.*T + 4.724*10^-4.*T.^2;
M3=0.157.*(T + 64.993).^2;
mu0=(1+M1.*S+M2.*S.^2).*((4.2844*10^-5) + (M3.^-1));

a1 = 9.999 * 10^2;a2= 2.034 * 10^-2;a3=-6.162 * 10^-3;a4= 2.261 * 10^-5;a5= -4.657 * 10^-8;
b1=8.020 * 10^2;b2=-2.001;b3= 1.677 * 10^-2;b4= -3.060 * 10^-5;b5= -1.613 * 10^-5;
R1=a1+a2.*T+a3.*T.^2+a4.*T.^3+a5.*T.^4;
R2=b1.*S+b2.*S.*T+b3.*S.*T.^2+b4.*S.*T.^3+b5.*S.^2.*T.^2;
rho0=R1+R2;
nu=mu0/rho0;
 
C = 5.328 - 9.76 * 10^-2.*S + 4.04 * 10.^-4.*S.^2;
D = -6.913 * 10^-3 + 7.351 * 10^-4.*S - 3.15 * 10^-6.*S.^2;
E = 9.6 * 10^-6 - 1.927 * 10^-6*S + 8.23 * 10^-9*S.^2;
F = 2.5 * 10^-9 + 1.666 * 10^-9*S - 7.125 * 10^-12*S.^2;
c0 = C + D.*(T - 273.15) + E.*(T - 273.15).^2 + F.*(T - 273.15).^3;
K1=(343.5 +0.037*S)/(T+273.15);
K2=(T+273.15)/(647+0.03*S);
KK=log10(240+0.0002*S) + 0.434*(2.3-K1)*(1-K2)^0.333;
K=10.^KK;
DT=K./(rho0.*c0);DS=0.01.*DT;
Pt = nu./DT;Ps = nu./DS;Pts=(Pt+Ps)/2;
etta = (nu^3/epsilon)^(1/4);
mux =sqrt(((mu^2) + (tan(gamma))^2)/(1 + (tan(gamma))^2));
muy =sqrt(((mu^2) + (tan(gamma))^2)/(1 + (mu^2)*(tan(gamma))^2));
deltaq=(3/2)*C1^2*((etta*q)^(4/3)) + C1^3*(etta*q)^2;

A=((alpha.^2).*chiT.*mux.*muy)./(4*pi*w.^2.*(epsilon.^(1/3)).*(q.^(11./3)));
B=1+C1.*(etta.^(2/3)).*(q.^(2./3));
D1=(w.^2).*exp(-C0.*deltaq./(C1.^2.*Pt));
D2=dr.*exp(-C0.*deltaq./(C1.^2.*Ps));
D3=w.*(dr+1).*exp(-C0.*deltaq./(2.*C1.^2.*Pts));
Phi = q.*(C0.*A.*B.*(D1+D2-D3))./(mux*muy);

fun = matlabFunction(Phi,'Vars',[q,theta]);

Error using sym/matlabFunction>checkVars
Variable names must be valid MATLAB variable names.

Error in sym/matlabFunction (line 158)
vars = checkVars(funvars,opts);

Pho=((integral2(fun,0,Inf,0,2*pi)).^(0.5));

Answer 1

[Torsten]

You did the same mistake as in the code I corrected before. "fun" must return a scalar if called by a value pair (q,theta). You function "fun" returns an array of length 5.

Maybe you mean

for i=1:5
  fun = matlabFunction(Phi(i));
  Pho(i)=((integral2(fun,0,Inf,0,2*pi)).^(0.5));
end

[Athira T Das]

@TorstenChanged my code accordingly by using loop, but the integration is unsucessful.

As you corrected my former question by including a third variable is not possible in this integration.

[Athira T Das]

I changed the code with

clc;close all;clear all;
mu=2;gamma=90;C0=0.72;C1=2.35;epsilon=10^-1;chiT=10^-7;
mux =sqrt(((mu^2) + (tan(gamma))^2)/(1 + (tan(gamma))^2));
muy =sqrt(((mu^2) + (tan(gamma))^2)/(1 + (mu^2)*(tan(gamma))^2));

syms q theta 
kx=q*cos(theta)/mux;
ky=q*sin(theta)/muy;

S = [32.4 33 35 34.3 32];
T=[23 25 24 22 21];
z=[11 12 13 14 15 ];

t=T;SP=S;
p=z.*9.8.*.1045;
long=90;lat=8;

alpha =  0.02;
betac =0.222;
H = diff(T)./diff(S);
H = [H(1) H];
w = (alpha.*H)./(betac);

clear dr
for i=1:length(w)
   
    if abs(w(i))>=1
        dr(i) = abs(w(i)) + (abs(w(i)).^0.5)*(abs(w(i))-1).^0.5;
    elseif  abs(w(i))< 0.5
        dr(i)=0.15.*abs(w(i));
    else  
        dr(i)=1.85.*abs(w(i))-0.85;
    end
    
end

%%
M1=1.541 + (1.998*10^-2).*T - (9.52*10^-5).*T.^2;
M2=7.974 - 7.561*10^-2.*T + 4.724*10^-4.*T.^2;
M3=0.157.*(T + 64.993).^2;
mu0=(1+M1.*S+M2.*S.^2).*((4.2844*10^-5) + (M3.^-1));

a1 = 9.999 * 10^2;a2= 2.034 * 10^-2;a3=-6.162 * 10^-3;a4= 2.261 * 10^-5;a5= -4.657 * 10^-8;
b1=8.020 * 10^2;b2=-2.001;b3= 1.677 * 10^-2;b4= -3.060 * 10^-5;b5= -1.613 * 10^-5;
R1=a1+a2.*T+a3.*T.^2+a4.*T.^3+a5.*T.^4;
R2=b1.*S+b2.*S.*T+b3.*S.*T.^2+b4.*S.*T.^3+b5.*S.^2.*T.^2;
rho0=R1+R2;
nu=mu0/rho0;

C = 5.328 - 9.76 * 10^-2.*S + 4.04 * 10.^-4.*S.^2;
D = -6.913 * 10^-3 + 7.351 * 10^-4.*S - 3.15 * 10^-6.*S.^2;
E = 9.6 * 10^-6 - 1.927 * 10^-6*S + 8.23 * 10^-9*S.^2;
F = 2.5 * 10^-9 + 1.666 * 10^-9*S - 7.125 * 10^-12*S.^2;
c0 = C + D.*(T - 273.15) + E.*(T - 273.15).^2 + F.*(T - 273.15).^3;
K1=(343.5 +0.037*S)/(T+273.15);
K2=(T+273.15)/(647+0.03*S);
KK=log10(240+0.0002*S) + 0.434*(2.3-K1)*(1-K2)^0.333;
K=10.^KK;
DT=K./(rho0.*c0);DS=0.01.*DT;
Pt = nu./DT;Ps = nu./DS;Pts=(Pt+Ps)/2;
etta = (nu^3/epsilon)^(1/4);


deltaa=(3/2)*C1^2*((etta)^(4/3))*(((mux*kx)^2)+((muy*ky)^2))^(2/3) + C1^3*(etta)^2*(((mux*kx)^2)+((muy*ky)^2))^(2);

A=((alpha.^2).*chiT.*mux.*muy)./(4*pi*w.^2.*(epsilon.^(1/3)).*((((mux*kx)^2)+((muy*ky)^2))^(2/3)));
B=1+C1.*(etta.^(2/3)).*((((mux*kx)^2)+((muy*ky)^2))^(1/3));
D1=(w.^2).*exp(-C0.*deltaa./(C1.^2.*Pt));
D2=dr.*exp(-C0.*deltaa./(C1.^2.*Ps));
D3=w.*(dr+1).*exp(-C0.*deltaa./(2.*C1.^2.*Pts));
Phi = q.*(C0.*A.*B.*(D1+D2-D3))./(mux*muy);

fun = matlabFunction(Phi,'Vars',[q,theta]);
Pho=((integral2(fun,0,Inf,0,2*pi)).^(0.5));

Error using integralCalc/finalInputChecks
Output of the function must be the same size as the input. If FUN is an array-valued integrand, set the 'ArrayValued' option to true.

Error in integralCalc/iterateScalarValued (line 315)
finalInputChecks(x,fx);

Error in integralCalc/vadapt (line 132)
[q,errbnd] = iterateScalarValued(u,tinterval,pathlen);

Error in integralCalc (line 75)
[q,errbnd] = vadapt(@AtoBInvTransform,interval);

Error in integral2Calc>@(xi,y1i,y2i)integralCalc(@(y)fun(xi*ones(size(y)),y),y1i,y2i,opstruct.integralOptions) (line 17)
innerintegral = @(x)arrayfun(@(xi,y1i,y2i)integralCalc( ...

Error in integral2Calc>@(x)arrayfun(@(xi,y1i,y2i)integralCalc(@(y)fun(xi*ones(size(y)),y),y1i,y2i,opstruct.integralOptions),x,ymin(x),ymax(x)) (line 17)
innerintegral = @(x)arrayfun(@(xi,y1i,y2i)integralCalc( ...

Error in integralCalc/iterateScalarValued (line 314)
fx = FUN(t);

Error in integralCalc/vadapt (line 132)
[q,errbnd] = iterateScalarValued(u,tinterval,pathlen);

Error in integralCalc (line 83)
[q,errbnd] = vadapt(@AToInfInvTransform,interval);

Error in integral2Calc>integral2i (line 20)
[q,errbnd] = integralCalc(innerintegral,xmin,xmax,opstruct.integralOptions);

Error in integral2Calc (line 7)
[q,errbnd] = integral2i(fun,xmin,xmax,ymin,ymax,optionstruct);

Error in integral2 (line 105)
Q = integral2Calc(fun,xmin,xmax,yminfun,ymaxfun,opstruct);

[Walter Roberson]

clc;close all;clear all;
mu=2;gamma=90;C0=0.72;C1=2.35;epsilon=10^-1;chiT=10^-7;
mux =sqrt(((mu^2) + (tan(gamma))^2)/(1 + (tan(gamma))^2));
muy =sqrt(((mu^2) + (tan(gamma))^2)/(1 + (mu^2)*(tan(gamma))^2));

syms q theta 
kx=q*cos(theta)/mux;
ky=q*sin(theta)/muy;

S = [32.4 33 35 34.3 32];
T=[23 25 24 22 21];
z=[11 12 13 14 15 ];

t=T;SP=S;
p=z.*9.8.*.1045;
long=90;lat=8;

alpha =  0.02;
betac =0.222;
H = diff(T)./diff(S);
H = [H(1) H];
w = (alpha.*H)./(betac);

clear dr
for i=1:length(w)
   
    if abs(w(i))>=1
        dr(i) = abs(w(i)) + (abs(w(i)).^0.5)*(abs(w(i))-1).^0.5;
    elseif  abs(w(i))< 0.5
        dr(i)=0.15.*abs(w(i));
    else  
        dr(i)=1.85.*abs(w(i))-0.85;
    end
    
end

%%
M1=1.541 + (1.998*10^-2).*T - (9.52*10^-5).*T.^2;
M2=7.974 - 7.561*10^-2.*T + 4.724*10^-4.*T.^2;
M3=0.157.*(T + 64.993).^2;
mu0=(1+M1.*S+M2.*S.^2).*((4.2844*10^-5) + (M3.^-1));

a1 = 9.999 * 10^2;a2= 2.034 * 10^-2;a3=-6.162 * 10^-3;a4= 2.261 * 10^-5;a5= -4.657 * 10^-8;
b1=8.020 * 10^2;b2=-2.001;b3= 1.677 * 10^-2;b4= -3.060 * 10^-5;b5= -1.613 * 10^-5;
R1=a1+a2.*T+a3.*T.^2+a4.*T.^3+a5.*T.^4;
R2=b1.*S+b2.*S.*T+b3.*S.*T.^2+b4.*S.*T.^3+b5.*S.^2.*T.^2;
rho0=R1+R2;
nu=mu0/rho0;

C = 5.328 - 9.76 * 10^-2.*S + 4.04 * 10.^-4.*S.^2;
D = -6.913 * 10^-3 + 7.351 * 10^-4.*S - 3.15 * 10^-6.*S.^2;
E = 9.6 * 10^-6 - 1.927 * 10^-6*S + 8.23 * 10^-9*S.^2;
F = 2.5 * 10^-9 + 1.666 * 10^-9*S - 7.125 * 10^-12*S.^2;
c0 = C + D.*(T - 273.15) + E.*(T - 273.15).^2 + F.*(T - 273.15).^3;
K1=(343.5 +0.037*S)/(T+273.15);
K2=(T+273.15)/(647+0.03*S);
KK=log10(240+0.0002*S) + 0.434*(2.3-K1)*(1-K2)^0.333;
K=10.^KK;
DT=K./(rho0.*c0);DS=0.01.*DT;
Pt = nu./DT;Ps = nu./DS;Pts=(Pt+Ps)/2;
etta = (nu^3/epsilon)^(1/4);


deltaa=(3/2)*C1^2*((etta)^(4/3))*(((mux*kx)^2)+((muy*ky)^2))^(2/3) + C1^3*(etta)^2*(((mux*kx)^2)+((muy*ky)^2))^(2);

A=((alpha.^2).*chiT.*mux.*muy)./(4*pi*w.^2.*(epsilon.^(1/3)).*((((mux*kx)^2)+((muy*ky)^2))^(2/3)));
B=1+C1.*(etta.^(2/3)).*((((mux*kx)^2)+((muy*ky)^2))^(1/3));
D1=(w.^2).*exp(-C0.*deltaa./(C1.^2.*Pt));
D2=dr.*exp(-C0.*deltaa./(C1.^2.*Ps));
D3=w.*(dr+1).*exp(-C0.*deltaa./(2.*C1.^2.*Pts));
Phi = q.*(C0.*A.*B.*(D1+D2-D3))./(mux*muy);

Phi = Phi(:)

Phi = 

limit(Phi, q, 0, 'left')

ans = 

limit(Phi, q, 0, 'right')

ans = 

The integral over q starts at 0 but there is an irreducible singularity there.

If you look at the display of the symbolic formulas, and go down to then with q being a linear factor for both terms of the addend, would be 0 when q is 0. is involved in the other σ terms as a linear factor, so they are all 0. And if you look at the formulas for the expressions, then all involve on the numerator and on the denominator, so they all involve 0 / 0

[Athira T Das]

@Walter Roberson I removed the singularity by modifying my equation;

clc;close all;clear all;

mu=2;gamma=pi/2;C0=0.72;C1=2.35;epsilon=10^-2;chiT=10^-6;
mux =sqrt(((mu^2) + (tan(gamma))^2)/(1 + (tan(gamma))^2));
muy =sqrt(((mu^2) + (tan(gamma))^2)/(1 + (mu^2)*(tan(gamma))^2));

syms q theta 
kx=q*cos(theta)/mux;
ky=q*sin(theta)/muy;
S=34;T=25;
w = -3;
dr = abs(w) + (abs(w).^0.5)*(abs(w)-1).^0.5;
%%
M1=1.541 + (1.998*10^-2).*T - (9.52*10^-5).*T.^2;
M2=7.974 - 7.561*10^-2.*T + 4.724*10^-4.*T.^2;
M3=0.157.*(T + 64.993).^2;
mu0=(1+M1.*S+M2.*S.^2).*((4.2844*10^-5) + (M3.^-1));

a1 = 9.999 * 10^2;a2= 2.034 * 10^-2;a3=-6.162 * 10^-3;a4= 2.261 * 10^-5;a5= -4.657 * 10^-8;
b1=8.020 * 10^2;b2=-2.001;b3= 1.677 * 10^-2;b4= -3.060 * 10^-5;b5= -1.613 * 10^-5;
R1=a1+a2.*T+a3.*T.^2+a4.*T.^3+a5.*T.^4;
R2=b1.*S+b2.*S.*T+b3.*S.*T.^2+b4.*S.*T.^3+b5.*S.^2.*T.^2;
rho0=R1+R2;
nu=mu0/rho0;
%%to calculate Prandtl numbers
C = 5.328 - 9.76 * 10^-2.*S + 4.04 * 10.^-4.*S.^2;
D = -6.913 * 10^-3 + 7.351 * 10^-4.*S - 3.15 * 10^-6.*S.^2;
E = 9.6 * 10^-6 - 1.927 * 10^-6*S + 8.23 * 10^-9*S.^2;
F = 2.5 * 10^-9 + 1.666 * 10^-9*S - 7.125 * 10^-12*S.^2;
c0 = C + D.*(T - 273.15) + E.*(T - 273.15).^2 + F.*(T - 273.15).^3;
K1=(343.5 +0.037*S)/(T+273.15);
K2=(T+273.15)/(647+0.03*S);
KK=log10(240+0.0002*S) + 0.434*(2.3-K1)*(1-K2)^0.333;
K=10.^KK;
DT=K./(rho0.*c0);DS=0.01.*DT;
Pt = nu./DT;Ps = nu./DS;Pts=(Pt+Ps)/2;
etta = (nu^3/epsilon)^(1/4);


deltaa=(3/2)*C1^2*((etta)^(4/3))*(((mux*kx)^2)+((muy*ky)^2))^(2/3) + C1^3*(etta)^2*(((mux*kx)^2)+((muy*ky)^2))^(2);
alpha=0.004;
A=((alpha.^2).*chiT.*mux.*muy)./(4*pi*w.^2.*(epsilon.^(1/3)).*((((mux*kx)^2)+((muy*ky)^2))^(2/3)));
B=1+C1.*(etta.^(2/3)).*((((mux*kx)^2)+((muy*ky)^2))^(1/3));
D1=(w.^2).*exp(-C0.*deltaa./(C1.^2.*Pt));
D2=dr.*exp(-C0.*deltaa./(C1.^2.*Ps));
D3=w.*(dr+1).*exp(-C0.*deltaa./(2.*C1.^2.*Pts));
z=4000;
lambda=530;
k = (2.*pi)./(lambda);
Xx=(((pi.^2).*(k.^2).*z)./3);
Phi = Xx.*q^3.*(C0.*A.*B.*(D1+D2-D3))./(mux*muy);
limit(Phi, q, 0, 'left')

ans = 

0

limit(Phi, q, 0, 'right')

ans = 

0

 
fun = matlabFunction(Phi,'Vars',[q,theta]);
Pho=((integral2(fun,0,Inf,0,2*pi)).^(0.5));

Warning: Inf or NaN value encountered.

Warning: The integration was unsuccessful.

[Torsten]

Plot a slice of your function at theta = 0.2, e.g., by adding the following lines to your code:

theta = 0.2;
q = 0:0.1:12;
plot(q,fun(q,theta))

Seems your function does not converge to 0 as q -> Inf, but blows up very fast.

Maybe you forgot a minus sign in some exponential.

[Athira T Das]

@Torsten It works finally!!!

Answer 2

[Walter Roberson]

syms q theta

q becomes a scalar symbolic variable

q=sqrt(qx^2+qy^2);

but now it is an expression

fun = matlabFunction(Phi,'Vars',[q,theta]);

and now you try to use it as a variable name in creating a function. The expression being turned into a function does not have any variable named q in it -- once q is turned into an expression, whever q is used, it gets expanded to the expression.

[Torsten]

Nevertheless, symvar(Phi) returns q and theta...

[Walter Roberson]

I see what you mean, and

fun = matlabFunction(Phi,'Vars',[sym('q'),theta]);

can be used -- it might not give the expected result since q means two different things, but it would at least pass that step.

But then it would run into the problem that you pointed out, that Phi is non-scalar.

[Torsten]

So you think the result of the integration will be different if one would have set

qi = sqrt(qx^2+qy^2);

instead of

q = sqrt(qx^2+qy^2);

and replaced all q's in the sequel by "qi" ?

That would be very bad on my opinion.

[Walter Roberson]

In my opinion, if q is an expression or function, then requesting to integrate over q is a request for a path integral. Which should probably be dealt with by a change of variables.

If a path integral is not what is desired, then in my opinion do not write the expressions in such a way that you need to sym() up a variable name because you have aliased a variable. If path integral is what is desired, then deal with it properly.

[Torsten]

Doesn't seem to influence the function to be integrated.

syms x
x = x^2;
f = sin(x);
fnum = matlabFunction(f,'Vars',sym('x'));
v = integral(fnum,0,2*pi)

v = 0.6421

syms y
yv = y^2;
g = sin(yv);
gnum = matlabFunction(g,'Vars',y);
w = integral(gnum,0,2*pi)

w = 0.6421

[Walter Roberson]

I would put it to you that asking to integrate

syms x
x = x^2
f = sin(x)

integrate f over x, is a request to integrate over the path x^2 -- that morally you are asking for

rather than

In my opinion, asking to integrate with respect to a variable should always refer to the current value of the variable, not with respect to some value the variable might previously have had.