Solving 2D Poisson Eq with mixed BC's in Python

Question

I am trying to numerically solve the Poisson's equation

$$ u_{xx} + u_{yy} = - \cos(x) \quad \text{if} - \pi/2 \leq x \leq \pi/2 \quad \text{0 otherwise} $$

The domain is the rectangle with vertices $(-π, 0), (-π,2), (π,0)$ and $(π,2)$. The boundary conditions are mixed:

Dirichlet: $u(x,0) = 0$ Neumann: $u_y(x,2) = 0$ Periodic: $u(-\pi,0) = u(\pi,0) = 0$.

I have the code written, but I can't figure out why my code isn't plotting the correct solution.

My solution:

https://i.stack.imgur.com/i2ced.jpg

Correct solution:

https://capture.dropbox.com/4hzsLwatNsc05Fhl

Here is my code:

import numpy as np
from scipy.sparse import diags, csc_matrix
from scipy.sparse.linalg import spsolve
import matplotlib.pyplot as plt
import seaborn as sns
import math

x_min = -np.pi #Left endpoint of x interval
x_max = np.pi  #Right endpoint of x interval
y_min = 0      #Left endpoint of y interval
y_max = 2      #Right endpoint of y interval

nx = 50                                     # Number of grid points in x
ny = 50                                     # Number of grid points in y
dx = (x_max - x_min)/nx                    # Spacing in x direction
dy = (y_max - y_min)/ny                    # Spacing in y direction
n = nx*ny                                  # Dimension of system
x = np.linspace(x_min, x_max-dx, nx)       # Grid points in x-direction
y = np.linspace(y_min+dy, y_max, ny)       # Grid points in y-direction
xg, yg = np.meshgrid(x, y)     

# Create the diagonals
ones = np.full(n, 1)

a=np.ones(ny-1,dtype='float')
b=np.array([0.0],dtype='float')
c=np.concatenate((a,b),axis=0)
upper_diagonal = np.tile(c,nx)

d=np.ones(ny-2,dtype='float')
e=np.array([2.0,0.0],dtype='float')
f=np.concatenate((d,e),axis=0)
lower_diagonal = np.tile(f,nx)

# Create the offsets
offsets = [0, -1, 1, ny, -ny, (nx-1)*ny, -(nx-1)*ny]

# Create the sparse matrix
A = diags([-2*ones/dx**2 -2*ones/dy**2,lower_diagonal/dx**2,upper_diagonal/dx**2,ones/dy**2,ones/dy**2,ones/dy**2,ones/dy**2], offsets, shape=(n, n), dtype=float)

# Construct the RHS
b = np.zeros((ny,nx),dtype=float)

for j in range(nx):
    
    #for i in range(ny):
    
    if abs(x[j]) <= np.pi/2:    

        b[:,j] = -np.cos( x[j] )

b1 = b.reshape(n,1)

#Solve the linear system using a sparse matrix solver
As = csc_matrix(A)
bs = csc_matrix(b1)
u = spsolve(As, bs).reshape(ny, nx)

#Plot the solution
fig = plt.figure(figsize=(8, 8))
ax = plt.axes(projection='3d')
surf = ax.plot_surface(xg, yg, u)
ax.set_xlabel('x')
ax.set_ylabel('y')
ax.set_zlabel('u');
plt.show()

The only mistake I can think of is in the construction of the matrix $A$. Because I have been able to write a similar code where I can solve a number of BVPs with all Dirichlet conditions. So I must not be taking into account the Neumann and periodic boundary conditions correctly into account. However, I think the matrices I generate are correct. For instance, here's the heatmap of the matrix when $n_y = 7, n_x = 4$:

https://i.stack.imgur.com/NdHit.jpg

For simplicity, I have shown the matrix prior to having divided each entry by $\delta x, \delta y$. The distribution of $-4, 1, 2$'s seems right to me, so I don't know what's going wrong. Suggestions?

Answer 1

Poisson's equation:

$$u_{xx} + u_{yy} = - \cos(x) \quad \text{if} - \pi/2 \leq x \leq \pi/2 \quad \text{0 otherwise}$$

Boundary conditions:

Dirichlet: $u(x,0)= u(-\pi,y)= u(\pi,y)=0$ Neumann: $u_y(x,2) = 0$

2D discretization of the Laplace operator: $$u_{xx} + u_{yy}\simeq\frac{u(x-h_x,y)+u(x+h_x,y)-2 u(x,y)}{h_x^2}+\frac{u(x,y-h_y)+u(x,y+h_y)-2 u(x,y)}{h_y^2}$$

Neumann BC with Backward difference approximation: $$u_y(x,y)\simeq \frac{-4 u(x,y-h_y)+u(x,y-2 h_y)+3 u(x,y)}{2 h_y}$$

Here is the python code:

import numpy as np
from scipy.sparse import diags, csc_matrix
from scipy.sparse.linalg import spsolve
import matplotlib.pyplot as plt
import seaborn as sns
import math

# Setup
x_min = -np.pi # Left endpoint of x interval
x_max = +np.pi # Right endpoint of x interval
y_min = 0      # Left endpoint of y interval
y_max = 2      # Right endpoint of y interval

nx = 31                                    # Number of grid points in x
ny = 21                                    # Number of grid points in y
hx = (x_max - x_min)/(nx-1)                # Spacing in x direction
hy = (y_max - y_min)/(ny-1)                # Spacing in y direction
n  = nx*ny                                 # Dimension of system

A = np.zeros((n, n), dtype='float')
b = np.zeros(n, dtype='float')

# Define Lambda functions
index = lambda row, col: row*ny+col
f = lambda x: np.cos(x) if np.abs(x)<=np.pi/2 else 0

# Index of equation
row = 0

# Set 2D Laplace operator
for i in range(1, nx-1):
     for j in range(1, ny-1):
         A[row, index(i-1, j)] += 1/hx**2
         A[row, index(i, j)]   -= 2/hx**2+2/hy**2
         A[row, index(i+1, j)] += 1/hx**2
         A[row, index(i, j-1)] += 1/hy**2
         A[row, index(i, j+1)] += 1/hy**2
         b[row] = -f(x_min+hx*(i-1))
         row += 1

# Set Neumann boundary condition
for i in range(1, nx-1):
    if row >= n: break
    A[row, index(i, ny-1)] += 3 #/(2*hy)
    A[row, index(i, ny-2)] -= 4 #/(2*hy)
    A[row, index(i, ny-3)] += 1 #/(2*hy)
    b[row] = 0
    row += 1

# Set Dirichlet boundary conditions
for i in range(0, nx):
    for j in range(0, ny):
        if i<=0 or j<=0 or i>=nx-1:
            if row >= n: break
            A[row, index(i, j)] = 1
            b[row] = 0
            row += 1
        
if row != n:
     raise Exception('Wrong number of equations')

plt.spy(A)

rank = np.linalg.matrix_rank(A)
print('rank=',rank,'row=', row)

if rank < row:
    raise Exception('Insufficient rank.')

# Solve the linear system
u = np.linalg.solve(A, b)
# np.savetxt('uuu.csv', u, delimiter=',')

error = np.matmul(A,u) - b
norm = np.linalg.norm(error)
print('error norm=', norm)

if norm > 1E-9:
    raise Exception('Error norm')

# Prepare visualization
x = np.linspace(x_min, x_max, nx)       # Grid points in x-direction
y = np.linspace(y_min, y_max, ny)       # Grid points in y-direction
xg, yg = np.meshgrid(x, y)     
uu = u.reshape(nx, ny).transpose()

# Plot the solution
fig = plt.figure(figsize=(8, 8))
ax = plt.axes(projection='3d')
surf = ax.plot_surface(xg, yg, uu)
ax.set_xlabel('x')
ax.set_ylabel('y')
ax.set_zlabel('u');
ax.azim = 45
plt.show()