matlab数值积分计算核,matlab数值积分计算求指导数学小木虫学术科研互动社区...

用数值积分!

程序如下:

function jifen

re &＃61; 2.8e-15;

Rsun &＃61; 6.96e8;          % 太阳半径

f&＃61;32e9;                 % 信号工作频率

c &＃61; 3e8;

lamba &＃61; c/f;

lo&＃61;2e7;                 % 外尺度大小

p &＃61; 11/3;

AU&＃61;1.496e11;

Lse&＃61;AU;

Lsp&＃61;1.5*AU;

v &＃61; 300e3;

omega &＃61; 2.7e-6;

i&＃61;0;

for m &＃61; 3:2:30

i &＃61; i&＃43;1;

Rhelio&＃61;m*Rsun;

temp1&＃61;4*pi^3*re^2*lamba^2*lo^(3-p)*pi^((p-1)/2)*gamma((p-1)/2);

temp2&＃61;sin(pi*p/4)*gamma(1/2)*gamma(p/2)*gamma((p-3)/2);

L&＃61;Lse*cos(asin(Rhelio/Lse))&＃43;Lsp*cos(asin(Rhelio/Lsp));

c&＃61;[Lse,Rhelio,lamba,v,Rsun,p,omega,L];

temp3&＃61;rombg(c,L,1e-6);

chi2(i)&＃61;temp1/temp2*temp3;

end

chi2

function y&＃61;f(z,c)

Lse&＃61;c(1);

Rhelio&＃61;c(2);

lamba&＃61;c(3);

v&＃61;c(4);

Rsun&＃61;c(5);

p&＃61;c(6);

omega&＃61;c(7);

L&＃61;c(8);

R&＃61;sqrt((z-Lse*cos(asin(Rhelio/Lse))).^2&＃43;Rhelio^2);

Rf&＃61;sqrt(lamba.*z.*(L-z)/L);

Ar&＃61;160./((R/Rsun).^1.5)&＃43;1;

theta&＃61;atan(-v./(omega.*R));

y&＃61;R.^(-4).*Ar.*Rf.^(p-2)./(1&＃43;(Ar.^2-1).*(sin(theta)).^2).^(0.5);

function temp3&＃61;rombg(c,L,er)

h&＃61;L;

T(1,1)&＃61;h*(f(0,c)&＃43;f(L,c))/2;

i&＃61;2;h&＃61;h/2;

T(i,1)&＃61;T(1,1)/2&＃43;h*f(L/2,c);

T(i,2)&＃61;(4*T(2,1)-T(1,1))/3;

while abs(T(i,i)-T(i-1,i-1))>er

i&＃61;i&＃43;1;h&＃61;h/2;

T(i,1)&＃61;T(i-1,1)/2&＃43;h*sum(f(h:2*h:L,c));

for k&＃61;2:1:i

T(i,k)&＃61;(4^(k-1)*T(i,k-1)-T(i-1,k-1))/(4^(k-1)-1);

end

end

temp3&＃61;T(i,i);

运行结果:

>> chi2 &＃61;

1.0e-058 *

Columns 1 through 7

0.1517    0.1500    0.1474    0.1441    0.1401    0.1355    0.1304

Columns 8 through 14

0.1249    0.1192    0.1132    0.1072    0.1011    0.0951    0.0892

>>&＃xff0c;