24小时热门版块排行榜    

查看: 519  |  回复: 1

firstlaker

金虫 (正式写手)

[求助] 请问这个程序如何优化呀,求数值积分的问题

如题。
unispec[w_,n_,d_,min_,max_,epsi10_,epsi1e_,epsi20_,epsi2e_]:=Module[{},
dkx=(max-min)/n;
kx=Table[0,{i,1,n+1}];
swkx=Table[0,{i,1,n+1}];
k0=w/c0;
Table[{
kx[[ind]]=min+(ind-1)*dkx;
kz10=If[Im[Sqrt[epsi10 k0^2-kx[[ind]]^2]]>=0,Sqrt[epsi10 k0^2-kx[[ind]]^2],-Sqrt[epsi10 k0^2-kx[[ind]]^2]];
kz20=If[Im[Sqrt[epsi20 k0^2-kx[[ind]]^2]]>=0,Sqrt[epsi20 k0^2-kx[[ind]]^2],-Sqrt[epsi20 k0^2-kx[[ind]]^2]];
kz1e=If[Im[Sqrt[epsi10 k0^2-kx[[ind]]^2 epsi10/epsi1e]]>=0,Sqrt[epsi10 k0^2-kx[[ind]]^2 epsi10/epsi1e],
-Sqrt[epsi10 k0^2-kx[[ind]]^2 epsi10/epsi1e]];
kz2e=If[Im[Sqrt[epsi20 k0^2-kx[[ind]]^2 epsi20/epsi2e]]>=0,Sqrt[epsi20 k0^2-kx[[ind]]^2 epsi20/epsi2e],
-Sqrt[epsi20 k0^2-kx[[ind]]^2 epsi20/epsi2e]];
kz3=Sqrt[k0^2-kx[[ind]]^2];
rs31=(kz3-kz10)/(kz3+kz10);
rp31=(epsi10 kz3-kz1e)/(epsi10 kz3+kz1e);
rs32=(kz3-kz20)/(kz3+kz20);
rp32=(epsi20 kz3-kz2e)/(epsi20 kz3+kz2e);
tpstemp=Abs[1-rs31 rs32 Exp[2 I kz3 d]];
tpptemp=Abs[1-rp31 rp32 Exp[2I kz3 d]];
etemp=Exp[-2Im[kz3]d];
If[kx[[ind]]<=k0,{trans=(1-Abs[rs31]^2)(1-Abs[rs32]^2)/tpstemp^2+(1-Abs[rp31]^2)(1-Abs[rp32]^2)/tpptemp^2;
swkx[[ind]]=kx[[ind]]trans/(4Pi^2);},
{trans=4 etemp Im[rs31]Im[rs32]/Abs[1-rs31 rs32 etemp]^2+4 etemp Im[rp31]Im[rp32]/Abs[1-rp31 rp32 etemp]^2;
swkx[[ind]]=kx[[ind]] trans/(4Pi^2)}];}
,{ind,1,n+1}];
swtot=(swkx[[1]]+4Sum[swkx[[k]],{k,2,n,2}]+2 Sum[swkx[[k]],{k,3,n-1,2}]+swkx[[n+1]])*dkx/3]

Qtd={};
For[i=1,i<=3,i=i+0.2,
d=10^i nm;
Print[10^i];
Print[DateString[]];
Table[{
wi=wstart+(ii-1)dw;
a1=(wi-10)/c0;
ae1=(wi+10)/c0;
ae2=Max[0.6/d,const1 wi/c0];
ae3=Max[5/d,const2 wi/c0];
Q1[[ii]]=h wi/(Exp[h wi/(kb T1)]-1);
Q2[[ii]]=h wi/(Exp[h wi/(kb T2)]-1);
twp[[ii]]=unispec[wi,nk,d,a0,a1,epsi01[wi],epsie1[wi],epsi02[wi],epsie2[wi]];
twe1[[ii]]=unispec[wi,nk,d,ae1,ae2,epsi01[wi],epsie1[wi],epsi02[wi],epsie2[wi]];
twe2[[ii]]=unispec[wi,nk,d,ae2,ae3,epsi01[wi],epsie1[wi],epsi02[wi],epsie2[wi]];
sw[[ii]]=(twp[[ii]]+twe1[[ii]]+twe2[[ii]])(Q1[[ii]]-Q2[[ii]])}
,{ii,1,nw}];
Qt=(sw[[1]]+4 Sum[sw[[k]],{k,2,nw-1,2}]+2Sum[sw[[k]],{k,3,nw-1,2}]+sw[[nw]])dw/3;
AppendTo[Qtd,Qt];
Print[DateString[]];
]
如上所示,其实整个程序就是一个双重的数值积分,所以这里需要有两个长循环,利用Simpson法求积分(两个长的Table语句,一个在子函数unispec里,一个在最后的结果里)。写的时候思路是做For循环,写完For换成Table而已.
我觉得这个运行得太慢了(结果应该没问题)。程序应该从哪里进行优化?
这里我尝试把Table换成 ParallelTable,算出来结果却不对了(都是0)。
回复此楼
已阅   回复此楼   关注TA 给TA发消息 送TA红花 TA的回帖

virtualzx

木虫 (著名写手)

直接用内置数值积分函数。比你自己写的快,精度还更高

发自小木虫IOS客户端
2楼2016-02-25 00:04:55
已阅   回复此楼   关注TA 给TA发消息 送TA红花 TA的回帖
相关版块跳转 我要订阅楼主 firstlaker 的主题更新
信息提示
请填处理意见