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