# Using finite differences to solve BVP

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.

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')

• Im trying to run this on matlab, but It is not running. Commented Apr 21, 2019 at 20:34
• You should use Octave
– VoB
Commented Apr 22, 2019 at 9:06
• What line of code should we change to make it work ln matlab? Commented Apr 22, 2019 at 9:53
• @ILoveMath Which line raises an error when you try to run it in Matlab?
– Ian
Commented Apr 23, 2019 at 19:09
• The definition of F and JF Commented Apr 23, 2019 at 22:16