# Solve a differential equation with finite difference method

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

• I checked your analytic solution for both options $f_\infty =\pm\sqrt{\dfrac{\mu}{2a}}$, and none of them satisfy the differential equation. – nicoguaro Sep 4 '14 at 20:11
• I tried to solve the equation in Maple, and the answer was this $$f \left( x \right) ={\frac {k_{{2}}\sqrt {\mu}}{\sqrt {a{k_{{2}}}^{2}- a+\mu}}{\rm sn} \left( {\frac { \left( \sqrt {-2\,a+2\,\mu}x+k_{{1}} \right) \sqrt {\mu}}{\sqrt {a{k_{{2}}}^{2}-a+\mu}}},{\frac {\sqrt {a} k_{{2}}}{\sqrt {-a+\mu}}} \right) } \enspace ,$$ where $\mathrm{sn}$ is the sine-like Jacobi Elliptic function. – nicoguaro Sep 4 '14 at 20:34
• This is a physics contest. $f_\infty$ is equal to 1... in this way $f(x)$ is a solution. The boundary condition are $|f(\pm\infty)| = f_\infty$ – apt45 Sep 4 '14 at 23:09

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 = u
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 - x
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 • $2g$ was $4a$ because $g=2a$. Sorry for this problem... but i stayed all the day doing calculus :) – apt45 Sep 4 '14 at 23:41
• What's the origin of your iteration? It seems like a Crank-Nicolson scheme – apt45 Sep 6 '14 at 16:56
• It is an explicit FDM. You can consider the iteration as a fictitious time. This scheme is based on Jacobi iteration, maybe SOR (Succesive over-relaxation)?? I'm not good with names. – nicoguaro Sep 6 '14 at 17:34