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]);
Pho=((integral2(fun,0,Inf,0,2*pi)).^(0.5));
回答（2 个）
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 个评论
Athira T Das
2023-4-24
编辑：Athira T Das
2023-4-24
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(:)
limit(Phi, q, 0, 'left')
limit(Phi, q, 0, 'right')
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
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
Torsten
2023-4-24
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
Walter Roberson
2023-4-23
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 个评论
Torsten
2023-4-23
Nevertheless, symvar(Phi) returns q and theta...
Walter Roberson
2023-4-23
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.
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
2023-4-23
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.
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)
syms y
yv = y^2;
g = sin(yv);
gnum = matlabFunction(g,'Vars',y);
w = integral(gnum,0,2*pi)
Walter Roberson
2023-4-24
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.
类别
在 帮助中心 和 File Exchange 中查找有关 Code Performance 的更多信息
2023-4-23
2023-4-25
Community Treasure Hunt
Find the treasures in MATLAB Central and discover how the community can help you!
Start Hunting!
Translated by 
