DSolve to obtain a tanh solution

We are given a simple ODE with BCs:

$$\xi^2 \frac{df^2}{dx^2} + f - f^3 = 0$$

$$f(x=0) = 0$$

$$\\f(x\to\infty) = 1$$

On paper this is quite easy to solve. One can obtain the solution

$$f(x) = \operatorname{tanh}\left(\frac{x}{\sqrt{2}\xi}\right).$$

The actual physics of the problem needs not a particular boundary at infinity, but just that far away from $$x=0$$ (or many $$\xi$$ lengths away), the value of $$f$$ tends to $$1$$ (here, $$f$$ has been normalised already).

This problem is a special (zero-field) case of a more general problem, so later I will try and generalise this approach to something that I don't know the analytic solution to easily. However, I am having trouble even obtaining this simple case in Mathematica. I am trying to use DSolve to find an analytic solution. I show three attempts.

Solving without the boundary conditions, to implement later

First, the standard approach

ZFeqn = xi^2 * f''[x] + f[x] - f[x]^3 == 0

sol = DSolve[{ZFeqn}, f[x], x]


However, this gives me a horrid result, in terms of Jacobi functions.

Solving with the boundary conditions

Secondary, with the BCs as well

ZFeqn = xi^2 * f''[x] + f[x] - f[x]^3 == 0

sol = DSolve[{ZFeqn, f[0] == 0, f[Infinity] == 1}, f[x], x]


This runs for a little while (1 minute) and then gives up. I have seen some previous posts discuss how to explicitly deal with boundary conditions at infinity, but to avoid that...

Implementing the boundary condition by hand

As is generally a good idea, we can manipulate first. We multiply by $$f'$$ and perform the integral, applying the BC at infinity to find the integration constant ($$1/2$$). We now solve the 1st order ODE with one BC

ZFeqn = xi^2 * (f'[x])^2 + f[x] - (1/2)*f[x]^3 - 1/2 == 0

sol = DSolve[{ZFeqn, f[0] == 0}, f[x], x]


However, even this finds a solution in terms of inverse elliptic functions - not the simple tanh function I expect!

Any further wisdom would be greatly appreciated.

EDIT: I have put this third equation into Maple, and have found the solution required... I am very interested to know why Mathematica couldn't/shouldn't be able to handle this!

EDIT2: Removed the screenshot of Maple since there was a typo in the factors. Easiest route is to square root each side of the equation and get the first order ODE in the form f'(x) rather than (f(x)')^2, which Mathematica seemingly struggles to solve. Problem is now solved; thanks to all of Michael, Ulrich and Alexei!

• Mathematica reproduces the Maple result thus: DSolve[xi^2*(f'[x]) + (1/2)*f[x]^2 - 1/2 == 0, f[x], x] // FullSimplify. The ODE you fed Maple does not seem equivalent to the first ZFeqn above, but maybe I'm blind to some obvious transformation. – Michael E2 May 10 at 15:11
• Dear Michael - thanks for your response. There's a small typo in your equation above (I see you've basically taken the square root of each side but have omitted the square roots of the 1/2 factor and xi^2). This is then equivalent to what has been fed into Maple, and is therefore the same as the third and final ZFeqn above. However, this answers my question of why the previous attempts were not finding the solutions that I expected - I forgot that MMA does not like equations of the type (f'[x])^2 = ...! – Brad May 10 at 15:32
• Actually I just copied the equation from the Maple output (1), so the error originates there, which was really my point. Or rather: the Maple code and Mathematica code do not seem equivalent. – Michael E2 May 10 at 15:55
• Whoops... the mistake was fully mine in there then; I've made another typo somewhere, but it's only an overall factor or so out... Regardless, after your happy reminder that MMA likes to solve equations in a particular form. I'll remove the image; good spot, thank you! – Brad May 10 at 15:59

This is a well-known problem of a rigid domain wall called after Zhirnov, who solved it first. The main point here is that the domain wall solution (that you are, in fact, looking for) corresponds to a separatrix. That is, it is a special trajectory. To select it, one needs to take the corresponding boundary conditions. It brings up the desired simple solution, while other trajectories are described by special functions.

Let us be more precise.

First, let us rescale the coordinate x so as to turn the coefficient \[Xi] into 1. From here on, we will assume this to be done.

Second, let us take the left-hand part of your initial equation:

expr1 = f''[x] + f[x] - f[x]^3;


multiply it by f'[x] and integrate over x. The result let us equate to an arbitrary constant, C1:

eq = Distribute[Integrate[#, x] &[f'[x] expr1]] == C1

(*   f[x]^2/2 - f[x]^4/4 + 1/2 Derivative[1][f][x]^2 == C1   *)


To make it better visible, I show it once more as an image:

Now to clearly see the origin of the special choice of the boundary condition, let us look at this problem from another point of view. Equations that we solve are equivalent to dynamic equations describing the motion of a material point with the mass m=1 in the potential field, U=U(f) with

U= f^2/2 - f^4/4


It is shown by the code

Show[{
Plot[(f^2/2 - f^4/4), {f, -1.5, 1.5},
AxesLabel -> {Style["f", 16, Italic], Style["U", 16, Italic]}],
Graphics[{Red, PointSize[0.03], Point[{0, 0}], Green,
PointSize[0.03], Point[{1, 0.26}]}]
}]


in the image below:

Here f plays the role of the coordinate of this moving material point and f' - of its velocity. The equation eq then represents its total energy.

Your boundary conditions correspond to such a trajectory of this point that starts from the minimum of the potential well (indicated by the red point) with some non-zero kinetic energy and completely stops at the top of the potential hill (the green point). At this point, it has the velocity, f' equal to zero. Otherwise, it does not stop. This point f=1 corresponds to x=Infinity. As we already mentioned, here f'=0. One finds

eq /. {f'[x] -> 0, f[x] -> 1}

(*  1/4 == C1  *)


Thus, the equation describing the separatrix solution is

eq2 = eq /. C1 -> 1/4

(*  f[x]^2/2 - f[x]^4/4 + 1/2 Derivative[1][f][x]^2 == 1/4  *)


Now one can solve it straightforwardly:

dsl=DSolve[eq2, f, x]

(*  {{f -> Function[{x}, (E^(Sqrt[2] x) - E^(2 C[1]))/(
E^(Sqrt[2] x) + E^(2 C[1]))]}, {f ->
Function[{x}, (1 - E^(Sqrt[2] x + 2 C[1]))/(
1 + E^(Sqrt[2] x + 2 C[1]))]}}   *)


There are two solutions. Let us transform them a bit:

Simplify[dsl[[1, 1, 2, 2]] // ExpToTrig]

(*  Tanh[x/Sqrt[2] - C[1]]   *)


and

Simplify[dsl[[2, 1, 2, 2]] // ExpToTrig]

(*  -Tanh[x/Sqrt[2] + C[1]]  *)


First, we see that two solutions correspond to the two-fold degeneration of the equation due to its mirror symmetry. Second, the integration constant, C[1], only describes a shift of the domain wall origin. To satisfy your boundary condition, you should put C[1]=0.

Have fun!

• Dear Alexei - thank you for your very detailed answer! Very clear and helpful, particularly in context with the answer from Ulrich! Thank you for your time – Brad May 10 at 15:22

If you take your third step "multiply ode by f'[x]" and integrate in the range {x,0,\[Infinity]} it follows

\[Xi]^2/2 (f'[\[Infinity]]^2-f'[0]^2)+f[\[Infinity]]^2/2-f[\[Infinity]]^4/4==0

This expression (Conservation of energy) is different to your result!

Assumption f'[\[Infinity]]==0, which must be true because f[\[Infinity]]==1, leads to f'[0]== 1/(Sqrt[2]xi)

Now you know the right ics and get

F = DSolveValue[{ZFeqn, f[0] == 0, f'[0] == 1/(Sqrt[2] xi) }, f , x]
(*Function[{x}, (-1 + E^((Sqrt[2] x)/xi))/(1 + E^((Sqrt[2] x)/xi))]*)


This corresponds to the simple result you're expecting!

• Dear Ulrich - thank you for your answer; very helpful. I am surprised that this approach is not equivalent to, say the picture detailed in Maple, which was much easier to simply type in. However, this provides a full solution from the original equation with just a trivial rewriting of the boundary conditions, which is certainly useful. What made you think of rewriting the BCs like this? I would've expected my approach to produce the same results? – Brad May 10 at 14:43
• @Brad Knowledge of "consevation of energy" often helps to transform a boundary problem to an initial value problem. By the way how did you derive the ode xi^2 * (f'[x])^2 + f[x] - (1/2)*f[x]^3 - 1/2 == 0? – Ulrich Neumann May 10 at 14:48
• Perhaps you could explain how this idea of CoE applies here? I'm familiar with the concept, and certainly agree with your maths, but I've never considered how the idea of reducing BCs being described in an energy conservation manner. The ODE you state was derived in a very similar way to your procedure actually. By using the idea that f -> 1 at infinity, we therefore have that f' -> 0 at infinity. Hence by performing the integral, we obtain the indefinite result \[Xi]^2 f'[x]^2+f[x]^2 - f[x]^4/2=const. Evaluating at infinity, (f=1, f' = 0) leads to this constant being 1/2, – Brad May 10 at 14:57
• @Brad Same idea as yours. You integrated in the range x...Infinity and I took 0...Infinity` – Ulrich Neumann May 10 at 15:06
• Ah of course it is. Of course, the approaches would have had to have been equivalent, but my surprise was more than MMA wasn't able to solve the 1st order ODE that I proposed in my final code snippet in the question above (involving the ODE that you have asked). Implicitly, it seems that I have already done some of the work in solving it! – Brad May 10 at 15:10