Project 10: Numerical Solutions of the 2D Elliptic Equation

Objective

To numerically solve the 2D Piosson equation using finite difference methods and analyze the results.

Problem Description

The 2D Laplacian equation, often appearing in contexts like heat conduction, electrostatics, and fluid dynamics, is given by:

- Δ u = f (x, y)

where $Δ u = \nabla^{2} u = \frac{\partial^{2} u}{\partial x^{2}} + \frac{\partial^{2} u}{\partial y^{2}}$ , $f (x, y)$ is a given function and when $f (x, y) = 0$ , it is the Laplacian equation.

Domain: Let’s consider a rectangular domain $Ω$ defined by the coordinates $(x, y)$ where $0 \leq x \leq 1$ and $0 \leq y \leq 1$ .
Boundary condition: we use Dirichlet boundary condition for this problem
1. Left Edge (x = 0): Let’s set $f_{1} (y)$ to be a linear function of $y$ . $f_{1} (y) = 2 \sin (2 π y)$
2. Right Edge (x = 1): For $f_{2} (y)$ , let’s use a constant function. $f_{2} (y) = 2 \sin (2 π y)$
3. Bottom Edge (y = 0): Let $g_{1} (x)$ be a sinusoidal function. $g_{1} (x) = \sin (2 π x)$
4. Top Edge (y = 1): Finally, for $g_{2} (x)$ , we can choose a quadratic function. $g_{2} (x) = \sin (2 π x)$

Step 1: Discretization and initialization

Discretize the domain $Ω$ into a grid of $(N + 1) \times (N + 1)$ points, where $N_{x}$ and $N_{y}$ are the number of grid points in the $x$ and $y$ directions, respectively. The grid spacing in the $x$ and $y$ directions are $Δ x = 1 / N$ and $Δ y = 1 / N$ , respectively. The grid points are indexed by $i$ and $j$ in the $x$ and $y$ directions, respectively, as shown in the figure below.

import numpy as np
N = 10
x = np.linspace(0,1,N+1) # N+1 points in x direction
y = np.linspace(0,1,N+1)
dx = x[1]-x[0]
dy = y[1]-y[0]

Step 3: Boundary Condition Function

Define the boundary condition functions $f_{1} (y)$ , $f_{2} (y)$ , $g_{1} (x)$ , and $g_{2} (x)$ .

def f1(y):
    return 
def f2(y):
    return 
def g1(x):
    return 
def g2(x):
    return 

Step 4: Numerical Solution

Use the finite difference method to discretize the 2D Laplacian equation and obtain the numerical solution $u (x, y)$ at $x = x_{i}$ and $y = y_{j}$ for $i = 1, 2, \dots, N - 1$ and $j = 1, 2, \dots, N - 1$ .

The Poisson Equation is discretized using $δ_{x}^{2}$ is the central difference approximation of the second derivative in the $x$ direction $\begin{matrix} (1) & δ_{x}^{2} = \frac{1}{h^{2}} (u_{i + 1 j} - 2 u_{i j} + u_{i - 1 j}), \end{matrix}$ and $δ_{y}^{2}$ is the central difference approximation of the second derivative in the $y$ direction $\begin{matrix} (2) & δ_{y}^{2} = \frac{1}{h^{2}} (u_{i j + 1} - 2 u_{i j} + u_{i j - 1}) . \end{matrix}$ This gives the Poisson Difference Equation, $\begin{matrix} (3) & - (δ_{x}^{2} u_{i j} + δ_{y}^{2} u_{i j}) = f_{i j} (x_{i}, y_{j}) \in Ω_{h}, \end{matrix}$ $\begin{matrix} (4) & u_{i j} = g_{i j} (x_{i}, y_{j}) \in \partial Ω_{h}, \end{matrix}$ where $u_{i} j$ is the numerical approximation of $U$ at $x_{i}$ and $y_{j}$ . Expanding the Poisson Difference Equation gives the five-point method, $\begin{matrix} (5) & - (u_{i - 1 j} + u_{i j - 1} - 4 u_{i j} + u_{i j + 1} + u_{i + 1 j}) = h^{2} f_{i j} \end{matrix}$ for $i = 1, . . ., N - 1$ and $j = 1, . . ., N - 1$ which can be written $\nabla_{h}^{2} u_{i j} = f_{i j}$

The understanding of this five-point method is important to construct the matrix $A$ and the right-hand side vector $b$ in the linear system $\begin{matrix} (6) & A u = b \end{matrix}$ $u$ is the vector of unknowns $u_{i j}$ , $A$ is the matrix of coefficients of the unknowns, and $b$ is the right-hand side vector.

This can be written as a system of $(N - 1) \times (N - 1)$ equations can be arranged in matrix form $\begin{matrix} (7) & A u = r, \end{matrix}$ , where $A$ is an $(N - 1)^{2} \times (N - 1)^{2}$ matrix made up of the following block tridiagonal structure

(\begin{array}{ccccccc} T & I & 0 & 0 & . & . & . \\ I & T & I & 0 & 0 & . & . \\ . & . & . & . & . & . & . \\ . & . & . & 0 & I & T & I \\ . & . & . & . & 0 & I & T \end{array}),

where $I$ denotes an $(N - 1) \times (N - 1)$ identity matrix and $T$ is the tridiagonal matrix of the form:

T = (\begin{array}{ccccccc} - 4 & 1 & 0 & 0 & . & . & . \\ 1 & - 4 & 1 & 0 & 0 & . & . \\ . & . & . & . & . & . & . \\ . & . & . & 0 & 1 & - 4 & 1 \\ . & . & . & . & 0 & 1 & - 4 \end{array}) .

N2=(N-1)*(N-1)
A=np.zeros((N2,N2))
## Diagonal            
for i in range (0,N-1):
    for j in range (0,N-1):           
        A[i+(N-1)*j,i+(N-1)*j]=-4

# LOWER DIAGONAL        
for i in range (1,N-1):
    for j in range (0,N-1):           
        A[i+(N-1)*j,i+(N-1)*j-1]=1   
# UPPPER DIAGONAL        
for i in range (0,N-2):
    for j in range (0,N-1):           
        A[i+(N-1)*j,i+(N-1)*j+1]=1   

# LOWER IDENTITY MATRIX
for i in range (0,N-1):
    for j in range (1,N-1):           
        A[i+(N-1)*j,i+(N-1)*(j-1)]=1        
        
        
# UPPER IDENTITY MATRIX
for i in range (0,N-1):
    for j in range (0,N-2):           
        A[i+(N-1)*j,i+(N-1)*(j+1)]=1

After we assemble the matrix $A$ in Eq. ( $6$ ), we need to construct the right-hand side vector $b$ in Eq. ( $6$ ). The right-hand side vector $b$ is a vector of length $(N - 1)^{2}$ , where $b_{i}$ is the right-hand side of the $i$ -th equation in Eq. ( $6$ ). The right-hand side vector $b$ is constructed by evaluating the function $f (x, y)$ at the grid points $(x_{i}, y_{j})$ for $i = 1, 2, \dots, N - 1$ and $j = 1, 2, \dots, N - 1$ .

# Construct the right-hand side vector b 
b = np.zeros(N2)

Step 5: Plot the Numerical Solution

To solve the system for $u$ invert the matrix $A$ $\begin{matrix} (8) & A u = r, \end{matrix}$ such that $\begin{matrix} (9) & u = A^{- 1} r . \end{matrix}$ Lastly, as $u$ is in vector it has to be reshaped into grid form to plot.

C=np.dot(np.linalg.inv(A),r-b)
w[1:N,1:N]=C.reshape((N-1,N-1))

fig = plt.figure(figsize=(8,6))
ax = fig.add_subplot(111, projection='3d');
# Plot a basic wireframe.
ax.plot_wireframe(X, Y, w,color='r');
ax.set_xlabel('x');
ax.set_ylabel('y');
ax.set_zlabel('w');
plt.title(r'Numerical Approximation of the Poisson Equation',fontsize=24,y=1.08);
plt.show();

Step 6: Analyze the Numerical Solution

We only study a simple boundary case, where $f_{1} (y) = f_{2} (y) = g_{1} (x) = g_{2} (x) = 0$ . In this case, the exact solution is $u (x, y) = 0$ . We can compute the error between the numerical solution and the exact solution.

error = np.linalg.norm(u)
print('Error = ', error)

Deliverables

Have a jupyter-notebook to show your code and results. Based on the results, prepare a 4-minute presentation to show your results and answer the question.

Reference: Finite Difference Methods for the Poisson Equation