Error in Double Integration problem
3 visualizzazioni (ultimi 30 giorni)
Mostra commenti meno recenti
%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 Commenti
Risposte (2)
Torsten
il 23 Apr 2023
Spostato: Torsten
il 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 Commenti
Torsten
il 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
il 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 Commenti
Torsten
il 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
il 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
![](https://www.mathworks.com/matlabcentral/answers/uploaded_files/1364578/image.png)
rather than
![](https://www.mathworks.com/matlabcentral/answers/uploaded_files/1364583/image.png)
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.
Vedere anche
Categorie
Scopri di più su Numbers and Precision in Help Center e File Exchange
Community Treasure Hunt
Find the treasures in MATLAB Central and discover how the community can help you!
Start Hunting!