# Stiff second order ODE

I am trying to solve numerically $f''=f^{3}-f$, whose exact solution is $\tanh(x)$. The problem is that numerical solution fails if I come closer to $\tanh$ plateau. "StiffnessSwitching" method doesn't help.

fff = NDSolve[{f''[x] == f[x]^3 - f[x], f[-3] == 1, f[3] == -1}, f[x],
{x, -3, 3}, Method -> "StiffnessSwitching"]
N[Tanh[3]]

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


What am I doing wrong?

• Are you sure the solution is Tanh? FullSimplify[(f''[x] - f[x]^3 + f[x]) /. {f -> Function[{x}, Tanh[x]]}] = -Sech[x]^2 Tanh[x] – mattiav27 May 29 '16 at 12:05
• Direct Calculation : $Tanh''[x]=-2 Sech[x]^2 Tanh[x]=Tanh[x]^3-Tanh[x]$Divide on tanh u see $Sech[x]^2=1-Tanh[x]^2$ u get identity. – satoru May 29 '16 at 12:45
• @Artes, likely a toy problem, and his actual problem has a qualitatively similar solution. – J. M. is away May 29 '16 at 13:52
• @Artes DSolve gives me a solution (in terms of JacobiSN), only when I discard the boundary conditions, in which case I am left with two unknown constants that are not easy to evaluate. Have you done better? Thanks. – bbgodfrey May 29 '16 at 16:37
• Generally speaking, an NDSolve::ndsz error in a (nonsingular, non-stiff) BVP is (most likely) coming from a singularity developing in the initial conditions chosen by the built-in shooting method used by NDSolve. Explicitly setting "StartingInitialConditions" can help, when possible. So can manually implementing the shooting method. (This should be better known/advertised. It's come up several times.) – Michael E2 May 30 '16 at 6:32

The equation is not stiff, despite the claim by Mathematica. It can be solved by using the "Shooting" Method.

NDSolve[{f''[x] == 2 (f[x]^3 - f[x]), f[-3] == 1, f[3] == -1}, f, {x, -3, 3},
Method -> {"Shooting", "StartingInitialConditions" -> {f[-3] == 1, f'[-3] == 0}}][[1, 1]];
Plot[f[x] /. %, {x, -3, 3}, AxesLabel -> {f, x}]


This looks like a plot of -Tan[x] but is not precisely that because of the boundary conditions.

Addendum - Sample Solution for Question in Comment

satoru describes a much more difficult problem in a comment below. A sample solution to it is

xm = 10;
NDSolve[{f''[x] == f[x]^3 - f[x] + α f[x] (1 - Tanh[x])^2,
f[-xm] == Sqrt[1 - 2 α], f[xm] == -1} /. α -> 10^-6, f, {x, -xm, xm},
WorkingPrecision -> 30, Method -> {"Shooting",
"StartingInitialConditions" -> {f[xm] == -1, f'[xm] == 0}}][[1, 1]];
Plot[f[x] /. %, {x, -xm, xm}, AxesLabel -> {f, x}]


Note how the Tanh-like solution is shifted to the left even for tiny α. Solutions for larger a probably can be obtained by centering the integration range on the value of x for which f[x] == Sqrt[1 - 2 α] - 1, and obtaining an estimate for that quantity by extrapolating results from smaller values of α. See, for instance, the automated process in my answer here.

• Wow, this realy works! But the new question is what to do if i dont know exact solution and cant use precise boundary conditions? Real problem is $f''=f^{3}-f+\alpha f(1-Tanh[x])^{2}$ I know how solution behaves on infinities, external field part becomes constant and i have equation i wrote above. So $f[-\infty]=-\sqrt{1-2\alpha}Tanh[\sqrt{1-2\alpha}(-\infty)]=\sqrt{1-2\alpha}$ and $f[\infty]=-1$ In my first example i used 3 instead of infinity as some big number, Tanh[3] is near one, and this dont work, how to plug infinities in boundary conditions? – satoru May 29 '16 at 14:23
• @satoru I was careless in my answer. Actually, boundary conditions using Tanh[]are not needed. I shall correct my answer in a few minutes. The question you pose in your comment is much harder. Usually, choosing some large value instead of infinity is sufficient for the boundary conditions. What value of alpha would you like? – bbgodfrey May 29 '16 at 15:09
• @bbgodfrey i began from small alpha (and sign before external part must be -), but something is wrong : a = 5 \[Alpha] = 0.1 s = NDSolve[{f''[x] == 2 *(f[x]^3 - f[x]) - \[Alpha]*f[x]*(1 - Tanh[x])^2, f[-a] == -Sqrt[1 + 2 \[Alpha]] Tanh[-a*Sqrt[1 + 2 \[Alpha]]], f[a] == -Tanh[a]}, f, {x, -a, a}, Method -> {"Shooting", "StartingInitialConditions" -> {f[-a] == -Sqrt[1 + 2 \[Alpha]] Tanh[-a*Sqrt[1 + 2 \[Alpha]]], f'[-a] == 0}}] Plot[Evaluate[f[x] /. s], {x, -a, a}, PlotRange -> All] stiffness report – satoru May 29 '16 at 16:34
• I also noticed if a use too large a (about 10) even with alpha=0 and exact boundary conditions solution becomes wrongly dislocated from x=0 – satoru May 29 '16 at 16:39
• @satoru Use WorkingPrecision -> 30 to solve the a == 10 problem, although the code runs much slower. I am not surprised that α == 0.1 fails. Try α == 10^-6 and you will already see a large shift. I shall post the latter answer, but solving your problem more generally will be quite difficult. – bbgodfrey May 29 '16 at 16:49
ClearAll["Global*"];
Remove["Global*"];

sol = NDSolve[{f''[x] == f[x]^3 - f[x], f[-3] == 1, f[3] == -1}, f[x], {x, -3, 3},
Method -> {"Shooting", "StartingInitialConditions" -> {f'[0] == 0, f[0] == 1}}];


UPDATE:

Maples symbolic solution(The symbolic solution is quite long and complex) ,converted to numeric form.Maple gives 5 solutions other than Mathematica.

Maples JacobiSN[z,k] is equal to Mathematica JacobiSN[z,k^2].

With 20 digits precision.

maple1 = (-1.5279997865913399554 - 0.58636062345487262098*I)*
JacobiSN[(0.94894046904339260202 -
0.94416766565747605255*I)*(0.70710678118654752440*x +
1.6985155326343786960 -
0.20793350722877537103*I), (-0.50021641106619559220 -
1.1156113783217666782*I)^2];
maple2 = (-0.57150532368340334446 -
0.87200104339127158745*I)*
JacobiSN[(1.5911818830635583675 -
0.31319690342130914617*I)*(0.70710678118654752440*x +
1.7471798211151876014 +
0.44795429731942814652*I), (-0.24192870056640322887 -
0.59564049424232431474*I)^2];
maple3 = (-1.5279997865913399554 -
0.58636062345487262098*I)*
JacobiSN[(0.94894046904339260202 -
0.94416766565747605255*I)*(0.70710678118654752440*x +
1.6985155326343786960 -
0.20793350722877537103*I), (-0.50021641106619559220 -
1.1156113783217666782*I)^2];
maple4 = (-1.5279997865913399554 -
0.58636062345487262098*I)*
JacobiSN[(0.94894046904339260202 -
0.94416766565747605255*I)*(0.70710678118654752440*x +
1.6985155326343786960 -
0.20793350722877537103*I), (-0.50021641106619559220 -
1.1156113783217666782*I)^2];
maple5 = (-1.5279997865913399554 +
0.58636062345487262098*I)*
JacobiSN[(0.94894046904339260202 +
0.94416766565747605255*I)*(0.70710678118654752440*x +
1.6985155326343786960 +
0.20793350722877537103*I), (-0.50021641106619559220 +
1.1156113783217666782*I)^2];


.

Boundary conditions check:

  Re[{maple1, maple2, maple3, maple4, maple5}] /. x -> -3 // N
(*{1., 1., 1., 1., 1.}*)
Re[{maple1, maple2, maple3, maple4, maple5}] /. x -> 3 // N
(*{-1., -1., -1., -1., -1.}*)


.

  Plot[{Re[maple1], Re[maple2], Re[maple2], Re[maple2], Re[maple2],
f[x] /. sol}, {x, -3, 3},
PlotLegends -> {"maple1", "maple2", "maple3", "maple4", "maple5",
"NDSOLVE"},
PlotStyle -> {Red, {Green, Dashing[{0.2, 0.05}],
Thickness[0.01]}, {Blue, Thickness[0.01],
Dashing[{0.3, 0.1}]}, {Black, Thickness[0.01],
Dashing[{0.1, 0.1}]}, Yellow, Brown}, AxesLabel -> {x, f[x]}]


• How you converted the long symbolic expression into numeric in maple? – zhk May 29 '16 at 14:47
• @MMM. I use a evalf command. – Mariusz Iwaniuk May 29 '16 at 17:11
• Have you checked the residual error of Maple's solution (with precision higher than MachinePrecision)? – xzczd Jan 31 '17 at 3:28
• @xzczd .It seems that everything is ok.Yes i checked.Boundary conditions ok. – Mariusz Iwaniuk Jan 31 '17 at 10:20

In versiion 12.0 you can use the nonlinear FEM solver:

sol = NDSolveValue[{f''[x] == f[x]^3 - f[x], f[-3] == 1, f[3] == -1},
f[x], {x, -3, 3}, Method -> "FiniteElement"]
Plot[sol, {x, -3, 3}, AxesLabel -> {f, x}]