how to solve numerically 2nd order ode with boundary conditions at both ends?

Question

my differential equation has the following form:

$$\frac{d^2\phi}{dx^2}=f(x,\phi',\phi)$$ and it is defined for $0<x<1$.

EDIT: I forgot to give $f(x,\phi',\phi)$ (I think its not that important anyway): $$\frac{d^2\phi}{dx^2}=axI_{1/2}(\frac{\phi}{x})$$,where $a$ is a constant and $I_{1/2}$ is the Fermi-Dirac integral $I_{k}=\int_0^\infty \frac{y^k dy}{1+\exp(y-x)}$. That can be analytically approximated.

In addition I have the following initial values: $$\phi(0)=c, \phi(1)=\phi'(1)=d$$, where $c$ is a constant and $d$ may be adjusted, but there is a first guess for that since $\phi(r)$ usually becomes linear close to $x=1$, see the next image (cannot post more than 2 links yet).

The solution should look similar to this: typical profile of the function $\phi(x) (solid line) $

I tried to solve it numerically using the runge-kutta method (2nd order). Since I don't have an initial value for $\phi'(0)$ I think I have to start from the right at $x=1$ and move to the left.

My idea so far: $$\phi'_n=\phi_{n+1}'- \frac12*h*[f(x_{n+1},\phi'_{n+1})+(x_{n+1}-h)*f(x_{n+1}-h,\phi_{n+1}'-h*f(x_{n+1},\phi'_{n+1}))]$$ $$\phi_n=\phi_{n+1}-\phi'_{n+1}*h-\frac12h^2[f(x_{n+1},\phi'_{n+1})]$$

Here $h$ is the step-size (=1/1000).

However, my result looks like that: $\phi(x)$ (red) and $\phi'(x)$ green. $\phi(x)$ should change most near $x=0$. Of course, in my case it is not, because I don't now how to include the additional boundary-condition $\phi(x=0)=c$ into my code.

Any hints?

Answer 1

There are lots of methods for ODE BVPs. For convenience let's assume that the problem is

$$y'=f(t,y),g(y(0))=0,h(y(1))=0$$

where $y$ is an unknown $\mathbb{R}^2$-valued function of $t$ and $f,g,h$ are given functions. In your case you can construct $f$ in the usual way of encoding second order ODEs as first order systems: $y'_1=y_2,y'_2$ is your second derivative. Then your BCs are generated by $g(x)=x_1-c$, and $h(x)=x_1-x_2$.

Conveniently, this set of BCs completely specifies the value of $y_1(0)$. Because of this, you can solve the BVP by simulating the IVPs:

$$y'=f(t,y),y_1(0)=c,y_2(0)=s$$

where $s$ is a parameter. You then simulate these problems to compute $y_1(1)$ and $y_2(1)$. You then choose a new value of $s$, hopefully to reduce $|y_1(1)-y_2(1)|$. Thus you have reduced the problem of solving the BVP to the problem of 1D root finding: you have a function $r(s)$ which maps $s$ to $y_1(1)-y_2(1)$ where $y'=f(t,y),y_1(0)=a,y_2(0)=s$, and you want to solve $r(s)=0$ for $s$.

In this problem Newton's method is hopeless: $r$ is much too complicated to explicitly differentiate. Bisection is viable but will require some "blind searching" to find a bracket. A secant method is also viable.

There are more sophisticated methods out there; for example, Matlab's bvp4c and bvp5c are collocation-type methods. These will typically see better performance than these "shooting" type methods, but their implementation is more complicated.

Just because I had fun with this question, here's some example code:

function x=secant(f,tol,x0,x1)
    y0=f(x0);
    y1=f(x1);
    while abs(y1)>tol
        slope=(y1-y0)/(x1-x0);
        x2=x1-y1/slope;
        y2=f(x2);
        x0=x1;
        x1=x2;
        y0=y1;
        y1=y2;
    end
    x=x1;
end
function [t,y]=bvp_ode45_shoot(odefun,left_bcfun,right_bcfun,xspan,tol,s0,s1)
    function r=simulate_check_rightbc(s)
        [~,z]=ode45(odefun,xspan,left_bcfun(s));
        r=right_bcfun(z(end,:));
    end
    s=secant(@simulate_check_rightbc,tol,s0,s1);
    [t,y]=ode45(odefun,xspan,left_bcfun(s));
end

and a test case:

odefun=@(t,y) [y(2);-y(1)];
left_bcfun=@(s) [0;s];
right_bcfun=@(r) r(1)-1;
xspan=[0 1];
tol=1e-2;
s0=1;
s1=2;
[t,y]=bvp_ode45_shoot(odefun,left_bcfun,right_bcfun,xspan,tol,s0,s1);

The code implements a shooting method using the 1D secant method for the search. The test case solves $y''=-y$ subject to $y(0)=0,y(1)=1$ beginning with $y'(0)=1$ and $y'(0)=2$ as initial guesses. One can compare the results of this program to the explicit solution $y=\frac{\sin(x)}{\sin(1)}$.