问题简介:
已知ahki(一维数组,长度为n,此处n=200),aH和atau(一维数组,长度都是29),常数param1,param2,sum1,函数sumFun;
要求eps(一维数组长度为n,但已知eps的第一、二个元素eps[[1]]=1,eps[[2]]=0.9966353966646951`)。
eps的计算过程是这样的:从第三个元素开始,其值计算与该值前面的元素有关(所以程序复杂度就上来了),具体的,
eps[[3]]=param1*eps[[2]]-param2*eps[[1]]-sumFun[2,1]*(eps[[2]]-eps[[1]])/sum1;
eps[[4]]=param1*eps[[3]]-param2*eps[[2]]-sumFun[3,1]*(eps[[2]]-eps[[1]])/sum1-sumFun[3,2]*(eps[[3]]-eps[[2]])/sum1;...
下面是程序:
akhi = {0, 2.3287936864188195`*^-6, ...2947`,
0.0015648884595468354`, 0.001582400475354751`};见https://mmaqa.com/?qa=blob&qa_blobid=2628598349081426233
aH = {0.001430085925000597`, 0.0016675838954522528`,
0.0019445237518651497`, 0.0022674557075536023`,
0.002644017786250091`, 0.003083116477520201`,
0.0035951373940762837`, 0.004192190913477146`,
0.004888398614142006`, 0.00570022728041296`, 0.006646878376194832`,
0.007750742197153878`, 0.009037927434105591`,
0.010538878629833178`, 0.012289095928226376`,
0.014329973160626298`, 0.016709769425524348`,
0.019484714892418913`, 0.02272017089502634`, 0.0264913657698098`,
0.030881308242713845`, 0.03596425066700329`, 0.04171792332650145`,
0.047571148303899205`, 0.04982598082764393`, 0.031665129676599436`,
0.008140592433938151`, 0.0018935256812937901`,
0.00045598357538307243`};
atau = {1.`*^-20, 1.`*^-19, 1.`*^-18, 1.`*^-17, 1.`*^-16, 1.`*^-15,
1.`*^-14, 1.`*^-13, 1.`*^-12, 1.`*^-11, 1.`*^-10, 1.`*^-9, 1.`*^-8,
1.`*^-7, 1.`*^-6, 0.00001`, 0.0001`, 0.001`, 0.01`, 0.1`, 1.`,
10.`, 100.`, 1000.`, 10000.`, 100000.`, 1.`*^6, 1.`*^7, 1.`*^8};
sumFun[n_,
j_] := (4*
Total[aH[[#]]*
Exp[-(akhi[[n + 1]] - akhi[[j + 1]])/atau[[#]]] & /@
Range[29]] -
2 aH[[1]]*Exp[-(akhi[[n + 1]] - akhi[[j + 1]])/atau[[1]]] -
2 aH[[29]]*
Exp[-(akhi[[n + 1]] - akhi[[j + 1]])/atau[[29]]])*2.302*0.125/4;
sum1 = 0.061087175858476464`;
param1 = 0.9697185699822564`;
param2 = -0.026916826682438735`;
解1.Nest (用了链表结构)
eps0 = Nest[{#,
param1*(res = Flatten[#])[[-1]] - param2*res[[-2]] -
Sum[sumFun[Length[res], i]*Differences[res][[i]]/sum1, {i, 1,
Length[res] - 1}]} &, {1, 0.9966353966646951`}, 200] //
Flatten; // AbsoluteTiming
返回{6.128351, Null}
解2.Do
eps1 = Flatten[{1, 0.9966353966646951`, Table[0, {200}]}];
Do[eps1[[n + 1]] =
param1*eps1[[n]] - param2*eps1[[n - 1]] -
Sum[sumFun[n, i]*(eps1[[i + 1]] - eps1[[i]])/sum1, {i, 1,
n - 1}], {n, 2, 201}]; // AbsoluteTiming
Norm[eps1 - eps0]
返回
{5.870336, Null}
0.
这个程序应该是规模n(此处n=200)的平方次复杂度(规模为n时,约算n^2/2次sumFun函数),我要算的规模n=30000,如果这样算肯定很久。
问题1. Do语句能不能优化?比如能不能Compile?(Compile我不懂,问错了的话见谅),我不会用,试了下,列在下面,结果不对(Do语句没有执行,结果还是返回原定义的eps2),怎么体验到原生的汇编速度。
eps2 = Flatten[{1, 0.9966353966646951`, Table[0, {200}]}];
Compile[{{n, _Integer}, {i, _Integer}},
Do[eps1[[n + 1]] =
param1*eps1[[n]] - param2*eps1[[n - 1]] -
Sum[sumFun[n, i]*(eps1[[i + 1]] - eps1[[i]])/sum1, {i, 1,
n - 1}], {n, 2, 201}]]; // AbsoluteTiming
Norm[eps0 - eps2]
Norm[Flatten[{1, 0.9966353966646951`, Table[0, {200}]}] - eps2]
返回
{0., Null}
14.1098
0.
问题2. 因为要多次调用sumFun函数,应该最耗时的就是n^2/2次sumFun函数的计算了吧,这个能怎么优化么?整个程序任何位置可以改进的都可以,谢谢大家!