This program computes an approximate solution of Dirichlet's problem for the Poisson equation on the unit square using the 5-point Laplacian. Thus it finds a grid function \(u\) such that \[ \Delta_h u_h = f \text{ on interior grid points}, \quad u_h = g \text{ on boundary grid points.} \]
from numpy import pi, sin, cos, exp
import numpy as np
import scipy
from scipy.sparse.linalg import spsolve
Supply the functions \(f\) and \(u\). \(u\) need only be given on the boundary (i.e., \(g\)), but if it is known and given at all grid points we can compute the error below.
# Example 0: u is a polynomial of degree < 4 in x and y, so the method is exact at the nodes
#def u(x, y):
# return x**3*(y**3 - 3*y**2 + 2)/2 - 100*(x**2-x)*y
#def f(x, y):
# return 3*x*y**3 - 9*x*y**2 + 3*x**3*y -3*x**3 - 200*y + 6*x
# Example 1: u = sin pi x sin pi y
def u(x, y):
return sin(pi*x)*sin(pi*y)
def f(x, y):
return -2*pi**2*sin(pi*x)*sin(pi*y)
# Example 2: f = 0, u is known only on the boundary
#def u(x, y):
# return .02*sin(8*pi*x*y)
#def f(x, y):
# return np.zeros(x.shape)
Construct the mesh. Grid points are: \[ (X_{ij}, Y_{ij}) = (x_i, y_j), \quad 0\le i, j\le n. \] Interior grid points are: \[ (X^0_{ij}, Y^0_{ij}) = (X_{i+1,j+1}, Y_{i+1,j+1}),\quad 0\le i,j\le n-2. \]
n = 8
h = 1./n
x = np.linspace(0., 1., n+1)
y = np.linspace(0., 1., n+1)
[X, Y] = np.meshgrid(x, y, indexing='ij')
Evaluate \(f\) at interior grid points
fh0 = f(X[1:n, 1:n], Y[1:n, 1:n])
Evaluate \(u\) at grid points (values at interior grid points will be ignored).
uh = u(X, Y)
Form the 5-point Laplacian matrix (times \(h^2\)) as a sparse matrix using CSR (compressed sparse row) storage.
I = scipy.sparse.eye(n-1) # (n-1) x (n-1) identity matrix
e = np.ones(n-1) # vector (1, 1, ..., 1) of length n-1
e0 = np.ones(n-2) # vector (1, 1, ..., 1) of length n-2
A = scipy.sparse.diags([e0, -4*e, e0], [-1, 0, 1]) # (n-1) x (n-1) tridiagonal matrix
J = scipy.sparse.diags([e0, e0], [-1, 1]) # same with zeros on diagonal
Lh = scipy.sparse.kron(I, A) + scipy.sparse.kron(J, I) # the desired matrix
Lh = Lh.tocsr() # convert it to CSR
Create the right hand side by multiplying \(f\) by \(h^2\), adjusting for boundary values, and then reshaping the (n-1) x (n-1) matrix as an vector of length \((n-1)^2\).
fh0 = h**2*fh0
fh0[:, 0] -= uh[1:n, 0]
fh0[:, n-2] -= uh[1:n, n]
fh0[0, :] -= uh[0, 1:n]
fh0[n-2, :] -= uh[n, 1:n]
F = np.reshape(fh0, (n-1)**2)
Solve the linear system to find the vector \(U\) corresonding to the interior values of the grid function \(u_h\), reshape it into an \((n-1)\times(n-1)\) matrix, and insert it into the grid function \(u_h\).
U = spsolve(Lh, F)
uh[1:n, 1:n] = np.reshape(U, (n-1, n-1))
Plot the computed solution.
# plot the computed solution
from mpl_toolkits.mplot3d import Axes3D
from matplotlib import cm
from matplotlib.ticker import LinearLocator, FormatStrFormatter
import matplotlib.pyplot as plt
fig = plt.figure()
ax = fig.gca(projection='3d')
surf = ax.plot_surface(X, Y, uh, rstride=1, cstride=1, color='yellow')
plt.show()
If the true solution was given above, compute the error.
utrue = u(X, Y)
err = abs(uh - utrue).max()
print("Absolute error in max norm = {:10.3e}".format(err))
print("Relative error in max norm = {:6.3f}%".format(100*err/abs(utrue).max()))