0
$\begingroup$

I'm trying to numerically solve a system of differential equations which contains functions defined using NIntegrate:

Sa[y_] := (2*π*y)/(1 - Exp[-2 π*y]);
MSa[z_?NumericQ] := 2/√π NIntegrate[Sa[(z/u)^(1/2)]*u^(1/2)*Exp[-u], {u, 0, +∞}];

Sb[y_] := (2^10*π)/3*y^5/(1 + y^2)^2*Exp[-4 y*ArcCot[y]]/(1 - Exp[-2 π*y]);
MSb[z_?NumericQ] :=  2/√π NIntegrate[Sb[(z/u)^(1/2)]*(u^(1/2)*Exp[z])/(Exp[z + u] - 1)*Exp[-u], {u, 0, +∞}];

Yeq[x_?NumericQ] := c3 x^(3/2) Exp[-x];

fion[z_?NumericQ] :=  2^7/3 NIntegrate[η/(1 + η^2)^2*Exp[-4 η*ArcCot[η]]/(1 - Exp[-2 π*η])*1/(Exp[z (1 + 1/η^2)] - 1), {η, 0, +∞}];

odes := {Y'[z] == (-c1*MSa[z])/z^2 (Y[z]^2 - Yeq[4*z/α^2]^2) - (c1*MSb[z])/z^2 Y[z]^2 + c2*z*fion[z] (S[z] + T[z]), S'[z] == (c1*MSb[z])/4 z^2*Y[z]^2 - c2 z (1 + fion[z]) S[z], T'[z] == (3*c1*MSb[z])/(4*z^2)*Y[z]^2 - c2 z (ca + fion[z]) T[z], Y[(5*α^2)/4] == Yeq[5], S[(5*α^2)/4] == 1/2*Yeq[10 - (5*α^2)/4], T[(5*α^2)/4] == 3/2*Yeq[10 - (5*α^2)/4]}

sol = NDSolve[odes /. {m -> 1000, α -> 137^-1, Subscript[M, pl] -> 2.44*10^18}, {Y, S, T}, {z, (5*137^-2)/4, 10}]

where c1, c2, c3 and ca are fixed constants given by:

c1 := 2.13 Subscript[M, pl] m^-1 α^4
c2 := 0.003 Subscript[M, pl] α m^-1
c3 := 0.0026
ca := 4 (π^2 - 9) * α/(9*\[Pi])

But I get the error

NDSolve::ndsz: At z == 0.00006659917949810859`, step size is effectively zero; singularity or stiff system suspected. >>

Following the suggestion of @Michael E2 the integrations runs without errors, however the results is pretty crazy:

LogLogPlot[Y[z] /. sol, {z, 10^-3, 10}, PlotRange -> {{10^-3, 10}, {10^-15, 1}}]

enter image description here

$\endgroup$
  • 1
    $\begingroup$ You seem to gave us the wrong code. f[x] is not defined in your code but used in the NDSolve $\endgroup$ – Julien Kluge Sep 29 '16 at 14:40
  • $\begingroup$ You are right, I've just edit the post now should be correct. $\endgroup$ – Cervantes Sep 29 '16 at 14:48
  • $\begingroup$ Okay now a "c3" is undefined for the initial constants. $\endgroup$ – Julien Kluge Sep 29 '16 at 14:58
  • $\begingroup$ c3 it's just a constant. $\endgroup$ – Cervantes Sep 29 '16 at 15:14
  • $\begingroup$ Then give us examples for values for c1, c2, c3, ca please. $\endgroup$ – Julien Kluge Sep 29 '16 at 15:46
1
$\begingroup$

The relative sizes of the initial values for the derivatives and for the variables, suggests you need (at the least) a high WorkingPrecision, as well as "StiffnessSwitching", taking a hint from the OP's error.

IVs = Rule @@@ odes[[-3 ;;]];
params = {m -> 1000, α -> 137^-1, Subscript[M, pl] -> 2.44*10^18};
{Y'[z], S'[z], T'[z]} /. 
    First@Solve[odes[[;; 3]], {Y'[z], S'[z], T'[z]}] /. 
   z -> (5*α^2)/4 /. IVs /. params
(*  {1.76281*10^23, -6.61822*10^26, -1.98517*10^27}  *)

{Y[z], S[z], T[z]} /. z -> (5*α^2)/4 /. IVs /. params
(*  {0.000195865, 1.86648*10^-6, 5.59944*10^-6}  *)

The following integrates without error:

c1 := 2.13`100 * Subscript[M, pl] m α^4
c2 := 0.003`100 * Subscript[M, pl] α m^5
c3 := 0.0026`100;
ca := 4 (π^2 - 9)*α/(9*π)

sol = NDSolve[odes /. params, {Y, S, T}, {z, (5*137^-2)/4, 10}, 
   PrecisionGoal -> 8, AccuracyGoal -> 8, WorkingPrecision -> 32, 
   Method -> "StiffnessSwitching"];
$\endgroup$
  • $\begingroup$ Now it runs without errors, but still the results is pretty crazy (see the edit on the question). $\endgroup$ – Cervantes Sep 30 '16 at 8:37
  • $\begingroup$ @Cervantes Try LogLogPlot[Abs@Y[z] /. sol, {z, 10^-3, 10}, PlotRange -> {{10^-3, 10}, {10^-15, 1}}]: i.stack.imgur.com/MDAFf.png. The logarithm of a negative number is complex and won't be plotted. The solution is oscillatory, which appears as a sequence of humps, in case you think the humps are crazy. $\endgroup$ – Michael E2 Sep 30 '16 at 10:00

Your Answer

By clicking “Post Your Answer”, you agree to our terms of service, privacy policy and cookie policy

Not the answer you're looking for? Browse other questions tagged or ask your own question.