Finite Difference Matrix

Many matrices in scientific computing contain mostly zeros, particularly those arising from the discretistaion of partial differential equations (PDEs). Here we will construct a sparse matrix using scipy.sparse that is derived from the finite difference discretistaion of the Poisson equation. In 1D, Poisson equation is

uxx=f(x) for 0≤x≤1u_{xx} = f(x)\text{ for }0 \le x \le 1

The central FD approximation of uxxu_{xx} is:

uxx≈u(x+h)−2u(x)+u(x−h)h2u_{xx} \approx \frac{u(x + h) - 2u(x) + u(x-h)}{h^2}

We will discretise uxx=0u_{xx} = 0 at NN regular points along xx from 0 to 1, given by x1x_1, x2x_2:

          +----+----+----------+----+> x
          0   x_1  x_2    ... x_N   1

Using this set of point and the discretised equation, this gives a set of NN equations at each interior point on the domain:

vi+1−2vi+vi−1h2=f(xi) for i=1...N\frac{v_{i+1} - 2v_i + v_{i-1}}{h^2} = f(x_i) \text{ for } i = 1...N

where vi≈u(xi)v_i \approx u(x_i).

To solve these equations we will need additional equations at x=0x=0 and x=1x=1, known as the boundary conditions. For this example we will use u(x)=g(x)u(x) = g(x) at x=0x=0 and x=1x=1 (also known as a non-homogenous Dirichlet bc), so that v0=g(0)v_0 = g(0), and vN+1=g(1)v_{N+1} = g(1), and the equation at x1x_1 becomes:

vi+1−2vi+g(0)h2=f(xi)\frac{v_{i+1} - 2v_i + g(0)}{h^2} = f(x_i)

and the equation at xNx_N becomes:

g(1)−2vi+vi−1h2=f(xi)\frac{g(1) - 2v_i + v_{i-1}}{h^2} = f(x_i)

We can therefore represent the final NN equations in matrix form like so:

1h2[−211−21⋱⋱⋱1−211−2][v1v2⋮vN−1vN]=[f(x1)f(x2)⋮f(xN−1)f(xN)]−1h2[g(0)0⋮0g(1)] \frac{1}{h^2} \begin{bmatrix} -2 & 1 & & & \\ 1 & -2 & 1 & & \\ &\ddots & \ddots & \ddots &\\ & & 1 & -2 & 1 \\ & & & 1 & -2 \end{bmatrix} \begin{bmatrix} v_1 \\ v_2 \\ \vdots \\ v_{N-1}\\ v_{N} \end{bmatrix} = \begin{bmatrix} f(x_1) \\ f(x_2) \\ \vdots \\ f(x_{N-1}) \\ f(x_N) \end{bmatrix} - \frac{1}{h^2} \begin{bmatrix} g(0) \\ 0 \\ \vdots \\ 0 \\ g(1) \end{bmatrix}

The relevant sparse matrix here is AA, given by

A=[−211−21⋱⋱⋱1−211−2] A = \begin{bmatrix} -2 & 1 & & & \\ 1 & -2 & 1 & & \\ &\ddots & \ddots & \ddots &\\ & & 1 & -2 & 1 \\ & & & 1 & -2 \end{bmatrix}

As you can see, the number of non-zero elements grows linearly with the size NN, so a sparse matrix format is much preferred over a dense matrix holding all N2N^2 elements!

Additional Reading

For more on the Finite Difference Method for solving PDEs:

K. W. Morton and D. F. Mayers. Numerical Solution of Partial Differential Equations: An Introduction. Cambridge University Press, 2005.