5
$\begingroup$

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. *)
$\endgroup$
2
  • 2
    $\begingroup$ 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)$. $\endgroup$
    – yarchik
    Commented May 22 at 8:37
  • $\begingroup$ @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)$ ? $\endgroup$
    – ydd
    Commented May 22 at 15:48

3 Answers 3

4
$\begingroup$

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)  *)
$\endgroup$
4
$\begingroup$

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}]

enter image description here

$\endgroup$
3
$\begingroup$

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}]

enter image description here

$\endgroup$
1
  • $\begingroup$ 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? $\endgroup$
    – ydd
    Commented May 22 at 15:10

Your Answer

By clicking “Post Your Answer”, you agree to our terms of service and acknowledge you have read our privacy policy.

Not the answer you're looking for? Browse other questions tagged or ask your own question.