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')