Error in Double Integration problem

Question

Athira T Das 2023-4-23

0
链接

此问题的直接链接

https://ww2.mathworks.cn/matlabcentral/answers/1951833-error-in-double-integration-problem

评论： Athira T Das 2023-4-25

%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));

0 个评论
显示 -2更早的评论隐藏 -2更早的评论

请先登录，再进行评论。

请先登录，再回答此问题。

Answer 1

Torsten 2023-4-23

0
链接

此回答的直接链接

https://ww2.mathworks.cn/matlabcentral/answers/1951833-error-in-double-integration-problem#answer_1221513

移动：Torsten 2023-4-23

在 MATLAB Online 中打开

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

6 个评论
显示 4更早的评论隐藏 4更早的评论

Athira T Das 2023-4-24

编辑：Athira T Das 2023-4-24

在 MATLAB Online 中打开

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 2023-4-24

在 MATLAB Online 中打开

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

%%

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 2023-4-24

在 MATLAB Online 中打开

@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 2023-4-24

在 MATLAB Online 中打开

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 2023-4-25

@Torsten It works finally!!!

请先登录，再进行评论。

Answer 2

Walter Roberson 2023-4-23

0
链接

此回答的直接链接

https://ww2.mathworks.cn/matlabcentral/answers/1951833-error-in-double-integration-problem#answer_1221538

在 MATLAB Online 中打开

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.

6 个评论
显示 4更早的评论隐藏 4更早的评论

Torsten 2023-4-24

在 MATLAB Online 中打开

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 2023-4-24

在 MATLAB Online 中打开

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.

请先登录，再进行评论。