Integral problem numerical error

1 次查看（过去 30 天）
Yalcin
Yalcin 2014-2-26
Hi, I am trying to solve a very complex equation numerically. I have following problem. Does anyone solve my problem? Thanks in advance.. My codes and error are in following.
function Fig1_pc_Gauss_scint_vs_norm_source_size clc
global kpp; global ksi; global s1x; global s1y; global s2x; global s2y;
PrT=7; %Prandtl numbers for temperature PrS=700; %Prandtl numbers for salinity beta=0.72; %Obukhov-Corrsin constant Q=2.75;
dist=100; %distance
KT=1; %Thermal diffusivity KS=1; %salt diffusivity XT=1e-8; %dissipation rate of mean squared temperature lambda=0.55e-6; %wavelength w=-3; %relative strength of temperature and salinity eta=1e-3; %kolmogorov microscale length eps=1e-7;
alphas=0.03;
px=0; py=0; %receiver coordinate
teta=KT/KS;
A=2.6e-4; B=1.75e-4;
k=2*pi/lambda;
j=1;
for rx=0:0.001:0.05
ry=rx;
r=(rx*rx+ry*ry)^0.5;
rt(j)=r*100;
coeff1=-1*pi*k*k*dist*beta*A*A*XT*(w-1)*(w-1)*((eps)^(-1/3))/(w*w*(w*w*teta+1-w*(1+teta))); cf1=1/(lambda*lambda*dist*dist);
dr1=@(s1x,s2x,s1y,s2y,kpp,ksi)((w*w*teta*(((exp(-3*beta*(eta^(4/3))/(2*PrT)*(kpp.^(4/3))-(Q*beta*eta*eta/PrT)*kpp.*kpp)).*(kpp.^(-8/3))))+... (((exp(-3*beta*(eta^(4/3))*(PrT+PrS)/(4*PrT*PrS)*(kpp.^(4/3))-(Q*beta*eta*eta*(PrT+PrS)/(2*PrT*PrS))*kpp.*kpp)).*(kpp.^(-8/3))))-... w*(1+teta)*(((exp(-3*beta*(eta^(4/3))/(2*PrS)*(kpp.^(4/3))-(Q*beta*eta*eta/PrS)*kpp.*kpp)).*(kpp.^(-8/3))))+... w*w*teta*Q*(eta^(2/3))*(((exp(-3*beta*(eta^(4/3))/(2*PrT)*(kpp.^(4/3))-(Q*beta*eta*eta/PrT)*kpp.*kpp)).*(kpp.^(-2))))+... Q*(eta^(2/3))*(((exp(-3*beta*(eta^(4/3))/(2*PrS)*(kpp.^(4/3))-(Q*beta*eta*eta/PrS)*kpp.*kpp)).*(kpp.^(-2))))-... w*(1+teta)*Q*(eta^(2/3))*(((exp(-3*beta*(eta^(4/3))*(PrT+PrS)/(4*PrT*PrS)*(kpp.^(4/3))-(Q*beta*eta*eta*(PrT+PrS)/(2*PrT*PrS))*kpp.*kpp)).*(kpp.^(-2))))).*... ((1-besselj(0,(kpp.*(((ksi*px+(1-ksi).*s1x-(1-ksi).*s2x).*(ksi*px+(1-ksi).*s1x-(1-ksi).*s2x)+(ksi*py+(1-ksi).*s1y-(1-ksi).*s2y).*(ksi*py+(1-ksi).*s1y-(1-ksi).*s2y)).^0.5))))));
innerIntegral1=@(s1x,s1y,s2x,s2y)(integral2(@(kpp,ksi)dr1(s1x,s1y,s2x,s2y,kpp,ksi),0,inf,0,1));
fc=@(s1x,s1y,s2x,s2y)((exp((-s1x.*s1x-s1y.*s1y)/(2*alphas*alphas))).*(exp((-s2x.*s2x-s2y.*s2y)/(2*alphas*alphas))).*... (exp((1i)*pi*((s1x-px).*(s1x-px)+(s1y-py).*(s1y-py)-(s2x-(px+rx)).*(s2x-(px+rx))-(s2y-(py+ry)).*(s2y-(py+ry)))/(lambda*dist))).*... (exp(coeff1*(innerIntegral1))));
innerIntegral2 = @(s1x)(@(y)integral3(@(s1y,s2x,s2y)fc(y,s1y,s2x,s2y),-inf,inf,-inf,inf,-inf,inf));
fcr(j) = (integral(innerIntegral2,-inf,inf))*cf1
j=j+1;
end
hold on plot(rt,fcr1,'--k','LineWidth',2); hold off
xlabel( '\it\bfr\rm (cm)','FontSize',24,'FontSize',24,'FontName','Times New Roman'); set(gca,'FontSize',24,'FontName','Times New Roman'); ylabel('|<\itu\rm(\bfp\rm)\itu*\rm(\bfp\rm+\bfr\rm)>|','FontSize',24,'FontName','Times New Roman');set(gca,'FontSize',24,'FontName','Times New Roman');
Error using integralCalc/finalInputChecks (line 511) Input function must return 'double' or 'single' values. Found 'function_handle'.
Error in integralCalc/iterateScalarValued (line 315) finalInputChecks(x,fx);
Error in integralCalc/vadapt (line 133) [q,errbnd] = iterateScalarValued(u,tinterval,pathlen);
Error in integralCalc (line 104) [q,errbnd] = vadapt(@minusInfToInfInvTransform,interval);
Error in integral (line 89) Q = integralCalc(fun,a,b,opstruct);
Error in figure01_r_eksen_dongu_L_v1 (line 72) fcr(j) = (integral(innerIntegral2,-inf,inf))*cf1

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

回答（0 个）

类别

在 Help Center 和 File Exchange 中查找有关 MATLAB 的更多信息

标签

Community Treasure Hunt

Find the treasures in MATLAB Central and discover how the community can help you!

Start Hunting!

Translated by