Laplacian matrix in two space dimensions¶
AMath 585, Winter Quarter 2020 at the University of Washington. Developed by R.J. LeVeque and distributed under the BSD license. You are free to modify and use as you please, with attribution.
These notebooks are all available on Github.
This notebook illustrates how to convert a 2D grid of points to a vector as needed in solving a linear system, and then how to set up the matrix corresponding to the 2D 5-point Laplacian using Kronecker products.
First see Grids2D.html for the discussion on setting up grid functions in 2D and plotting them.
%matplotlib inline
from pylab import *
from scipy import sparse
Make a little function to plot grids:
def plot_grid(X,Y):
plot(X,Y,'k')
plot(X.T,Y.T,'k')
#plot(X[1:-1,1:-1],Y[1:-1,1:-1],'ro')
axis('scaled')
Define a sample 2D function:¶
Using code from Grids2D.html:
mx = 9
x = linspace(0,1,mx+2)
print('x = ',x)
my = 3
y = linspace(5,6,my+2)
print('y = ',y)
X,Y = meshgrid(x,y,indexing='ij')
print('X = \n', X)
print('Y = \n', Y)
print()
U = X + 10*Y
print('U = \n', U)
Converting between 2D grid functions and vectors¶
We defined U as a 2D array but to solve a linear system for the elements of U we need to reshape it into a 1D vector. One can do this in the natural row-wise or column-wise order, or some other order. The solution of the linear system should be the same in any case (assuming the right hand side is also reshaped in the same way), but the structure of the matrix will change.
The reshape command will reshape arrays, by default in the column-wise ordering:
xvec_colwise = reshape(X, (mx+2)*(my+2))
print('xvec_colwise = \n', xvec_colwise)
yvec_colwise = reshape(Y, (mx+2)*(my+2))
print('yvec_colwise = \n', yvec_colwise)
For row-wise ordering, use order='F'. This is "Fortran ordering" while the default is "C ordering", i.e. if you don't specify order the default is order='C'. The naming comes from the way arrays are laid out in memory when programming in the Fortran or C languages.
xvec_rowwise = reshape(X, (mx+2)*(my+2), order='F')
print('xvec_rowwise = \n', xvec_rowwise)
yvec_rowwise = reshape(Y, (mx+2)*(my+2), order='F')
print('yvec_rowwise = \n', yvec_rowwise)
uvec_rowwise = reshape(U, (mx+2)*(my+2), order='F')
print('uvec_rowwise = \n', uvec_rowwise)
To go back from a vector to a 2D array, we can again use reshape (and make sure to specify the same order as used in the other direction).
UU = reshape(uvec_rowwise, (mx+2, my+2), order='F')
print('UU = \n', UU)
print('Max difference between U and UU: ', abs(U-UU).max())
Kronecker product matrices¶
If grid arrays are converted to vectors using the natural row-wise or column-wise order, then finite difference operators such as the discrete Laplacian can be expressed as matrices that are the Kronecker product of 1-dimensional discretization matrices.
First here's a small example of Kronecker product matrices:
A = array([[1,2],[3,4.]])
B = diag((1,10,100))
print('A = \n',A)
print('B = \n',B)
print('kron(A,B) = \n',kron(A,B))
print('kron(B,A) = \n',kron(B,A))
Matrix for 2D Laplacian¶
Suppose we use the natural row-wise ordering. For clarity assume $h = \Delta x = \Delta y = 1$ so that $1/h^2 = 1$ and these factors don't appear (but more generally remember to include them).
Then in the case mx=3 and my=2, for example, we want to construct a matrix that is $2\times 2$ block tridiagonal where each block is a $3\times 3$ matrix:
where
\begin{align*} D_x = \left[ \begin{array}{ccc|ccc} -2 & 1 & 0 & 0 & 0 & 0 \\ 1 &-2 & 1 & 0 & 0 & 0 \\ 0 & 1 &-2 & 0 & 0 & 0 \\ \hline 0 & 0 & 0 &-2 & 1 & 0 \\ 0 & 0 & 0 & 1 &-2 & 1 \\ 0 & 0 & 0 & 0 & 1 &-2 \\ \end{array} \right] \quad\text{and}\quad D_y = \left[ \begin{array}{ccc|ccc} -2 & 0 & 0 & 1 & 0 & 0 \\ 0 & -2 & 0 & 0 & 1 & 0 \\ 0 & 0 & -2 & 0 & 0 & 1 \\ \hline 1 & 0 & 0 & -2 & 0 & 0 \\ 0 & 1 & 0 & 0 & -2 & 0 \\ 0 & 0 & 1 & 0 & 0 & -2\\ \end{array} \right] \end{align*}These can be constructed as:
mx = 3
my = 2
em = ones(mx)
em1 = ones(mx-1)
Ax = sparse.diags([em1, -2*em, em1], [-1, 0, 1], shape=(mx,mx)).toarray()
Iy = eye(my)
Dx = kron(Iy,Ax)
print('Dx = \n',Dx)
Ix = eye(mx)
em1 = ones(my-1)
Ay = sparse.diags([em1, -2*em, em1], [-1, 0, 1], shape=(my,my)).toarray()
Dy = kron(Ay,Ix)
print('Dy = \n',Dy)
A = Dx + Dy
print('A = \n', A)
For a slightly larger case:
mx = 4
my = 3
em = ones(mx)
em1 = ones(mx-1)
Ax = sparse.diags([em1, -2*em, em1], [-1, 0, 1], shape=(mx,mx)).toarray()
Iy = eye(my)
Dx = kron(Iy,Ax)
print('Dx = \n',Dx)
Ix = eye(mx)
em1 = ones(my-1)
Ay = sparse.diags([em1, -2*em, em1], [-1, 0, 1], shape=(my,my)).toarray()
Dy = kron(Ay,Ix)
print('Dy = \n',Dy)
A = Dx + Dy
print('A = \n', A)
Function to make the discrete Laplacian matrix¶
Now we include the 1/dx**2 and 1/dy**2 terms.
def make_A(mx,my,dx,dy):
em = ones(mx)
em1 = ones(mx-1)
Ax = sparse.diags([em1, -2*em, em1], [-1, 0, 1], shape=(mx,mx)).toarray()
Iy = eye(my)
Dx = kron(Iy,Ax) / dx**2
Ix = eye(mx)
em1 = ones(my-1)
Ay = sparse.diags([em1, -2*em, em1], [-1, 0, 1], shape=(my,my)).toarray()
Dy = kron(Ay,Ix) / dy**2
A = Dx + Dy
return A
Test by applying to a function:¶
The discrete Laplacian matrix as set up above can be applied only to the interior points, and effectively assumes $u \equiv 0$ on the boundary of the rectangle since these points are not included.
So we choose a test function $u(x,y)$ that is zero on the boundary and quadratic in $x$ and $y$ so that the 5-point Laplacian should be exact.
$$ u(x,y) = x(2-x)y(1-y) \quad\text{for}~0\leq x \leq 2, ~~0\leq y \leq 1. $$ax = 0.
bx = 2.
ay = 0.
by = 1.
mx = 19
my = 9
dx = (bx-ax)/(mx+1)
dy = (by-ay)/(my+1)
A = make_A(mx,my,dx,dy)
x = linspace(0,2,mx+2)
y = linspace(0,1,my+2)
X,Y = meshgrid(x,y,indexing='ij')
U = (X-ax)*(bx-X)*(Y-ay)*(by-Y)
# interior points:
xint = x[1:-1]
yint = y[1:-1]
Uint = U[1:-1, 1:-1]
Xint = X[1:-1, 1:-1]
Yint = Y[1:-1, 1:-1]
# reshape as vectors in row-wise ordering:
xvec_rowwise = reshape(Xint, mx*my, order='F')
yvec_rowwise = reshape(Yint, mx*my, order='F')
uvec_rowwise = reshape(Uint, mx*my, order='F')
print('shape of Uint: ', Uint.shape)
print('uvec_rowwise has length %i' % len(uvec_rowwise))
print('shape of A: ', A.shape)
plot_grid(X,Y)
contourf(X,Y,U)
colorbar(shrink=0.6)
title('Origina u(x,y) on full domain');
Apply the matrix and then reshape result as 2D array:¶
ulap_rowwise = dot(A,uvec_rowwise)
Ulap_int = reshape(ulap_rowwise, (mx,my), order='F')
Compare to exact Laplacian:¶
Ulap_true = -2*(X-ax)*(bx-X) - 2*(Y-ay)*(by-Y)
print('Maximum error is %.4e' % abs(Ulap_int- Ulap_true[1:-1,1:-1]).max())
Plot $\nabla^2 u(x,y)$¶
Note that this is defined only at the interior points so the contourf function colors only in the interior of this region:
plot_grid(X,Y)
contourf(Xint,Yint,Ulap_int)
colorbar(shrink=0.6)
axis([ax,bx,ay,by]);
Solving the Poisson problem¶
The matrix $A$ defined above could be used in solving a Poisson problem $\nabla^2u(x,y) = f(x,y)$.
This is not set up here -- it would require incorporating the boundary conditions into the right hand side along with the function values $f(x_i,y_j)$.