# Singularity or stiff system suspected in ODE that shouldn't have a singularity in the solution range

I am trying to numerically solve the following ODE where I know the value at p[1000]:

diffEq = 2.2*10^9 + 7.6 p[v] - 4.*10^-8 p[v]^2 +  8.*10^-17 p[v]^3 == -v p'[v]
bc = p[1000] == 101.3 * 10^3;

NDSolve[{diffEq , bc} , p , {v, 5 , 1000}]


But I get the message

NDSolve::ndsz: At v == 815.3628963499186, step size is effectively zero; singularity or stiff system suspected.


If we put diffEq in normal form, we see there is indeed a singularity, but it shouldn't happen until v goes to 0 and I'm only solving in the region {v, 5 , 1000}:

Solve[diffEq, p'[v]]

{{p'[v] -> (-2.2*10^9 - 7.6 p[v] +
4.*10^-8 p[v]^2 - 8.*10^-17 p[v]^3)/v}}


v is the only thing in the denominator when diffEq is in normal form, so I'm confused as to what else could be causing the singularity/stiffness.

Some things I've tried so far:

I did find this related answer but the shooting method doesn't work here (though I may be misusing it, I can't say I understand completely what this does although I am making some progress understanding it here):

NDSolve[{diffEq, bc}, p, {v, 5, 1000},
Method -> {"Shooting", "StartingInitialConditions" -> {bc}}]

(*NDSolve::ndsz: At v == 815.3628963499186, step size is effectively zero; singularity or stiff system suspected.*)


And the suggestion of using ParametricNDSolve doesn't seem to change anything either:

sol = ParametricNDSolve[diffEq && p[1000] == p0, p, {v, 5, 1000}, p0]
bcs = Quiet[FindRoot[(p[p0][1000] == 101.3*10^3) /. sol, {p0, 100^3}],
ParametricNDSolve::ndsz];
p[p0] /. bcs /. sol
(*returns an interpolating function with Domain {{815., 1000.}}, but I still can't go below v = 815. *)

• Your ode can be analysed analytically. It is an autonomous equation, therefore, the first step is to change dependent and independent variables and to analyse the ode for $v(p)$. Commented May 22 at 8:37
• @yarchik Oh I think I understand! But to make sure I do, do you mean I change my differential equation to look like $$\frac{1}{a+b p+c p^2+d p^3}=-\frac{v'(p)}{v}$$ And now solve for $v(p)$ ?
– ydd
Commented May 22 at 15:48

If we can solve $$\int_{v_0}^{v_s} {dv \over v} = \int_{p_0}^\infty {dp \over A + Bp + Cp^2 +Dp^3} \,,$$ we will obtain the value $$v_s$$ at which the IVP $${dp \over dv} = {A + Bp + Cp^2 +Dp^3 \over v} \,,\quad p(v_0)=p_0\,,$$ has singularity.

sepvars = Times @@@ Apply[Power,
Lookup[
GroupBy[p'[v] /. First@Solve[diffEq, p'[v]] // FactorList,
FreeQ[p]],
{True, False}
],
{2}]
singEQ = Integrate[First[%], {v, 1000, vs}, Assumptions -> 0 < vs < 1000] ==
NIntegrate[1/Last[%], {p[v], 101.3*10^3, Infinity},
MaxRecursion -> 20]
Solve[singEQ]
(*
{1./v, -2.2*10^9 - 7.6 p[v] + 4.*10^-8 p[v]^2 - 8.*10^-17 p[v]^3}
-6.90776 + 1. Log[vs] == -0.204122

{{vs -> 815.3628806015575}}
*)


This is the approximate value of v where NDSolve quits:

NDSolveValue[{diffEq, bc}, p["Domain"], {v, 5, 1000}]
(*  {{815.3628963499186, 1000.}}   *)


We can get closer by increasing precision:

NDSolveValue[{diffEq, bc}, p["Domain"], {v, 5, 1000},
PrecisionGoal -> 12]
(*  {{815.3628806060308, 1000.}}  *)


Update: Numerical test for the order of the pole

The limit of $$(v-v_s)p'(v)/p(v)$$ is $$-k$$, where $$k$$ is the order of the pole at $$v=v_s$$. Here we see the order is $$k=1/2$$:

sol = NDSolveValue[{diffEq, bc}, p, {v, 5, 1000}, PrecisionGoal -> 12];
Plot[sol'[v] (v - sol["Domain"][[1, 1]])/sol[v],
{v,
First[sol@"Domain"] . {1 - 1*^-8, 1*^-8},
First[sol@"Domain"] . {0.99, 0.01}}]


For completeness:

$$k=0$$: Some other kind of singularity, like log, log-log, powers of log, etc. Example:

sol = NDSolveValue[{p'[v] == 1/v, p[1] == 1}, p, {v, 0, 1}]


$$-k>0$$: Weak singularity with $$p(v) \sim (v-v_s)^{-k}$$. Example:

sol = NDSolveValue[{p'[v] == 1/(3 p[v]^2), p[1] == 1}, p, {v, 0, 1}]
(*  -k = 1/3, p[v] = v^(1/3)  *)


This ODE also can be solved symbolically.

DSolve[Rationalize[diffEq, 0], p, v];
Rule @@ Reduce[%[[1]] /. v -> 1000 /. p[1000] -> Rationalize[101.3*10^3], C[1]];
SolveValues[%%[[1]] /. %, Log[v]][[1]] /. p[v] -> x

(* 312499999999999975000000 ((3 Log[2] + 3 Log[5] +
156249999999999987500000 RootSum[687499999999999945000000000000000 +
2374999999999999810000000 #1 - 12499999999999999 #1^2 +
25000000 #1^3 &, Log[101300 - #1]/(
1187499999999999905000000 - 12499999999999999 #1 +
37500000 #1^2) &])/312499999999999975000000 -
1/2 RootSum[687499999999999945000000000000000 +
2374999999999999810000000 #1 - 12499999999999999 #1^2 +
25000000 #1^3 &, Log[x - #1]/(
1187499999999999905000000 - 12499999999999999 #1 + 37500000 #1^2) &]) *)


For completeness, plot v, the exponential of this expression, against x to display the solution,

ps = Quiet@ParametricPlot[{Exp[%], x}, {x, 101300, 10^12},
ScalingFunctions -> {None, "Log"}, AspectRatio -> 1,
PlotRange -> All, PlotPoints -> 10000, PlotStyle -> {Dashed, Red}];


And compare it with the numerical solution.

NDSolveValue[{diffEq, bc}, p[v], {v, 5, 1000}];
pn = LogPlot[%, {v, 815.3628963499186, 1000}, PlotRange -> All];

Show[pn, ps, AxesLabel -> {v, p[v]}, LabelStyle -> {12, Bold, Black}]


If you integrate your equation backwards from 1000 you see that it will diverge near 815:

diffEq =
2.2*10^9 + 7.6  p[v] - 4.*10^-8  p[v]^2 +
8.*10^-17  p[v]^3 == -v  p'[v]
bc = p[1000] == 101.3*10^3;
x1 = 816; x2 = 1000;
sol[x_] = p[x] /. NDSolve[{diffEq, bc}, p, {v, x1, x2}][[1]]
Plot[sol[x], {x, x1, x2}]
`

• Thanks for answering. I am still scratching my head as to why it diverges near $v = 815$. If I start at $v =1000$ and plug in my boundary value, I will get a finite value for $p'(v)$. Since the derivative is just a cubic divided by v: $$p'(v)=-\frac{a+b p(v)+c p(v)^2+d p(v)^3}{v}$$ And $p(v)$ starts as a finite number, $p'(v)$ should stay finite, right? So If I keep moving down from $v = 1000$ in small steps of $dv$, then $dv ~p'(v)$ stays finite, meaning that $p(v - dv)$ stays finite. There is no singularity in $p'(v)$ until $v = 0$ so I'm still confused why $p(v)$ explodes at 815?
– ydd
Commented May 22 at 15:10