# Message Boards

0
|
2164 Views
|
0 Replies
|
0 Total Likes
View groups...
Share
 Actually, I want do numerical integration, I use the Simpson method, the integration in on k and w, so there are two loops. I use Table to replace the for loop, but the speed almost the same. I also write the Matlab code, it's much fast. this result is right, I want to speed up it , how can I do? 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[2 I kz3 d]]; etemp = Exp[-2 Im[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/(4 Pi^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/(4 Pi^2)}];}, {ind, 1, n + 1}]; swtot = (swkx[[1]] + 4 Sum[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}] + 2 Sum[sw[[k]], {k, 3, nw - 1, 2}] + sw[[nw]]) dw/3; AppendTo[Qtd, Qt]; Print[DateString[]]; ]  Attachments: