求助,做一个复杂函数的数值积分,涉及到变上限问题
mathematica吧
全部回复
仅看楼主
level 1
皮皮cutest 楼主
代码运行后,不报错但就是一直运行不出结果,我自己的电脑跑了将近5小时仍无输出,调试了许久也无结果。是因为函数太复杂的缘故吗,降低积分精度可行吗或者又有其他的隐含问题待解决呢?
代码如下:
I10[r_?NumberQ, td_?NumberQ, T_?NumberQ, t_?NumberQ] :=
NIntegrate[
Exp[-R^2]*R^2*Sin[theta]*Exp[t], {phi, 0, 2*Pi}, {theta, 0,
Pi}, {R, 0, (2*T + td)/2 - t}];
I1[r_?NumberQ, td_?NumberQ, T_?NumberQ] :=
NIntegrate[I10[r, td, T, t], {t, T - r/2, (2*T + td)/2}];
I20[r_?NumberQ, td_?NumberQ, T_?NumberQ, t_?NumberQ] :=
NIntegrate[
Exp[-R^2 + r^2 + 2*r*R*Cos[theta]]*R^2*Sin[theta]*Exp[t], {phi, 0,
2*Pi}, {theta, 0, Pi}, {R, 0, (2*T - td)/2 - t}];
I2[r_?NumberQ, td_?NumberQ, T_?NumberQ] :=
NIntegrate[I20[r, td, T, t], {t, T - r/2, (2*T - td)/2}];
Cx[r_, td_, T_,
t_] := -(((2*T + td)/2 - t)^2 +
r^2 - ((2*T - td)/2 - t)^2)/(2*((2*T + td)/2)*r);
Cy[r_, td_, T_,
t_] := (((2*T - td)/2 - t)^2 +
r^2 - ((2*T + td)/2 - t)^2)/(2*((2*T - td)/2)*r);
I30[r_?NumberQ, td_?NumberQ, T_?NumberQ, t_?NumberQ] :=
NIntegrate[
Exp[-R^2]*R^2*Sin[theta]*Exp[t], {phi, 0, 2*Pi}, {theta,
Pi - ArcCos[Cx[r, td, T, t]], Pi}, {R, 0, (2*T + td)/2 - t}];
I3[r_?NumberQ, td_?NumberQ, T_?NumberQ] :=
NIntegrate[I30[r, td, T, t], {t, -Infinity, T - r/2}];
I40[r_?NumberQ, td_?NumberQ, T_?NumberQ, t_?NumberQ] :=
NIntegrate[
Exp[-R^2 + r^2 + 2*r*R*Cos[theta]]*R^2*Sin[theta]*Exp[t], {phi, 0,
2*Pi}, {theta, 0, Pi - ArcCos[Cx[r, td, T, t]]}, {R,
0, (2*T - td)/2 - t}];
I4[r_?NumberQ, td_?NumberQ, T_?NumberQ] :=
NIntegrate[I40[r, td, T, t], {t, -Infinity, T - r/2}];
I50[r_?NumberQ, td_?NumberQ, T_?NumberQ, t_?NumberQ, theta_?NumberQ] :=
NIntegrate[
Exp[-R^2]*R^2*Sin[theta]*Exp[t], {phi, 0, 2*Pi}, {R,
0, (((2*T - td)/2 - t)*Cx[r, td, T, t])/Cos[theta]}];
I51[r_?NumberQ, td_?NumberQ, T_?NumberQ, t_?NumberQ] :=
NIntegrate[
I50[r, td, T, t, theta], {theta, 0,
Pi - ArcCos[Cx[r, td, T, t]]}];
I5[r_?NumberQ, td_?NumberQ, T_?NumberQ] :=
NIntegrate[I51[r, td, T, t], {t, -Infinity, T - r/2}];
I60[r_?NumberQ, td_?NumberQ, T_?NumberQ, t_?NumberQ, theta_?NumberQ] :=
NIntegrate[
Exp[-R^2 + r^2 + 2*r*R*Cos[theta]]*R^2*Sin[theta]*Exp[t], {phi, 0,
2*Pi}, {R, 0, -(((2*T + td)/2 - t)*Cy[r, td, T, t])/Cos[theta]}];
I61[r_?NumberQ, td_?NumberQ, T_?NumberQ, t_?NumberQ] :=
NIntegrate[
I60[r, td, T, t, theta], {theta, Pi - ArcCos[Cy[r, td, T, t]], Pi}];
I6[r_?NumberQ, td_?NumberQ, T_?NumberQ] :=
NIntegrate[I61[r, td, T, t], {t, -Infinity, T - r/2}];
P[r_, td_, T_] :=
Exp[-(I1[r, td, T] + I2[r, td, T] + I3[r, td, T] + I4[r, td, T] +
I5[r, td, T] + I6[r, td, T])];
Gx[r_, td_] := (E^((1/4)*(-(2*r*td) - r*(r + 2) - td^2))*(r -
td)*(r + td)*
(Sqrt[Pi]*
E^((1/4)*(r + td + 1)^2)*(r^4 - r^2*(td*(2*td + 3) + 3) +
td*(td + 1)*(td*(td + 2) + 7))*Erfc[(1/2)*(r + td + 1)] -
2*(r - td)*((r - 1)*r - td*(td + 2)) - 10*td));
Gy[r_, td_] := (E^((1/4)*((-r)*(2 + r) - 2*r*td - td^2))*(r - td)*(r +
td)*(2*(-r^3 + r^2*(1 + td) - r*td*(1 + td) +
td*(5 + td*(2 + td))) +
E^((1/4)*(1 + r + td)^2)*
Sqrt[Pi]*(r^4 - r^2*(3 + td) - td*(1 + td)*(7 + td*(2 + td)))*
Erfc[(1/2)*(1 + r + td)]));
f[k_?NumberQ, td_?NumberQ, T_?NumberQ] :=
NIntegrate[
Exp[2*T]*Cos[k*td]*P[r, td, T]*SphericalBesselJ[2, k*r]*Gx[r, td]*
Gy[r, td]*k/(r^6), {r, td, Infinity}];
integratedDelta[k_] :=
NIntegrate[f[k, td, T], {td, 0, Infinity}, {T, -Infinity, Infinity},
Method -> "MonteCarlo"];
integratedDelta[1]
2025年03月19日 14点03分 1
level 1
皮皮cutest 楼主
实际上 我还需要画一个代码中integratedelta关于k的一个曲线图 预计需要算很多次k不同取值,但现在一次也算不出来,很迷茫
2025年03月19日 14点03分 2
level 1
皮皮cutest 楼主
按理说Nintegrate不应该有一个理论二分尝试次数吗 达到阈值应该会有反馈,但现在情况是运行很久也无任何信息更新。。。
2025年03月19日 14点03分 3
level 1
皮皮cutest 楼主
2025年03月19日 14点03分 4
level 8
你这是一个三重嵌套的数值积分,如果把定义域离散成N个点的话,系统相当于要做N^3次独立的NIntegrate才能得到最后的结果,耗时是可想而知的。
所以没有信息更新再正常不过了。
2025年03月20日 02点03分 5
谢谢您的回复! 那这样的话岂不是只能降低离散点个数N的值来加快进程吗?
2025年03月20日 02点03分
吧务
level 15
……大致看了看,不止是多重嵌套,貌似还是无限域还振荡,够呛,就算能算出来,算不算得准都是问题。(蒙德卡洛法大概率是要跪。)关于这个可以参看SE帖子《Numerical inverse Laplace-Hankel transform》(编号:113240),这才两重,就已经算得人欲仙欲死了。总之,这类问题硬算积分并不可取,建议结合背景想想别的法子。
2025年03月28日 18点03分 6
1