5
$\begingroup$

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?

Improve this question
$\endgroup$
12
  • 1
    $\begingroup$ 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] $\endgroup$
    – mattiav27
    May 29 '16 at 12:05
  • $\begingroup$ 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. $\endgroup$
    – satoru
    May 29 '16 at 12:45
  • 2
    $\begingroup$ @Artes, likely a toy problem, and his actual problem has a qualitatively similar solution. $\endgroup$
    – J. M.'s torpor ♦
    May 29 '16 at 13:52
  • 1
    $\begingroup$ @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. $\endgroup$
    – bbgodfrey
    May 29 '16 at 16:37
  • 2
    $\begingroup$ 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.) $\endgroup$
    – Michael E2
    May 30 '16 at 6:32
10
$\begingroup$

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

enter image description here

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

enter image description here

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.

Improve this answer
$\endgroup$
9
  • $\begingroup$ 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? $\endgroup$
    – satoru
    May 29 '16 at 14:23
  • $\begingroup$ @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? $\endgroup$
    – bbgodfrey
    May 29 '16 at 15:09
  • $\begingroup$ @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 $\endgroup$
    – satoru
    May 29 '16 at 16:34
  • $\begingroup$ 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 $\endgroup$
    – satoru
    May 29 '16 at 16:39
  • $\begingroup$ @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. $\endgroup$
    – bbgodfrey
    May 29 '16 at 16:49
3
$\begingroup$ 
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. enter image description here

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

enter image description here

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

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

enter image description here

Improve this answer
$\endgroup$

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.