2
$\begingroup$

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? grap 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];
}
Improve this question
$\endgroup$
  • 1
    $\begingroup$ I checked your analytic solution for both options $f_\infty =\pm\sqrt{\dfrac{\mu}{2a}}$, and none of them satisfy the differential equation. $\endgroup$ – nicoguaro♦ Sep 4 '14 at 20:11
  • $\begingroup$ 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. $\endgroup$ – nicoguaro♦ Sep 4 '14 at 20:34
  • $\begingroup$ 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$ $\endgroup$ – apt45 Sep 4 '14 at 23:09
3
$\begingroup$

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 Comparison of solutions

Improve this answer
$\endgroup$
  • $\begingroup$ $2g$ was $4a$ because $g=2a$. Sorry for this problem... but i stayed all the day doing calculus :) $\endgroup$ – apt45 Sep 4 '14 at 23:41
  • $\begingroup$ What's the origin of your iteration? It seems like a Crank-Nicolson scheme $\endgroup$ – apt45 Sep 6 '14 at 16:56
  • $\begingroup$ 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. $\endgroup$ – nicoguaro♦ Sep 6 '14 at 17:34

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.