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:
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$:
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?
