Using finite differences to solve BVP

Question

I have the following ODE

$$ u'' = -(1 + e^u), u(0)=0, u(1)=1$$

Using a second order accurate finite difference I obtain

$$ -(1+e^{u_i}) \approx \frac{ u_{i+1} - 2 u_i + u_{i-1} }{h^2} $$

and $u_0 = 0$ and $u_N=1$. Let's say we use mesh size $N$ so we will obtain $N-2$ nonlinear equations for unknowns $u_2,...,u_{N-1}$. Im looking for some help to implement this in Matlab as it would be quite laborious to do it by hand.

Answer 1

After the discretization, you end up with $-(1+e^{u_i}) = \frac{ u_{i+1} - 2 u_i + u_{i-1} }{h^2}$, which is equivalent to find the zeros of the vector valued function $F(u_i)=(1+e^{u_i})+ \frac{ u_{i+1} - 2 u_i + u_{i-1} }{h^2}$ where $2\leq i\leq N-1$

What you are missing is to use a Newton's method in order to solve this system of non-linear equations.

Consider the tridiagonal matrix which arises from the discretization of the second derivative, say it $A$. (But pay attentio at the boundary conditions, you should modify the first and last row)

Then, you can write $\frac{ u_{i+1} - 2 u_i + u_{i-1} }{h^2}$ as a matrix-vector multiplication $A \cdot \vec{u}$.

With this idea in mind, we rewrite the function $F$ by using that matrix:

---

$F(u)=A\cdot u + \mathbf{1}+e^{\mathbf{u}}$ (note that the RHS is a (N-2)-vector)

Now you just need to use Newton method for the non-linear system, and so you need to compute the Jacobian, which, in this particular case, is

$JF(u)=A+diag(e^u)$

Now just apply Newton's routine and you'll get your numerical solution.

|EDIT|

The following runable Octave code show the correct numerical solution (cfr with the previous question)

clear all
close all

m=201;
h=1/(m-1);
x=linspace(0,1,m)';
A = toeplitz(sparse([1,2],[1,1],[-2,1]/(h^2),m,1));

F=@(u) [u(1);(A*u+ones(m,1)+exp(u))(2:m-1);u(m)-1];
JF=@(u) [[1,zeros(1,m-1)];(A+diag(exp(u)))(2:m-1,1:m);[zeros(1,m- 
1),1]];
u0=ones(m,1); %Starting guess for newton's method
res=-JF(u0)\F(u0);
tol=h^2/100;
while (norm(res,inf)>tol)
     u0+=res;
        res=-JF(u0)\F(u0);
end
u0+=res;

plot(x,u0,'b-o')

---

|EDIT|$^2$

Note that the following way to impose boundary conditions is in some sense naïf, there are better ways

In order to write a runable MatLab code, I modify just the first and last rows of the function that we have to set to zero.

First of all define the vector of ones $b=\mathcal{1}$, so the system reads:

$Au+b+e^u=0$ (again, with $e^u$ I mean the vector $[e^{u_1},\ldots,e^{u_n}]$

And, I want to have $u_1=0$. To this aim, I impose the first row of the $A$-matrix to have just a $1$ in the (1,1)-position: in this way only $u_1$ survives.

I do this with the command A(1,1:2)=[1,0]

Since I have only a dependence on $u$ on the exponential term, I multiply $u$ by an identity matrix $I$, with $I(1,1)=0$, in order to neglect the term $e^{u_1}$. I do this with the command I(1,1)=0

So far, the first line is $u_1+ b(1) + e^0=0$, or better:

\begin{align} u_1+ b(1) + 1=0 \end{align}

Now, I give the first component of $b$ the value I need in order to have $u_1=0$, which is simply $-1$ ! I do this with the command b(1)=-1

The first line of the system reads now $u_1 +1-1=0$, hence $u_1=0$, as wanted.

The same reasoning applies to the last row, you can see how in the following code.

clear all
close all

m=201;
h=1/(m-1);
x=linspace(0,1,m)';
A = toeplitz(sparse([1,2],[1,1],[-2,1]/(h^2),m,1));
I=speye(m); %initialize sparse matrix of size m
b=ones(m,1);
%Boundary conditions
A(1,1:2)=[1,0];
I(1,1)=0; 
b(1)=-1;

A(m,m-1:m)=[0,1];
I(m,m)=0;
b(m)=-2;

%Functions to set to zero and Newton's loop
F=@(u) A*u+b+exp(I*u);
JF=@(u) A+I*diag(exp(I*u));

u0=ones(m,1); %Starting guess for newton's method
res=-JF(u0)\F(u0);
tol=h^2/100;
while (norm(res,inf)>tol)
    u0=u0+res;
    res=-JF(u0)\F(u0);
end
u0=u0+res;
plot(x,u0,'b-o')