24小时热门版块排行榜    

查看: 664  |  回复: 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 的主题更新
最具人气热帖推荐 [查看全部] 作者 回/看 最后发表
[基金申请] 我就是申请一个面上项目而已,这评审意见是按照杰青的条件评的吧? +3 gouxfjh 2026-08-28 3/150 2026-08-29 01:01 by jurkat.1640
[基金申请] 投票:  有多少人是今天查系统知道结果的? +16 爱看书的可乐 2026-08-26 18/900 2026-08-29 00:28 by winsaint
[基金申请] 国自然评审意见 +11 wangmingqi 2026-08-28 15/750 2026-08-29 00:17 by fangyl2005
[基金申请] 为什么资助数各大高校都创新高,自己申请怎么就这么难 +12 Kittylucky 2026-08-27 13/650 2026-08-29 00:04 by superceng
[基金申请] 国自然面上复盘~欢迎讨论 (金币+15) +13 晴天加油 2026-08-26 14/700 2026-08-28 18:36 by huagongfeihu
[基金申请] 有没有仍没收到信息的 +3 德尚中行 2026-08-27 3/150 2026-08-28 15:24 by rongshuxia
[教师之家] 导师吐槽:我怎么摊上了这么个极品研究生! +8 苏东坡二世 2026-08-23 8/400 2026-08-28 14:41 by lkforward
[基金申请] 2026年叶企孙基金 +3 bud_bud 2026-08-27 6/300 2026-08-28 11:53 by bud_bud
[基金申请] 基金未中,这种答复是模板吗? +5 zhaosm1982 2026-08-27 6/300 2026-08-27 16:00 by lfy8008
[基金申请] 看板上这么多中的,有点像50人群里49个人都是骗子的那种感觉…… +5 a089 2026-08-26 6/300 2026-08-27 14:05 by jonewore
[基金申请] 范进中举一文的中心思想 +9 炎黄贵胄 2026-08-22 10/500 2026-08-26 15:40 by semaglutide
[基金申请] 能否退出参与的面上项目解除限项 +23 koalala 2026-08-24 26/1300 2026-08-26 14:29 by 宝贝虫子
[基金申请] 国合里面能看到了 +7 一怀馨秋 2026-08-26 7/350 2026-08-26 11:23 by zhaosm1982
[基金申请] 国际合作可查了,中了面上 (EPI+1)(金币+50) +18 Ldrop2023 2026-08-26 18/900 2026-08-26 11:15 by cmrandy
[基金申请] 项目信息和经费信息在系统里都可以看到了 +6 wittyboy 2026-08-26 14/700 2026-08-26 10:55 by wittyboy
[基金申请] 牛来!米来!面来! +8 beefly 2026-08-26 8/400 2026-08-26 08:37 by xuzhipiao
[基金申请] 某些机构,以效率低为荣,以效率低作为存在感 +9 yuleib84 2026-08-25 10/500 2026-08-25 17:14 by alexon
[基金申请] 如果此刻你正在为国基感到焦虑,不妨来听听这首《基金之外》 +8 scalable 2026-08-24 8/400 2026-08-25 12:52 by jnhyjjm
[基金申请] 没有任何消息-是不是就凉了 +9 图啦图啦 2026-08-24 10/500 2026-08-25 11:59 by 南海小哥
[教师之家] 跳槽后在研项目怎么办? +5 简单化xn 2026-08-22 10/500 2026-08-23 12:38 by 简单化xn
信息提示
请填处理意见