9
$\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$
12
$\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 used of an undocumented syntax of ListPlot here. See this post for more info.

$\endgroup$
10
$\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$
    – Nalt
    May 3 at 4:17
  • $\begingroup$ You're welcome, @Nalt! :) $\endgroup$
    – Michael E2
    May 3 at 12:06

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.