# Finite difference method: numerical solution is not approximate the exact solution.

Given boundary value problem $$$$\label{MKBdir} y''+p(x)y'+q(x)y=r(x),\quad a\leq x\leq b,\quad y(a)=\alpha\text{ and } y(b)=\beta.$$$$

Now I want to solve it numerically using finite difference method (central difference).

We discretize independent variable $$x$$ ($$x_1) with step size $$h$$. Thus, \begin{align}\label{maomao} y_i''+p(x_i)y_i'+q(x_i)y_i=r(x_i). \end{align} with boundary value $$y_1=\alpha\text{ and } y_n=\beta$$.

Substituting central difference formula we have \begin{align*} \dfrac{y_{i+1}-2y_{i}+y_{i-1}}{h^2}+p(x_i)\left(\dfrac{y_{i+1}-y_{i-1}}{2h}\right) +q(x_i)y_i=r(x_i) \end{align*} or equivalently \begin{align*} \left(\dfrac{1}{h^2}-\dfrac{1}{2h}p(x_i)\right)y_{i-1}+ \left(q(x_i)-\dfrac{2}{h^2}\right)y_i +\left(\dfrac{1}{h^2}+\dfrac{1}{2h}p(x_i)\right)y_{i+1}=r(x_i). \end{align*} Let: \begin{align*} a_i&=\dfrac{1}{h^2}-\dfrac{1}{2h}p(x_i)\\ b_i&=q(x_i)-\dfrac{2}{h^2}\\ c_i&=\dfrac{1}{h^2}+\dfrac{1}{2h}p(x_i). \end{align*} We have: \begin{align}\label{shock} a_i y_{i-1}+ b_i y_i +c_iy_{i+1}=r(x_i). \end{align} Now we rewrite the equation for $$i=2,3,\ldots,n-1$$. For $$i=1$$ and $$i=n$$ not need to rewrite since given the boundary value, i.e. $${y_1=\alpha}\text{ and }{y_n=\beta}.$$ For $$i=2$$, \begin{align*} &a_2 {y_{1}}+ b_2 y_2+c_2 y_{3}=r(x_2)\\ \iff& b_2 y_2+c_2 y_{3}=r(x_2)-a_2 {y_{1}} \end{align*} Consider that the value of $$y_1$$ given from boundary value, i.e. $${y_1=\alpha}$$. So, \begin{align*} \boxed{b_2 y_2+c_2 y_{3}=r(x_2)-a_2 {\alpha}.} \end{align*} For $$i=3$$, \begin{align*} \boxed{a_3 y_{2}+ b_3 y_3+c_3 y_{4}=r(x_3).} \end{align*} For $$i=4$$, \begin{align*} \boxed{a_4 y_{3}+ b_4 y_4+c_4 y_{5}=r(x_4).} \end{align*} $$\vdots$$

For $$i=n-3$$, \begin{align*} \boxed{a_{n-3} y_{n-4}+ b_{n-3} y_{n-3}+c_{n-3} y_{n-2}=r(x_{n-3}).} \end{align*} For $$i=n-2$$, \begin{align*} \boxed{a_{n-2} y_{n-3}+ b_{n-2} y_{n-2}+c_{n-2} y_{n-1}=r(x_{n-2}).} \end{align*} For $$i=n-1$$, \begin{align*} &a_{n-1} y_{n-2}+ b_{n-1} y_{n-1}+c_{n-1} {y_{n}}=r(x_{n-1})\\ \iff&a_{n-1} y_{n-2}+ b_{n-1} y_{n-1}=r(x_{n-1})-c_{n-1} {y_{n}}. \end{align*} Consider that the value of $$y_n$$ given from boundary value, i.e. $${y_n=\beta}$$. Thus, \begin{align*} \boxed{ a_{n-1} y_{n-2}+ b_{n-1} y_{n-1}=r(x_{n-1})-c_{n-1} {\beta}.} \end{align*} We get tridiagonal system of linear equation with $$n-2$$ equations and $$n-2$$ variables, i.e:

\begin{align} \begin{bmatrix} b_2&c_2&0&0&0&\cdots&0&0&0&0\\ a_3&b_3&c_3&0&0&\cdots&0&0&0&0\\ 0&a_4&b_4&c_4&0&\cdots&0&0&0&0\\ \vdots&\vdots&\vdots&\vdots&\vdots&\ddots&\vdots&\vdots&\vdots&\vdots\\ 0&0&0&0&0&\cdots&a_{n-3}&b_{n-3}&c_{n-3}&0\\ 0&0&0&0&0&\cdots&0&a_{n-2}&b_{n-2}&c_{n-2}\\ 0&0&0&0&0&\cdots&0&0&a_{n-1}&b_{n-1} \end{bmatrix} \begin{bmatrix} y_2\\ y_3\\ y_4\\ \vdots\\ y_{n-3}\\ y_{n-2}\\ y_{n-1} \end{bmatrix} = \begin{bmatrix} r(x_2)-a_2\alpha\\ r(x_3)\\ r(x_4)\\ \vdots\\ r(x_{n-3})\\ r(x_{n-2})\\ r(x_{n-1})-c_{n-1}\beta \end{bmatrix} . \end{align} Say we have system of linear equation $$AY=B$$, we can find the solution using inverse method: $$Y=A^{-1}B$$.

Now I have boundary value problem:

$$y''+\dfrac{4}{x}y'-\dfrac{2}{x^2}y=-\dfrac{2}{x^2}\ln x$$ for $$1\leq x\leq e$$ with Dirichlet boundary value $$y(1)=\frac{7}{2}$$ and $$y(e)=2e^{-\frac{3}{2}}\cosh\left( \frac{\sqrt{17}}{2}\right)+\frac{5}{2}$$. Given the exact solution $$y=x^{\frac{-3+\sqrt{17}}{2}}+x^{\frac{-3-\sqrt{17}}{2}}+\dfrac{3}{2}+\ln x.$$

Now I want to solve it numerically using central difference method with number of subinterval $$100$$. I solve using MATLAB code as follows:

clear all;
clc;
N=100;
ta=1;
tb=exp(1);
h=(tb-ta)/N;
x=ta:h:tb;
n=length(x);
for i=1:n
a(i)=1/h^2-2/(h*x(i));
b(i)=-2/x(i)^2-2/h^2;
c(i)=1/h^2+2/(h*x(i));
end
%Coefficient matrix
A=zeros(n-2,n-2);
A(1,1)=b(2);
A(1,2)=c(2);
for i=2:n-3
A(i,i-1)=a(i+1);
A(i,i)=b(i+1);
A(i,i+1)=c(i+1);
end
A(n-2,n-3)=a(n-1);
A(n-2,n-2)=b(n-1);
%RHS matrix
B=zeros(n-2,1);
B(1,1)=-2/x(2)^2*log(x(2))-a(2)*7/2;
for i=2:n-3
B(i,1)=-2/x(i+1)^2*log(x(i+1));
end
B(n-2,1)=-2/x(n-1)^2*log(x(n-1))-c(n-1)*2*exp(-3/2)*cosh(sqrt(17)/2)+5/2;

ynum(2:n-1)=A\B;

ynum(1)=7/2;

ynum(n)=2*exp(-3/2)*cosh(sqrt(17)/2)+5/2;

fprintf('  i     xi        ynum_i        yeks_i         error\n');
for i=1:n
yeks(i)=x(i)^((-3+sqrt(17))/2)+x(i)^((-3-sqrt(17))/2)+3/2+log(x(i));
error(i)=abs(ynum(i)-yeks(i));
fprintf('%3d  %5.2f  %10.10f  %10.10f  %10.10f\n',i,x(i),ynum(i),yeks(i),error(i));
end

subplot(2,1,1);
plot(x,ynum,'h','linewidth',1,'color','r','markerfacecolor','c');
hold on;
plot(x,yeks,'linewidth',1,'color','k');
hold on;
title(sprintf('Numerical and exact solution using h=%5.3f',h));
xlabel('x');
ylabel('y');
legend('numerical','exact');
grid on;
subplot(2,1,2);
plot(x,error,'-','color','r','markerfacecolor','g','markersize',5);
grid on;
title('Error');
xlabel('x');
ylabel('y');


After I run the program, I get the result as follows.

Why the exact solution and the numerical solution the error is too high?

I have checked my program over a hour and no mistake in coefficient matrix and constant matrix, all is correct. But the numerical solution is not approximate the exact solution.

Are there mistake in my program?

• Is n equal to N+1? Better use x=linspace(ta,tb,N+1) to get exactly what you want. You could make the code more compact using array operations like p = 4 ./ x; q = -2 ./ x.^2; a = 1/h^2-p/(2*h); etc. Commented Nov 3, 2021 at 16:13
• I have tried to check, n is equal to N+1. >> n==N+1 ans = logical 1 Commented Nov 4, 2021 at 1:20
• I found the mistake in my code. B(n-2,1)=-2/x(n-1)^2*log(x(n-1))-c(n-1)*2*exp(-3/2)*cosh(sqrt(17)/2)+5/2; should be B(n-2,1)=-2/x(n-1)^2*log(x(n-1))-c(n-1)*( 2*exp(-3/2)*cosh(sqrt(17)/2)+5/2 );. Commented Nov 4, 2021 at 1:22

You can avoid many index errors if you construct the occurring array variables via array operations. Here one can compress the code to

  N=100;
ta=1;
tb=exp(1);
h=(tb-ta)/N;
x=linspace(ta,tb,N+1);
% y''+4/x*y'-2/x^2*y = -2/x^2*ln(x)
yeks = x.^(-3/2+sqrt(17)/2)+x.^(-3/2-sqrt(17)/2)+3/2+log(x);
n=length(x);
p = 4 ./ x;
q = -2 ./ x.^2;
a = 1/h^2-p ./ (2*h);
b = q-2/h^2;
c = 1/h^2+p ./ (2*h);

%Coefficient matrix for the equations for x(2) .. x(N)
A=spdiags([ a(3:N+1); b(2:N); c(1:N-1) ]', [-1 0 1], N-1, N-1);
%RHS matrix
B=-2 ./ x.^2.*log(x);
B = B(2:N)';
B(1) -= a(2)*yeks(1);
B(N-1) -= c(N)*yeks(N+1);

ynum = [ yeks(1); A\B; yeks(N+1) ];


One can use several diag commands to construct the tridiagonal matrix from its sub-diagonals, or use the sparse banded matrix constructor spdiags as above. This gives an error level of $$10^{-5}$$ and the plot

• Thank you very much for your help. Your code is more simply. But I have a question. What the mean of -= in B(1) -= a(2)*yeks(1); B(N-1) -= c(N)*yeks(N+1);? If I try to run in MATLAB 2017a, the code gives an error. Commented Nov 4, 2021 at 1:25
• Error: File: add.m Line: 60 Column: 9 The expression to the left of the equals sign is not a valid target for an assignment. Commented Nov 4, 2021 at 1:25
• This was tested with octave. It can be that newer versions restrict the access to array data and require the addressing as column vector, B(k,1) instead of B(k). // += is "add to", -= is "subtract from". Commented Nov 4, 2021 at 6:08