Error in Double Integration problem
5 views (last 30 days)
Show older comments
%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));
0 Comments
Answers (2)
Torsten
on 23 Apr 2023
Moved: Torsten
on 23 Apr 2023
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 Comments
Torsten
on 24 Apr 2023
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.
Walter Roberson
on 23 Apr 2023
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 Comments
Torsten
on 24 Apr 2023
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
on 24 Apr 2023
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.
See Also
Categories
Find more on Calculus in Help Center and File Exchange
Community Treasure Hunt
Find the treasures in MATLAB Central and discover how the community can help you!
Start Hunting!