0
$\begingroup$

My question is about a reason of the inability of Mathematica to integrate the given system of differential equations, that depends on parameters, in a stable way: for some values of the parameters the integration proceeds without issues, whereas for slightly different values it fails, giving nderr (error test failure). These values do not lead to singularity or other problems, so I don't understand. I would be very grateful for your help.

Consider the datasets Table2, Table3. The simple input with them, i.e., Table2={{...},...}, may be viewed here (I cannot insert it here since the question body symbol limit is exceeded). Excuse me for the inconvenience.

Let us then define the following interpolations:

RatenTopWeak1[t_] = 
  Interpolation[Table2[[All, {1, 2}]], InterpolationOrder -> 1][t];
RatepTonWeak1[t_] = 
  Interpolation[Table2[[All, {1, 3}]], InterpolationOrder -> 1][t];
RatenpTodg1[t_] = 
  Interpolation[Table2[[All, {1, 4}]], InterpolationOrder -> 1][t];
RatedgTonp1[t_] = 
  Interpolation[Table2[[All, {1, 5}]], InterpolationOrder -> 1][t];
RatedDiss1[t_, YN_, tauN_] = ( Exp[-t/tauN] YN)/
   tauN Interpolation[Table2[[All, {1, 6}]], InterpolationOrder -> 1][
    t];
RatepTonPi1[t_, YN_, tauN_] = ( Exp[-t/tauN] YN)/
   tauN Interpolation[Table2[[All, {1, 7}]], InterpolationOrder -> 1][
    t];
RatenTopPi1[t_, YN_, tauN_] = ( Exp[-t/tauN] YN)/tauN*
   Interpolation[Table2[[All, {1, 8}]], InterpolationOrder -> 1][t];
\[Eta]Bsm[t_] = 
  Interpolation[Table2[[All, {1, 9}]], InterpolationOrder -> 1][t];
tauNint[mN_] = Interpolation[Table3, InterpolationOrder -> 1][mN];
YNv[mN_, tauN_] = Min[0.01, 1.2*10^10 tauNint[mN]/tauN]

I then define the followins system of equations:

SystemDynamics[mN_, tauN_,tmax_] := Block[{mass = mN, lifetime = tauN},
  tstart = 0.1;
  DynamicsnToy[t_, Xn_, 
    Xd_] = -Xn*(RatenTopWeak1[t] + 
       RatenTopPi1[t, YNv[mass, lifetime], lifetime]) + (1 - Xn - 
       Xd) (RatepTonWeak1[t] + 
       RatepTonPi1[t, YNv[mass, lifetime], lifetime]) - 
    Xn (1 - Xn - Xd) RatenpTodg1[t] + 
    Xd*1/\[Eta]Bsm[t] RatedgTonp1[t] + 
    Xd*RatedDiss1[t, YNv[mass, lifetime], lifetime];
  DynamicsdToy[t_, Xn_, Xd_] = 
   Xn (1 - Xn - Xd) RatenpTodg1[t] - 
    Xd*1/\[Eta]Bsm[t] RatedgTonp1[t] - 
    Xd*RatedDiss1[t, YNv[mass, lifetime], lifetime];
  SolutionNuclearToy = 
   NDSolve[{D[Xn[t], t] == DynamicsnToy[t, Xn[t], Xd[t]], 
     D[Xd[t], t] == DynamicsdToy[t, Xn[t], Xd[t]], Xn[tstart] == 0.5, 
     Xd[tstart] == 0}, {Xn,(*Xp,*)Xd}, {t, tstart, tmax}]
  ]

Here, mN, tauN vary in limits $$ 0.05 < mN < 3, \quad 0.01 < tauN < 10^7, $$ and the preferable value of tmax is tmax = 10^9.

It turns out that it fails to be evaluated: when launching, say, SystemDynamics[1, 0.02, 10^9] I get the message

NDSolve::nderr: Error test failure at t == 1.236971642471891`; unable to continue. Other values of tauN lead to a similar error with a similar error point. Probably, the error is related to the fact that RatedgTonp1, RatenpTodg1 turn on suddenly at this time. 

After adding several options to NDSolve:

Method -> {"BDF", "MaxDifferenceOrder" -> 2}, PrecisionGoal -> 11, AccuracyGoal -> 12, InterpolationOrder -> 3, Compiled -> True, MaxSteps -> 10^6

the system started working, however, in a buggy way: some values of tauN work, whereas some don't. For example, SystemDynamics[1, 1.3, 10^9] and SystemDynamics[1, 1.5, 10^9] work, whereas SystemDynamics[1, 2, 10^9] does not, displaying the same nderr error.

Could you please tell me how to fix this?

$\endgroup$
1
  • 1
    $\begingroup$ Method -> StiffnessSwitching helps, one needs to carefully check whether the solution is reliable or not though. BTW, if you have to include the data, consider using the SE uploader. $\endgroup$
    – xzczd
    Sep 10 '20 at 2:56
2
$\begingroup$

This is not a complete solution but a fair enhancement of Yours. This seems to be solvable by the method BFD.

Now the Mathematica documentation is far from beginner-friendly. So have a look at NDSolve IDA Method. This states that BFD is part of IDA. IDA fits best to the given problem of Your question.

SystemDynamics[mN_, tauN_, tmax_, idaopts___] := 
 Module[{mass = mN, lifetime = tauN, time, steps}, tstart = 0.1;
  DynamicsnToy[t_, Xn_, 
    Xd_] = -Xn*(RatenTopWeak1[t] + 
       RatenTopPi1[t, YNv[mass, lifetime], lifetime]) + (1 - Xn - 
       Xd) (RatepTonWeak1[t] + 
       RatepTonPi1[t, YNv[mass, lifetime], lifetime]) - 
    Xn (1 - Xn - Xd) RatenpTodg1[t] + 
    Xd*1/\[Eta]Bsm[t] RatedgTonp1[t] + 
    Xd*RatedDiss1[t, YNv[mass, lifetime], lifetime];
  DynamicsdToy[t_, Xn_, Xd_] = 
   Xn (1 - Xn - Xd) RatenpTodg1[t] - 
    Xd*1/\[Eta]Bsm[t] RatedgTonp1[t] - 
    Xd*RatedDiss1[t, YNv[mass, lifetime], lifetime];
  time = First[
    Timing[SolutionNuclearToy = 
      NDSolve[{D[Xn[t], t] == DynamicsnToy[t, Xn[t], Xd[t]], 
        D[Xd[t], t] == DynamicsdToy[t, Xn[t], Xd[t]], 
        Xn[tstart] == 0.5, Xd[tstart] == 0}, {Xn,(*Xp,*)Xd}, {t, 
        tstart, tmax}, 
       Method -> {"MethodOfLines", 
         "DifferentiateBoundaryConditions" -> False, 
         Method -> {IDA, idaopts}}]]]; 
  steps = Length[
     First[
      InterpolatingFunctionCoordinates[
       First[u /. SolutionNuclearToy]]]] "Steps"; {time, steps, 
   SolutionNuclearToy}]

sol = Table[
  SystemDynamics[mNt, 0.01, 300, 
   "ImplicitSolver" -> {"Newton", 
     "LinearSolveMethod" -> "Band"}], {mNt, 0.051, 3, 10^9}]

gives solutions for all parameters. And all nearly the same:

output across the range of tauN

But already

sol = Table[
  SystemDynamics[mNt, 1, .1, 
   "ImplicitSolver" -> {"Newton", 
     "LinearSolveMethod" -> "Band"}], {mNt, 0.051, 3, 0.3}]

fails again.

Have a look at the solutions for those parameters:

enter image description here

The problem is not the long time integration, but a short time. It seems as is the very same problem for other methods that the boundary conditions are not properly matched for this data.

The patterns are quite similar. A typical solution overall time looks this way:

Xn and Xd

short time, long time behavior at the lower boundary of tauN

To me it seems that the high time upper boundary might not suit the complete parameter field.

sol = Table[
  SystemDynamics[mNt, 0.011, 10^9, "ImplicitSolver" -> "GMRES"], {mNt,
    0.051, 3, 0.3}]
Row[{Plot[Table[sol[[i, 3, 1, 1, 2]][t], {i, 10}], {t, 0, 10^9}, 
   ImageSize -> 300, PlotRange -> All], 
  Plot[Table[sol[[i, 3, 1, 2, 2]][t], {i, 10}], {t, 0, 10^9}, 
   ImageSize -> 300, PlotRange -> Full]}]

long integration time

For shorter integration times this happens:

Row[{Plot[Table[sol[[i, 3, 1, 1, 2]][t], {i, 10}], {t, 0, 10^4}, 
   ImageSize -> 300, PlotRange -> {All, {-0.02, 0.15}}], 
  Plot[Table[sol[[i, 3, 1, 2, 2]][t], {i, 10}], {t, 0, 10^4}, 
   ImageSize -> 300, PlotRange -> Full]}]

shorter integration times

That is happening in the longer one too, but the resolution of the plot does show it. It causes the problems for higher tauN values to diverges rapidly instead of reaching this steady state.

$\endgroup$
1
  • $\begingroup$ Thanks. So your guess is that the problem is for the initial conditions, and choosing them in the proper way will help to avoid the error. I will try to resolve this problem. $\endgroup$ Sep 10 '20 at 13:11

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.