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?
Best Answer
You can avoid many index errors if you construct the occurring array variables via array operations. Here one can compress the code to
One can use several
diagcommands to construct the tridiagonal matrix from its sub-diagonals, or use the sparse banded matrix constructorspdiagsas above. This gives an error level of $10^{-5}$ and the plot