24小时热门版块排行榜    

查看: 663  |  回复: 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 的主题更新
最具人气热帖推荐 [查看全部] 作者 回/看 最后发表
[基金申请] 国自然评审意见 +10 wangmingqi 2026-08-28 12/600 2026-08-28 22:06 by water_still
[基金申请] 投票:  有多少人是今天查系统知道结果的? +16 爱看书的可乐 2026-08-26 18/900 2026-08-28 17:09 by winsaint
[基金申请] 为什么资助数各大高校都创新高,自己申请怎么就这么难 +11 Kittylucky 2026-08-27 12/600 2026-08-28 16:26 by lyxiaomc
[基金申请] 有没有仍没收到信息的 +3 德尚中行 2026-08-27 3/150 2026-08-28 15:24 by rongshuxia
[基金申请] 基金不中,共勉 +11 eulota 2026-08-26 11/550 2026-08-28 14:22 by 火星超人xi
[基金申请] 面上合作单位盖章 +5 ssyjh 2026-08-27 5/250 2026-08-27 20:50 by gdfollow
[基金申请] 怎么看青基中了没有啊 +5 叶九微 2026-08-26 5/250 2026-08-27 10:35 by l_zh2008
[文学芳草园] 梦想 +3 myrtle 2026-08-26 3/150 2026-08-27 10:01 by angelyueyi
[基金申请] 为什么 国际(地区)合作与交流项目 没有放榜? 10+3 majunge000 2026-08-26 11/550 2026-08-27 08:42 by 北京莱茵编辑
[基金申请] 范进中举一文的中心思想 +9 炎黄贵胄 2026-08-22 10/500 2026-08-26 15:40 by semaglutide
[基金申请] 2026年8月25日国自然放榜前突然收到列入评审专家邮件,有关系吗? +25 木水思豆 2026-08-25 28/1400 2026-08-26 14:53 by draco1987
[基金申请] 出来了 +9 trojank 2026-08-26 9/450 2026-08-26 14:25 by 宝贝虫子
[基金申请] 系统进不去 +4 yanglien 2026-08-26 5/250 2026-08-26 11:10 by wenfengw83
[基金申请] 项目信息和经费信息在系统里都可以看到了 +6 wittyboy 2026-08-26 14/700 2026-08-26 10:55 by wittyboy
[基金申请] 国合可查了 +3 paperzjh 2026-08-26 3/150 2026-08-26 10:41 by LemmonTr
[基金申请] 国合现在查不到了吗? +10 chengyan1220 2026-08-24 20/1000 2026-08-26 08:57 by peasantsprig
[基金申请] 牛来!米来!面来! +8 beefly 2026-08-26 8/400 2026-08-26 08:37 by xuzhipiao
[基金申请] 在坚冰还盖着北海的时候,我看到了怒放的梅花。 (金币+10) +6 ziyangfang 2026-08-25 9/450 2026-08-25 20:26 by huagongfeihu
[基金申请] 明天应该可查了!? +6 chengyan1220 2026-08-23 6/300 2026-08-25 19:45 by zfd97
[基金申请] 没有任何消息-是不是就凉了 +9 图啦图啦 2026-08-24 10/500 2026-08-25 11:59 by 南海小哥
信息提示
请填处理意见