# Two problems with NDSolve when using Method -> Projection

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["DifferentialEquationsNDSolveProblems"];
Needs["DifferentialEquationsNDSolveUtilities"];
Needs["DifferentialEquationsInterpolatingFunctionAnatomy"];


(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] • try integrating over a smaller time interval? Aug 22 '14 at 14:42
• Yes the idea of decreasing the time of evaluation worked. Thanks Chris Aug 22 '14 at 16:23
• Also sometimes increasing WorkingPrecision can help with overflows. Aug 22 '14 at 20:04
• 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.) May 13 '20 at 17:29

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] Plot[x[t] /. g, {t, 0, 29.2}] Plot[x'[t] /. g, {t, 0, 29.2}] Plot[x''[t] /. g, {t, 0, 29.2}] Plot[u[t] /. g, {t, 0, 29.2}] Plot[u'[t] /. g, {t, 0, 29.2}] Plot[y[t] /. g, {t, 0, 29.2}] Plot[z[t] /. g, {t, 0, 29.2}] 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["DifferentialEquationsNDSolveProblems"];
Needs["DifferentialEquationsNDSolveUtilities"];

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

Needs["DifferentialEquationsInterpolatingFunctionAnatomy"];