代码运行后,不报错但就是一直运行不出结果,我自己的电脑跑了将近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]
代码如下:
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]