Discretize the domain into a finite difference grid. That is, a set of discrete points upon which we will find the solution.
| : : |
O O O : O : O O O
| i-1 i i+1 |
x=0 x=R
Approximate the exact derivatives in the ODE by finite difference approximations (FDS)
Substitute the FDAs into the ODE to get the finite difference equation (FDE). $$\frac{y_{i+1} - 2y_i + y_{i-1}}{\Delta x^2} + P_i\frac{y_{i+1}-y_{i-1}}{2\Delta x} + Q_iy_i = F_i.$$
Solve the system of algebraic equations.
The Thomas Algorithm needs arrays $l$, $a$, $u$, and $b$.
For the above matrix $n=7$ so $i$ varies from $0$ to $n-1=6$, and our unknowns are $y_1$ to $y_{n-2} = y_{5}$. Our interior points are $y_2$ to $y_{n-3}=y_4$.
Solve the following with Dirichlet boundary conditions. $$y^{\prime\prime} + Py^{\prime} + Qy = F.$$
import numpy as np
import matplotlib.pyplot as plt
%matplotlib inline
from thomas import *
ngrd = 21
Ldom = 1.0
yL = 0.5
yR = 0.0
P = np.full(ngrd, 4.0)
Q = np.zeros(ngrd)
F = np.full(ngrd, -10.0)
x = np.linspace(0, Ldom, ngrd)
dx = x[1] - x[0]
l = 1.0/dx**2 - P/2/dx
a = Q - 2/dx**2
u = 1.0/dx**2 + P/2/dx
b = F
b[1] -= l[1]*yL
b[ngrd-2] -= u[ngrd-2]*yR
y = thomas(a[1:-1], u[1:-1], l[1:-1], b[1:-1])
y = np.insert(y,0,yL)
y = np.append(y,yR)
plt.rc('font', size=14)
plt.plot(x,y, 'bo-')
plt.xlabel('x')
plt.ylabel('y');
a = [1,2,3,4,5]
a[1:-1]
[2, 3, 4]