Solve a Two-point ODE Boundary Value Problem (BVP) with Finite Difference Method (FDM)
ODE-BVP (see Example 7.8 in page 277)
Using ODE theory (method of undetermined coefficients), we can find its exact solution:
with
c1=(3-exp(-2))/(exp(2)-exp(-2));c2=(exp(2)-3)/(exp(2)-exp(-2));
ufun=@(x) c1*exp(2*x)+c2*(exp(-2*x));%exact solution
ffun=@(x) zeros(size(x)); %right-hand size
Following our Lecture of Section 7.2, we first define a uniform mesh of [0,1]:
Here we divide it into (m+1) sub-intervals.
m=5; h=1/(m+1);%step size
x=linspace(0,1,m+2)'; %divide [0,1] into (m+2) points
figure(1)
plot(x,ufun(x),'o--'); %plot the exact solution
xlabel('x'); ylabel('u(x)');
Then we will discretize the ODE BVP by the Central FDM to get a sparse linear system.
Since the ODE itself holds at each interior point which gives
, with the boundary conditions
Now, using the central finite difference formula for approximating to get
.
Let and , then the above ODE at each point now reads
.
Upon dropping the truncation error term, we get a linear system of finite difference formulas (see lecture for the matrix formulation pattern)
,
where denotes the finite difference numerical approximation.
By formualting this into a sparse system of linear equations , where

we can solve the linear system for , which gives the numerical solution.
Under certain assuptions, we can show it converges to the exact solution in the sense that
as , (other vector norms can also be used)
where is the exact solution at the same grid points.
Below is the vectorized codes for building this system using sparse matrices (based on spdiags):
e=ones(m,1);
Ah=spdiags([e/h^2 -2*e/h^2-4 e/h^2],(-1:1),m,m);%tridiagonal matrix
fh=ffun(x(2:end-1));%define rhs over interior nodes
fh(1)=fh(1)-1/h^2;fh(end)=fh(end)-3/h^2;%enforce boundary condition
Once the system is constructed, we can for now solve it using backslash solver (sparse direct solver)
Uh=Ah\fh;%Solve the system by backslash solver, finite different solution
uh=ufun(x(2:end-1));%exact solution over interior nodes
maxerr=norm(Uh-uh,"inf") %this error should of O(h^2)
maxerr = 0.0045
The computed maxerr is of , which means it will be four time smaller if h is halved.
Plot the finite different solution versus the exact solution, see almost not difference
figure(2);hold on; plot(x(2:end-1),uh,'o--'); %plot the exact solution
plot(x(2:end-1),Uh,'x--r') %plot the finite difference solution
title('Exact solution vs finite different solution')
It is sometimes better to see their difference by plotting the point-wise errors (no errors at the boundary)
figure(3); plot(x(2:end-1),Uh-uh,'-s')
xlabel('x');title('point-wise error at interior nodes')
Often, one can use loglog plot of maxerr as a function of h to verify the order of accuracy.
Notice maxerr=C*h^2 implies log(maxerr)=log(C)+2*log(h), which is a log-log line with slope 2.
Below is a summary of the above codes within a loop of refined mesh sizes to generate an error log-log plot,
where you need to understand that h limits to zero when you look from right to the left.
figure(4)
mlist=2.^(5:10);
hlist=1./(mlist+1);
maxerr=zeros(size(mlist));
for k=1:length(mlist)
m=mlist(k); h=1/(m+1); x=linspace(0,1,m+2)'; e=ones(m,1);
Ah=spdiags([e/h^2 -2*e/h^2-4 e/h^2],(-1:1),m,m);%tridiagonal matrix
fh=ffun(x(2:end-1));%define rhs over interior nodes
fh(1)=fh(1)-1/h^2;fh(end)=fh(end)-3/h^2;%enforce boundary condition
Uh=Ah\fh;%Solve the system by backslash solver, finite different solution
uh=ufun(x(2:end-1));%exact solution over interior nodes
maxerr(k)=norm(Uh-uh,"inf"); %this error should of O(h^2)
end
p = polyfit(log(hlist),log(maxerr),1); slope=p(1) %Fix by a line to check the slop
slope = 1.9999
loglog(hlist,maxerr,'d-'); xlabel('h');ylabel('maxerr');%should give a straight line with slop 2