3
$\begingroup$

I tried to solve a system first order differential equations together with a constraint equation. I use the Method -> Projection to check if the constraint holds at each step. However, two problems occur:

The first problem is that I get the messages:

General::ovfl: Overflow occurred in computation. and General::unfl: Underflow occurred in computation,

I do not know how to handle these messages.

The second problem is that I get the message:

NDSolve::nlnum: "The function value {Overflow[], Overflow[], Overflow[], Overflow[]} is not a list of numbers with dimensions {4} at {t, u[t], x[t], y[t], z[t]} ...

I have attached the image of the code with the given explanations in (**) for convenience.

I would appreciate it if anyone would give me some help.

Needs["DifferentialEquations`NDSolveProblems`"];
Needs["DifferentialEquations`NDSolveUtilities`"];
Needs["DifferentialEquations`InterpolatingFunctionAnatomy`"];

(Defining constant parameters)

θ = 0.182; w = 0; ρi = 4.2; xi = 0.0; ui = 8.71; adot = \
-6.4; zi = 26.57; om = 11.72; ti = 0.0; ts = 10000;

(Defining an invariant quantity that must hold true during the \ dynamical evolution of the system)

 ina = -Exp[-3 x[t]]/(
 2 (2 om + 3)) (om/6*y[t]^2 - u[t]^2*z[t]^2 + u[t]*y[t]*z[t]) + 
 u[t]*Exp[3 x[t]] (-16*π*ρi*Exp[3 (1 + w) (xi - x[t])])

(Defining the system of first order differential equations to be \ solved numerically)

  eqs = {x'[t] + 
  Exp[-3 x[t]]/(
  2 (2 om + 3)*
  u[t]) (om/3*y[t] + 
   u[t]*z[t] + θ (y[t] - 2 u[t]*z[t]) z[t]) - (θ*
   Exp[3 x[t]]/u[t]) (-16 π*ρi*
   Exp[3 (1 + w) (xi - x[t])]) == 0, 
 u'[t] + Exp[-3 x[t]]/(2 (2 om + 3)) (y[t] - 2 u[t]*z[t]) + 
 6 θ*
 Exp[3 x[t]]*(-8 π (1 - w)*ρi*
   Exp[3 (1 + w) (xi - x[t])]) == 0, 
 y'[t] - Exp[
  3 x[t]]*(48 π*ρi*(1 - w)*Exp[3 (1 + w)*(xi - x[t])]) ==
 0, z'[t] - 
 Exp[-3 x[t]]/(2 (2 om + 3) u[t]) (y[t] - 2 u[t]*z[t])*z[t] + 
 Exp[3 x[t]]/u[t] (-16 π*ρi*Exp[3 (1 + w)*(xi - x[t])]) ==
 0}

(*Defing initial conditions: *)

ics = {u[ti] == ui, x[ti] == xi, z[ti] == zi, 
y[ti] == -(1/(3 zi*θ + om))
  3 (6 adot* E^(2 xi)*ui + ui*zi - 2 ui*zi^2*θ + 
  96 E^(6 xi)* π *θ *ρi + 
  4 adot* E^(2 xi) *ui *om + 
  64 E^(6 xi)* π*θ *ρi* om)}

(*Defining variables to be evaluated: *)

 vars = {x, x', x'', u, u', y, z}

(*Using the NDSolve Command: *)

g = NDSolve[{eqs, ics}, vars, {t, ti, ts}, 
Method -> "StiffnessSwitching", 
Method -> {"Projection", Method -> "ExplicitRungeKutta", 
"Invariants" -> ina == 0}, InterpolationOrder -> Automatic, 
WorkingPrecision -> MachinePrecision]

This image is a print screen of the code together with codes for convenience

$\endgroup$
4
  • 1
    $\begingroup$ try integrating over a smaller time interval? $\endgroup$
    – chris
    Aug 22 '14 at 14:42
  • $\begingroup$ Yes the idea of decreasing the time of evaluation worked. Thanks Chris $\endgroup$
    – user14750
    Aug 22 '14 at 16:23
  • $\begingroup$ Also sometimes increasing WorkingPrecision can help with overflows. $\endgroup$ Aug 22 '14 at 20:04
  • 2
    $\begingroup$ Basically, the solution develops a singularity. Normally, you get a NDSolve::ndsz singular or stiffness message, but in this case overflow happens before the problem is detected. The solution returned should be accepted, since the singularity is a feature of the IVP and not something that needs to be avoided. (WorkingPrecision -> 32 catches the singularity before overflow for me.) $\endgroup$
    – Michael E2
    May 13 '20 at 17:29
1
$\begingroup$

Excuse, for me, I am using Mathematica 12.0.0 at present, this code works brilliantly.

g = NDSolve[{eqs, ics}, vars, {t, ti, ts}, 
  Method -> "StiffnessSwitching", 
  Method -> {"Projection", Method -> "ExplicitRungeKutta", 
    "Invariants" -> ina == 0}, InterpolationOrder -> Automatic, 
  WorkingPrecision -> MachinePrecision]

Output

Plot[x[t] /. g, {t, 0, 29.2}]

Output

Plot[x'[t] /. g, {t, 0, 29.2}]

Output

Plot[x''[t] /. g, {t, 0, 29.2}]

Output

Plot[u[t] /. g, {t, 0, 29.2}]

Output

Plot[u'[t] /. g, {t, 0, 29.2}]

Output

Plot[y[t] /. g, {t, 0, 29.2}]

Output

Plot[z[t] /. g, {t, 0, 29.2}]

Output

In general, it is worth the openness not to use the option Workingpresicion and InterpolationOrder. Mathematic handles this internally better without in most cases with difficulties. StiffnessSwitching and RungeKutta are very good options to get convergence and a solution of value and regard. But leave open with order the RungeKutta method should be. Mathematica does the work internally optimal.

Many hope that equal step size does the work better, but that is not true. Mathematica's internal default stepsize control is most often better.

Here is the plot from the package

Needs["DifferentialEquations`NDSolveProblems`"];
Needs["DifferentialEquations`NDSolveUtilities`"];

StepDataPlot[g]

Output

This bigger step size in the middle and smaller at the end of the integration was used.

Most probable loading these two packages does the difference.

Another slight enhance may offer loading the package:

Needs["DifferentialEquations`InterpolatingFunctionAnatomy`"];

too.

$\endgroup$

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.