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

Question

Given boundary value problem \begin{equation}\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. \end{equation}

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

We discretize independent variable $x$ ($x_1<x_2<...<x_{n-1}<x_n$) 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?

Answer 1

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