14
$\begingroup$

Error when using NDSolve for $\epsilon y'' - y' + y^2 = 1$ with $0<x<1$ and $y(0) = \frac{1}{3}$, $y(1)=1$

My attempt:

Subscript[ϵ, 1] = 0.1;
Subscript[eqn, 1] = 
  NDSolve[{Subscript[ϵ, 1]*y''[x] - y'[x] + y[x]^2 == 1, 
    y[0] == 1/3, y[1] == 1}, y, {x, 0, 1}];

Error:

NDSolve::ndsz: At x == 0.9005611463004037`, step size is effectively zero; singularity or stiff system suspected.

However I see from this post https://math.stackexchange.com/a/3217456/697826 that it was numerically solved. What is wrong with my syntax?


Thank you to both of the submitted Answers. They are equally helpful!

$\endgroup$

3 Answers 3

18
$\begingroup$

Nothing wrong with your syntax. This is because

  1. The default setting of "Shooting" method doesn't work well on this problem.

  2. There's an arguable backslide of "Shooting" method around v11. See this and this post for more info.

If you're in v12.0 or higher, in which the nonlinear FEM solver is added, the simplest work-around is to turn to FiniteElement:

Subscript[ϵ, 1] = 0.1;
solfem = NDSolveValue[{Subscript[ϵ, 1]*y''[x] - y'[x] + y[x]^2 == 1, 
   y[0] == 1/3, y[1] == 1}, y, {x, 0, 1}, Method -> FiniteElement]

Plot[solfem[x], {x, 0, 1}]

enter image description here

If you're before v11, then the following should work:

Subscript[ϵ, 1] = 0.1;
solshoot = NDSolveValue[{Subscript[ϵ, 1]*y''[x] - y'[x] + y[x]^2 == 1, 
   y[0] == 1/3, y[1] == 1}, y, {x, 0, 1}, 
  Method -> {"Shooting", 
    "StartingInitialConditions" -> {y'[1] == 1, y[1] == 1}}]

solshoot // ListPlot

enter image description here

Notice I've made use of an undocumented syntax of ListPlot here. See this post for more info.

$\endgroup$
13
$\begingroup$

First, there are multiple solutions. Second, the one shown on math.SE has such a narrow range of feasible starting initial conditions fo $\epsilon = 1/10$ that the default shooting method in NDSolve is unlikely to find it, however the options are tweaked. It's probably best to manually implement the shooting method.

\[Epsilon] = 1/10;
psol = ParametricNDSolveValue[{\[Epsilon]*y''[x] - y'[x] + y[x]^2 == 
    1, y[0] == 1/3, y'[0] == yp}, y, {x, 0, 1}, {yp}]

ypsol = FindRoot[psol[yp][1] == 1, {yp, -0.94, -0.93}, 
  Method -> "Brent"]

(*  {yp -> -0.93371}  *)

sol = psol[yp] /. ypsol
ListLinePlot@sol

Here's another form of the shooting method, which allows one to figure out the interval to pass to the Brent method above:

Manipulate[
 ListLinePlot@
  NDSolveValue[{\[Epsilon]*y''[x] - y'[x] + y[x]^2 == 1, y[0] == 1/3, 
    y'[0] == yp}, y, {x, 0, 1}],
 {yp, -0.95, -0.5, Appearance -> "Labeled"}]
$\endgroup$
2
  • $\begingroup$ Thank you Michael! $\endgroup$
    – KZ-Spectra
    May 3, 2021 at 4:17
  • $\begingroup$ You're welcome, @Nalt! :) $\endgroup$
    – Michael E2
    May 3, 2021 at 12:06
3
$\begingroup$

Things are easy to solve, if you differentiate one time.

e1 = 1/10;
eqs0 = {e1*y''[x] - y'[x] + y[x]^2 == 1, y[0] == 1/3, y[1] == 1}
deq1 = D[eqs0[[1]], x]

ic2 = First@First@Solve[eqs0[[1]] /. x -> 1, y'[1]] /. y[1] -> 1 /. 
         Rule -> Equal

ysol = y /. 
  First@NDSolve[{deq1, ic2, y[0] == 1/3, y[1] == 1}, y, {x, 0, 1}]

{ysol[0], ysol[1]}
(*   {0.333333, 1.00001}   *)

Plot[ysol[x], {x, 0, 1}]

enter image description here

$\endgroup$

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.