24小时热门版块排行榜    

查看: 647  |  回复: 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 的主题更新
最具人气热帖推荐 [查看全部] 作者 回/看 最后发表
[基金申请] fileCode有新解读? +9 Tide man 2026-08-08 17/850 2026-08-08 14:56 by archvillain
[博后之家] 售SCI一区T0P文章,我:8O.55.1.O.54,科目全,可伽急 +3 HEQlVqMTIA7d 2026-08-07 6/300 2026-08-08 14:42 by oEVWOejN9taj
[公派出国] 售SCI文章,我:8O5.5.1.O.54,科目齐全,可+急 +3 HEQlVqMTIA7d 2026-08-07 6/300 2026-08-08 14:39 by oEVWOejN9taj
[教师之家] 售SCI一区T0P文章,我:8.O.55.1.O54,科目全,可伽急 +4 HEQlVqMTIA7d 2026-08-07 5/250 2026-08-08 14:27 by oEVWOejN9taj
[论文投稿] 售SCI一区T0P文章,我:8O.55.1.O.5.4,科目齐全,可+急 +3 HEQlVqMTIA7d 2026-08-07 5/250 2026-08-08 14:22 by oEVWOejN9taj
[基金申请] 售SCI-T0P文章,我:8O.5.5.1.O.54,科目齐全,可+急 +4 KXLV3nuBVBY7 2026-08-07 5/250 2026-08-08 14:07 by oEVWOejN9taj
[教师之家] 售SCI一区文章,我:8.O.55.1.O.54,科目齐全,可伽急 +3 KXLV3nuBVBY7 2026-08-07 5/250 2026-08-08 14:07 by oEVWOejN9taj
[基金申请] 国基金的申报应该改成非等额制,评价高的钱多评价低的钱少,但是增加资助率 +5 a089 2026-08-07 5/250 2026-08-08 14:04 by alian_214
[论文投稿] 售SCI文章,我:8O.5.5.1O.54,科目全,可十急 +3 KXLV3nuBVBY7 2026-08-07 6/300 2026-08-08 14:01 by oEVWOejN9taj
[基金申请] 关于代码变化问题,想知道的进来 +11 且听虎啸 2026-08-07 14/700 2026-08-08 12:50 by 布布和一二
[找工作] 售SCI-T0P文章,我:8O.5.5.1.O.54,科目齐全,可+急 +3 Vi50GxzrFcSG 2026-08-07 4/200 2026-08-08 10:39 by 3OOjAIS77qg2
[基金申请] 基金中了 +14 laoda193707 2026-08-06 14/700 2026-08-08 00:23 by 实验小白ha
[基金申请] 娱乐 +6 Tide man 2026-08-03 6/300 2026-08-07 22:40 by 铁帽子农民
[基金申请] 固定端突然变了,今天 +6 archvillain 2026-08-06 10/500 2026-08-07 16:03 by 医学老男孩
[基金申请] 听说今天filecode变了 +24 布布和一二 2026-08-06 47/2350 2026-08-07 16:02 by zhiyanjiang
[基金申请] filecode +14 等待解的谜 2026-08-06 19/950 2026-08-07 12:20 by wlwhappy
[基金申请] 影响面上的因素 +8 布布和一二 2026-08-05 11/550 2026-08-06 10:41 by 宝贝虫子
[基金申请] 有没有H口的?有收到消息的吗? +3 超级海虾 2026-08-04 3/150 2026-08-04 17:26 by 学教育滴
[基金申请] 面上提前没消息,有中的吗 +14 archvillain 2026-08-02 18/900 2026-08-04 14:42 by archvillain
[基金申请] 什么时候能放榜呀? +3 Jacob678 2026-08-03 3/150 2026-08-03 16:14 by gltch
信息提示
请填处理意见