Solve a differential equation with finite difference method

Question

I want to solve this equation $$ -\frac{1}{2}f''(x)+2a\ f(x)^3 = f(x)\mu $$ One exact solution (there are a lot of different kinds) of this equation is $f(x) = f_\infty \tanh(\sqrt{2a}f_\infty x) $ (where $2a f_\infty^2=\mu$ ) but I want to plot an approximation of it.

I use the finite difference method to write the derivative.

$f''(x) = \frac{f(x+\Delta x)-2f(x)+f(x-\Delta x)}{(\Delta x)^2} $

and then I have this formula $$ f(x+\Delta x) = 2f(x)-f(x-\Delta x) - 2\Delta x^2f(x)\mu + 4af(x)^3\Delta x^2 $$

I don't know what boundary conditions to use. Can you suggest me something?

I tried solving the equation in the range $[-4,4]$, initializing $f(x)$ in this way $$ f(x) = \begin{cases} -1\qquad -4 < x < 0\\ 0\qquad x=0 \\ 1 \qquad 0 < x < 4 \end{cases} $$ but I get a costant solution $$f(x) = \begin{cases} -1\qquad -4 < x < 0 \\ 1\qquad 0 < x < 4 \end{cases} $$

Improving the code the result is better I see that the correct formula is $$ F(x,\Delta x) = \frac{1}{2}\left[f^{old}(x+\Delta x) + f^{old}(x-\Delta x) \right] + 2a\Delta x^2f^{old}(x) - 2a\Delta x^2f^{old}(x)^3 $$ that is different from the formula suggested by nicoguaro but it is the same he used in his code. My question is: How did you approximate the second derivative? Why my iteration isn't good? Id i use a different step from $0.8$ (as nicoguaro did) i see that the solution goes to the function $sign(x)$. Why? Is the method inconsistent? I used the suggest of nicoguaro and this is my c++ code. Here $g=2a$

#include <iostream>
#include <cmath>
#include <fstream>
using namespace std;
double norm(double * array,int n)
{
    double tmp;
    double norm;
    for (int i = 0; i < n; ++i)
    {
        tmp = *(array + i);
        norm += tmp*tmp;
    }
    return sqrt(norm);
}

double norm(double * array1, double * array2,int n)
{
    double tmp;
    double norm;
    for (int i = 0; i < n; ++i)
    {
        tmp = *(array1 + i) - *(array2+i);
        norm += tmp*tmp;
    }
    return sqrt(norm);
}

double * FD(double *& u, double& g, double& dx, double& mu, int& n)
{
    double * tmp = new double[n];
    *(tmp+n-1) = *(u+n-1);
    *(tmp)=*(u);
    for (int i = 1; i < n-1; ++i)
    {
        tmp[i] = (u[i+1]+u[i-1])/2+dx*dx*u[i]*mu - g*dx*dx*u[i]*u[i]*u[i];
    }

    return tmp;

}
double * dot(double *  u, double& a,int& n)
{
    double * tmp = new double[n];
    for (int i = 0; i < n; ++i) tmp[i] = u[i]*a;
    return tmp;
}

double * sum(double * u1, double* u2, int n)
{
    double * tmp = new double[n];
    for (int i = 0; i < n; ++i) tmp[i] = u1[i]+u2[i];
    return tmp;
}

int main(int argc, char * argv[])
{
    std::cout.setf( std::ios::fixed, std:: ios::floatfield );
    std::fstream fs;
    double dx;
    double xmin,xmax,temp,f0;
    double precision;
    double v,c,k,g,mu;
    int stop = 0;
    cout << "xmin: ";
    cin >> xmin;
    cout << "xmax: ";
    cin >> xmax;
    cout << "passo: ";
    cin >> dx;
    cout << "precisione: ";
    cin >> precision;
    int n = (xmax-xmin)/dx;
    cin >> f0; //initial condition in f[xmin]
    cin >> g; // g = 2a
    cin >> c; // c = f_infinity

    mu = g*c*c; 

    double * fold = new double[n];
    double * fnew = new double[n];
    for (int i = 0; i < n; ++i)
    {
        if(i<n/2) fold[i]=f0;
        if(i>n/2) fold[i] = c;
        if(i==n/2) fold[i] = 0;
    }

    double alpha = 0.3;
    double x = 1-alpha;
while(stop==0)
{
    fnew = sum(dot(fold,alpha,n),dot(FD(fold,g,dx,mu,n),x,n),n);
    if(norm(fold,fnew,n)/norm(fold,n) < precision ) stop = 1;
    else for (int i = 0; i < n; ++i) fold[i] = fnew[i];
}

    fs.open("file.dat", std::fstream::in | std::fstream::out| fstream::app);
    for (int i = 0; i < n; ++i) fs << (xmin+dx*i) << "\t" << fold[i] << endl;

    fs.close();
    return 0;
}

the core of the code is

while(stop==0)
{
    fnew = sum(dot(fold,alpha,n),dot(FD(fold,g,dx,mu,n),x,n),n);
    if(norm(fold,fnew,n)/norm(fold,n) < precision ) stop = 1;
    else for (int i = 0; i < n; ++i) fold[i] = fnew[i];
}

Answer 1

I am turning my comments into this answer. First of all, the analytic expression that you proposed does not satisfy the differential equation. According to Maple, the solution is $$f \left( x \right) ={\frac {k_{{2}}\sqrt {\mu}}{\sqrt {a{k_{{2}}}^{2}- a+\mu}}{\rm sn} \left( {\frac { \left( \sqrt {2\mu - 2a}\ x+k_{{1}} \right) \sqrt {\mu}}{\sqrt {a{k_{{2}}}^{2}-a+\mu}}},{\frac {\sqrt {a} k_{{2}}}{\sqrt {\mu- a}}} \right) } \enspace ,$$ where $\rm{sn}$ is the sine-like Jacobi elliptic function. The parameters $k_1$ and $k_2$ depend on the boundary conditions, let's choose $k_1=k_2=1$ that gives $f(-4) \approx -1$ and $f(4) \approx 1$.

The formula that you used for the finite difference has a problem; the last term should not have the $2$ factor. I prefer another iteration, since yours just approximate the derivatives in one direction. My choose is this one $$ FD(x,\Delta x) \equiv \frac{1}{2}\left[f^{\rm{old}}(x + \Delta x) + f^{\rm{old}}(x - \Delta x)\right] + \mu \Delta x^2 f^{\rm{old}}(x) - 2\Delta x^2 a f^{\rm{old}}(x)^3 \enspace ,$$ and update your function approximation $$f^{\rm{new}}(x) = FD(x,\Delta x) \enspace ,$$ you can also do something like $$f^{\rm{new}}(x) = \alpha f^{\rm{old}}(x) + (1-\alpha) FD(x,\Delta x), \quad \alpha \in [0,1] \enspace .$$

This Python code use this approach (although it don't reach the analytic value).

import numpy as np
import matplotlib.pyplot as plt
from  scipy.special import ellipj

def FD_iter(u, a, mu, dx):
    n = len(u)
    un = np.zeros_like(u)
    un[0] = u[0]
    un[-1] = u[-1]
    alpha = 0.3
    for k in range(1, n-1):
        un[k] = u[k-1]/2 + u[k+1]/2 + dx**2*mu*u[k] - \
                2*dx**2*a*u[k]**3

    return alpha*u + (1 - alpha)*un

plt.close("all")    
n = 1001
a = 1.
mu = 2.
x = np.linspace(-4, 4, n)
u0 = np.tanh(x)
u0 = np.sign(u0)
plt.plot(x, u0, label="Initial")

# This BC lead to a solution with (approx) u(1) = -1  u(4)=1.
k1 = 1
k2 = 1   
v = np.sqrt(mu)*(k1 + np.sqrt(2*mu - 2*a)*x)/np.sqrt(a*k2**2 - a + mu)
m = np.sqrt(a)*k2/np.sqrt(mu - a)
amp = k2*np.sqrt(mu)/np.sqrt(a*k2**2 -a + mu)
z = v*0
ellipj(v, m, z)
z = amp*z
plt.plot(x, z, label="Analytic")


dx = x[1] - x[0]
N = 10000

tol = 1e-5
for k in range(N):
    u = FD_iter(u0, a, mu, dx)
    err = np.linalg.norm(u - u0)/np.linalg.norm(u0)
    if err<=tol:
        break
    u0 = u

plt.plot(x, u, label="Final")
plt.ylim([-1.1, 1.1])
plt.grid(True, color="b", alpha=0.4)
plt.legend()
plt.show()

The curves look like