I'm trying to solve an ODE system with the NDSolve method. This my ODE system, with the BC and the functions V and dV defined:

ODE system, with BC in (4) and the function V(B) and dV(B) defined in (5)

    s = NDSolve[{F'[
  t] == -(2 /Sqrt[3]) (F[t]^2/2 + a[t]^2*  V[B[ t]])^(1/2)*F[t] - 
a'[t] == a[t]/Sqrt[3] (F[t]^2/2 + a[t]^2*  V[B[ t]])^(1/2), 
B'[t] == F[t], F[tmin] == Subscript[v, 0], a[tmin] == a0 , 
B[tmin] == Subscript[B, 0] }, {F, a, B}, {t, tmin, tmax}, Method->"ExplicitRungeKutta"]

Where the solution functions are:

F(t), a(t), B(t)

I need the solutions for a big interval, like:

t = [-10^(10), -10^(-10)]


tmin = -10^(10)
tmax = -10^(-10)

But the system seems to be stiff. I'm trying to replace the "ExplicitRungeKutta" method for "StiffnessSwitching method", replacing the method with:

Method -> {"StiffnessSwitching", Method -> {"ExplicitRungeKutta", Automatic}}

But the system seems to remain stiff.

The complete code is:


V[B_] := 3;
dV[B_] :=  V'[B];

B0 = 1;
F0 = 0;
a0 = 10^(-10);

tmin = -10^(10);
tmax = -10^(-10);

  s = NDSolve[{F'[
  t] == -(2 /Sqrt[3]) (F[t]^2/2 + a[t]^2*  V[B[ t]])^(1/2)*F[t] - 
a'[t] == a[t]/Sqrt[3] (F[t]^2/2 + a[t]^2*  V[B[ t]])^(1/2), 
B'[t] == F[t], F[tmin] == F0, a[tmin] == a0 , B[tmin] == B0 }, {F,
 a, B}, {t, tmin, tmax},  Method->"ExplicitRungeKutta"]

Question: Is there any way to obtain the solutions of these equations without stiffness? What method should I use?

I'd appreciate much all the answers.

  • $\begingroup$ What if you just use Method -> "StiffnessSwitching" without invoking RK at all? $\endgroup$
    – J. M.'s torpor
    Jun 21 '16 at 23:03
  • 2
    $\begingroup$ (1) What's v0? (2) What if it's a singularity and not stiffness? $\endgroup$
    – Michael E2
    Jun 22 '16 at 1:02
  • $\begingroup$ Are tmin & tmax switched, or is tmin supposed to be greater than tmax? $\endgroup$
    – Michael E2
    Jun 22 '16 at 3:09
  • $\begingroup$ V0 was F0, sorry for the mistake. And both variables tmin and tmax were negative (another mistake writting the code here). $\endgroup$ Jun 22 '16 at 23:57
  • $\begingroup$ And the stiffnessSwitching method alone does not solve the (aparent) stiffness problem. $\endgroup$ Jun 22 '16 at 23:59

Let's read the error message:

NDSolve`Iterate::ndsz: At t == -1.*10^10, step size is effectively zero; singularity or stiff system suspected.

The "step size is effectively zero" means that in floating-point arithmetic t + dt is equal to t for the computed time step dt. If this is the problem, increasing working precision might help.

de = {F'[t] == -(2/Sqrt[3]) (F[t]^2/2 + a[t]^2*V[B[t]])^(1/2)*F[t] - 
   a'[t] == a[t]/Sqrt[3] (F[t]^2/2 + a[t]^2*V[B[t]])^(1/2), 
   B'[t] == F[t]};
ics = {F[tmin] == F0, a[tmin] == a0, B[tmin] == B0};
s = 
 NDSolve[{de, ics}, {F, a, B}, {t, tmin, tmax},
  PrecisionGoal -> 8, AccuracyGoal -> 8, WorkingPrecision -> 20,
  Method -> "StiffnessSwitching"]

Mathematica graphics

Note: The first of the OP's differential equations has

F'[t] == <..stuff..> * F[t]

because dV[B] is zero and the second terms goes away. Since F[t] starts out at F[tmin] == 0, that means F[t] is a constant 0 (zero). Thus the third differential equation B'[t] == F[t] and its initial condition imply B[t] is a constant 1 (one).


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.