Relaxation methods¶
Solve the linear ODE: $$ y'' + Py' + Qy = F$$ Dirichlet boundary conditions. Discretize the PDE using central differences: $$y'' \approx \frac{y_{i-1}-2y_i+y_{i+1}}{\Delta x^2},$$
$$y' \approx \frac{y_{i+1}-y_{i-1}}{2\Delta x}.$$
Insert into the ODE to get: $$\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_i y_i = F_i.$$
Now group terms: $$ \underbrace{\left(\frac{1}{\Delta x^2}-\frac{P_i}{2\Delta x}\right)}_{l_i}y_{i-1} + \underbrace{\left(Q_i-\frac{2}{\Delta x^2}\right)}_{a_i}y_i + \underbrace{\left(\frac{1}{\Delta x^2}+\frac{P_i}{2\Delta x}\right)}_{u_i}y_{i+1} = \underbrace{F_i}_{b_i}. $$
This gives a tridiagonal system of equations and is in the form
$$ Ay=b. $$
Set the $A$ matrix, and $b$ array, then solve for $y$ using the Thomas Algorithm.
Diagonal elements of $A$ are the $a_i$ coefficients. Upper diagonals are the $u_i$ coefficients. Lower diagonals are the $l_i$ coefficients.
$$
\begin{bmatrix}
a_2 & u_2 & 0 & 0 \\
l_3 & a_3 & u_3 & 0 \\
0 & l_4 & a_4 & u_4 \\
0 & 0 & l_5 & a_5
\end{bmatrix}
\cdot
\begin{bmatrix}
y_2 \\ y_3 \\ y_4 \\ y_5
\end{bmatrix}
=
\begin{bmatrix}
F_2-l_2y_1 \\ F_3 \\ F_4 \\ F_5-u_5y_6.
\end{bmatrix}
$$
The corresponding grid is
| |
* * * * * *
| |
1 2 3 4 5 6 --> 4 interior points (unknowns), 4 eqns.%matplotlib inline
import matplotlib.pyplot as plt
import numpy as np
def relax() :
npts = 21
Ldom = 1.0
yL = 0.5
yR = 0.0
P = 4.0
Q = 0.0
F = -10.0
nvar = npts - 2
x = np.linspace(0,Ldom,npts)
dx = x[1]-x[0]
l = (1.0/dx**2.0 - P/2.0/dx) * np.ones(nvar)
a = (Q - 2.0/dx**2.0) * np.ones(nvar)
u = (1.0/dx**2.0 + P/2.0/dx) * np.ones(nvar)
b = F * np.ones(nvar)
b[0] = b[0] - l[0]*yL
b[-1] = b[-1] - u[-1]*yR
y = thomas(a,u,l,b)
yall = np.hstack([yL,y,yR])
return x,yall
def thomas(a,d,c,b) :
n = np.size(a)
x = np.zeros(n)
for i in np.arange(1,n) :
a[i] = a[i] - c[i]/a[i-1]*d[i-1]
b[i] = b[i] - c[i]/a[i-1]*b[i-1]
x[-1] = b[-1]/a[-1]
for i in np.arange(n-2,-1,-1) :
x[i] = (b[i]-d[i]*x[i+1])/a[i]
return x
x,y = relax()
plt.plot(x,y,'o-')
plt.xlabel('x')
plt.ylabel('y')
Text(0, 0.5, 'y')
