Consider the boundary value problem

\begin{equation} \begin{aligned} u_{xx} + u_{yy} &= -\sin(\pi x) \sin(\pi y) \quad \text{for all } 0 < x < 2 \text{ and } 0 < y < 2 \\ u(x,0) &= 0 = u(x,2) \quad \text{for all } 0 \le x \le 2 \\ u(0,y) &= 0 = u(2,y) \quad \text{for all } 0 \le y \le 2 \end{aligned} \end{equation}

The exact solution to this BVP is

\begin{equation} u(x,y) = \frac{1}{2\pi^2} \sin(\pi x) \sin(\pi y) \end{equation}

(a) Solve this problem using the direct method (approximate the Laplacian with a 5-point stencil and solve the resulting system of equations)

\begin{equation} u_{xx} + u_{yy} \approx \frac{1}{\Delta x^2} \left( u_{i+1,j} + u_{i-1,j} + u_{i,j+1} + u_{i,j-1} - 4u_{i,j} \right) \end{equation}

where $\Delta x = \Delta y$ Use several different values of $\Delta x$ and make a log-log plot of the error to determine the order of accuracy.

#problem 4 a)
#Here I am just plotting it to see what the True solutions and the numerical solution looks like:
import numpy as np
import matplotlib.pyplot as plt
import scipy

#%matplotlib notebook

x0 = 0
xN = 2
y0 = 0
yN = 2

dx = 0.01
dy = dx

x = np.arange(x0, xN + dx/2, dx)
y = np.arange(y0, yN + dy/2, dy)

Nx = len(x)
Ny = len(y)
#boundary conditions
#u(x, -1)
def bot(x):
    return 0*x

#u(x,1)
def top(x):
    return 0*x

#u(-1,y)
def left(y):
    return 0*y

#u(1,y)
def right(y):
    return 0*y

def f(x,y):
    return -np.sin(np.pi*x)*np.sin(np.pi*y)

def true_sol(x, y):
    return 1/(2*np.pi**2)*np.sin(np.pi*x)*np.sin(np.pi*y)

N_total = (Nx - 2) * (Ny - 2)

A = scipy.sparse.dok_array((N_total, N_total))# make a sparse matrix instead
b = np.zeros((N_total, 1))

def point2ind(m, n):
    return (n - 1) * (Nx - 2) + m - 1

#filling matrix A and b, b needs subtract the boundaries if u(m,n) is next to a boundary and include the right side of the pde
#u_(m,n)
for n in range(1, Ny-1): #loop for n and m is slightly different because we got rid of boundary points
    for m in range(1, Nx-1): 
        k = point2ind(m, n)
        A[k, k] = -2 / dx ** 2 - 2/dy**2
        if m > 1:# good
            A[k, k - 1] = 1 / dx ** 2
        if n < Ny - 2: #good
            A[k, k + Nx - 2] = 1 / dy ** 2
        if m < Nx - 2:
            A[k, k + 1] = 1 / dx ** 2
        if n > 1:
            A[k, k - (Nx - 2)] = 1 / dy ** 2
        if n == 1:
            b[k] = b[k] - bot(x[m]) / dy ** 2
        if n == Ny-2:
            b[k] = b[k] - top(x[m]) / dy ** 2
        if m == 1:
            b[k] = b[k] - left(y[n]) / dx ** 2
        if m == Nx-2:
            b[k] = b[k] - right(y[n]) / dx ** 2
        
        b[k] = b[k] + f(x[m], y[n])
            
            
u = scipy.sparse.linalg.spsolve(A, b)

U_int = u.reshape((Ny-2, Nx-2))
U = np.zeros((Ny, Nx))
U[1:(Ny-1), 1:(Nx-1)] = U_int
U[0, :] = bot(x) #lower boundary
U[-1,:] = top(x)
U[:,0] = left(y)
U[:,-1] = right(y)


#These vectors are just to plot the boundary conditions
minusonex = np.ones_like(x)*0
onex = np.ones_like(x)*2
minusoney = np.ones_like(y)*0
oney = np.ones_like(y)*2

X, Y = np.meshgrid(x, y)

U_true = true_sol(X,Y)

ax = plt.axes(projection='3d')
# Plot the solution
ax.plot_surface(X, Y, U, alpha = 0.7, cmap='viridis')
# u(-1, y) 
ax.plot3D(minusoney, y, left(y), 'r')
# u(x, 1) 
ax.plot3D(x, onex, top(x), 'r')
# u(1, y) 
ax.plot3D(oney, y, right(y), 'r')
# u(x, -1) 
ax.plot3D(x, minusonex, bot(x), 'r')

ax.set_xlabel('x')
ax.set_ylabel('y')

plt.title('5-Point Laplacian Approximation')

plt.show()

ax = plt.axes(projection='3d')

# Plot the true solution
ax.plot_surface(X, Y, U_true, alpha = 0.7, cmap='viridis')
# u(-1, y) 
ax.plot3D(minusoney, y, left(y), 'r')
# u(x, 1) 
ax.plot3D(x, onex, top(x), 'r')
# u(1, y) 
ax.plot3D(oney, y, right(y), 'r')
# u(x, -1) 
ax.plot3D(x, minusonex, bot(x), 'r')

ax.set_xlabel('x')
ax.set_ylabel('y')

plt.title('True Solution')

plt.show()

#The plots above look alright!

#The plot below is what the solution of the pde looks like.

C:\Users\samip\anaconda3\Lib\site-packages\scipy\sparse\linalg\_dsolve\linsolve.py:214: SparseEfficiencyWarning: spsolve requires A be CSC or CSR matrix format
  warn('spsolve requires A be CSC or CSR matrix format',

Above is the numerical solution with $\Delta x = 0.01$ along side the true solution. By inspection, they look identical, so I can be confident my implementation of the numerical solution is correct.

Below, I will use the 5-point stencil numerical solution with values of $\Delta x$ of $\{2^{-5}, 2^{-6}, 2^{-7}, 2^{-8}, 2^{-9}\}$. For each value of $\Delta x$, I will find the maximum error between the true solution and the numerical solution. The errors will be plotted using a log-log plot to determine the order of accuracy of this numerical method.

#Here, I will copy the code from the previous cell and calculate
#the global error for various values of dx
import numpy as np
import matplotlib.pyplot as plt
import scipy
#boundary conditions
#u(x, -1)
def bot(x):
    return 0*x

#u(x,1)
def top(x):
    return 0*x

#u(-1,y)
def left(y):
    return 0*y

#u(1,y)
def right(y):
    return 0*y

def f(x,y):
    return -np.sin(np.pi*x)*np.sin(np.pi*y)

def true_sol(x, y):
    return 1/(2*np.pi**2)*np.sin(np.pi*x)*np.sin(np.pi*y)

def point2ind(m, n):
    return (n - 1) * (Nx - 2) + m - 1

dx_list = np.array([2**-5, 2**-6, 2**-7, 2**-8, 2**-9])
error_list5 = np.array([])

x0 = 0
xN = 2
y0 = 0
yN = 2

#loop over different values of dx
for dx in dx_list:
    dx = dx
    dy = dx

    x = np.arange(x0, xN + dx/2, dx)
    y = np.arange(y0, yN + dy/2, dy)

    Nx = len(x)
    Ny = len(y)

    N_total = (Nx - 2) * (Ny - 2)

    A = scipy.sparse.dok_array((N_total, N_total))# make a sparse matrix instead
    b = np.zeros((N_total, 1))


    #filling matrix A and b, b needs subtract the boundaries if u(m,n) is next to a boundary and include the right side of the pde
    #u_(m,n)
    for n in range(1, Ny-1): #loop for n and m is slightly different because we got rid of boundary points
        for m in range(1, Nx-1): 
            k = point2ind(m, n)
            A[k, k] = -2 / dx ** 2 - 2/dy**2
            if m > 1:# good
                A[k, k - 1] = 1 / dx ** 2
            if n < Ny - 2: #good
                A[k, k + Nx - 2] = 1 / dy ** 2
            if m < Nx - 2:
                A[k, k + 1] = 1 / dx ** 2
            if n > 1:
                A[k, k - (Nx - 2)] = 1 / dy ** 2
            if n == 1:
                b[k] = b[k] - bot(x[m]) / dy ** 2
            if n == Ny-2:
                b[k] = b[k] - top(x[m]) / dy ** 2
            if m == 1:
                b[k] = b[k] - left(y[n]) / dx ** 2
            if m == Nx-2:
                b[k] = b[k] - right(y[n]) / dx ** 2

            b[k] = b[k] + f(x[m], y[n])


    u = scipy.sparse.linalg.spsolve(A, b)

    U_int = u.reshape((Ny-2, Nx-2))
    U = np.zeros((Ny, Nx))
    U[1:(Ny-1), 1:(Nx-1)] = U_int
    U[0, :] = bot(x) #lower boundary
    U[-1,:] = top(x)
    U[:,0] = left(y)
    U[:,-1] = right(y)

    X, Y = np.meshgrid(x, y)
    U_true = true_sol(X,Y)

    error = np.max(np.abs(U-U_true))
    
    error_list5 = np.append(error_list5, error)
    
#print(error_list)

C:\Users\samip\anaconda3\lib\site-packages\scipy\sparse\linalg\_dsolve\linsolve.py:229: SparseEfficiencyWarning: spsolve requires A be CSC or CSR matrix format
  warn('spsolve requires A be CSC or CSR matrix format',

import pandas as pd
# Create DataFrame
df = pd.DataFrame({
    'Δx': dx_list,
    'Error': error_list5
})

# Round to 3 significant digits
formatted_df = df.applymap(lambda x: f'{x:.3e}')

# Print the DataFrame without the index
print(formatted_df.to_string(index=False))

# Plotting
plt.loglog(dx_list, error_list5, 'ro')
plt.title(r'Error vs $\Delta x$ Plot using 5-Point Laplacian')
plt.xlabel(r'$\Delta x$')
plt.ylabel('Max Error')
plt.show()

       Δx     Error
3.125e-02 4.071e-05
1.562e-02 1.017e-05
7.812e-03 2.543e-06
3.906e-03 6.358e-07
1.953e-03 1.589e-07

The printed values above the plot show a consistent pattern. Every time I a halve $\Delta x$, the error decreases by 1/4. So the error of the 5-Point Laplacian numerical solution is $O(\Delta x^2)$.

(b) Now solve this BVP using the following 9-point Laplacian:

\begin{equation} u_{xx} + u_{yy} \approx \nabla^2_9u = \frac{1}{6\Delta x^2} \left( 4u_{i-1,j} + 4u_{i+1,j} + 4u_{i,j-1} + 4u_{i,j+1} + u_{i-1,j-1} + u_{i-1,j+1} + u_{i+1,j-1} + u_{i+1,j+1} - 20u_{i,j} \right) \end{equation}

Use several different values of $\Delta x$ to solve the boundary value problem, then make a log-log plot of the errors to determine the order of accuracy.

#4 b) 9 point Laplacian
import numpy as np
import matplotlib.pyplot as plt
import scipy
#boundary conditions
#u(x, -1)
def bot(x):
    return 0*x

#u(x,1)
def top(x):
    return 0*x

#u(-1,y)
def left(y):
    return 0*y

#u(1,y)
def right(y):
    return 0*y

def f(x,y):
    return -np.sin(np.pi*x)*np.sin(np.pi*y)

def true_sol(x, y):
    return 1/(2*np.pi**2)*np.sin(np.pi*x)*np.sin(np.pi*y)

def point2ind(m, n):
    return (n - 1) * (Nx - 2) + m - 1

dx_list = np.array([2**-5, 2**-6, 2**-7, 2**-8, 2**-9])
error_list9 = np.array([])

x0 = 0
xN = 2
y0 = 0
yN = 2

#loop over different values of dx
for dx in dx_list:
    dx = dx
    dy = dx

    x = np.arange(x0, xN + dx/2, dx)
    y = np.arange(y0, yN + dy/2, dy)

    Nx = len(x)
    Ny = len(y)

    N_total = (Nx - 2) * (Ny - 2)

    A = scipy.sparse.dok_array((N_total, N_total))# make a sparse matrix instead
    b = np.zeros((N_total, 1))


    #filling matrix A and b, b needs subtract the boundaries if u(m,n) is next to a boundary and include the right side of the pde
    #u_(m,n)
    for n in range(1, Ny-1): #loop for n and m is slightly different because we got rid of boundary points
        for m in range(1, Nx-1): 
            k = point2ind(m, n)
            A[k, k] = -20 /(6* dx ** 2)
            if m > 1:# good
                A[k, k - 1] = 4 /(6*dx**2)
            if m < Nx - 2:
                A[k, k + 1] = 4 /(6*dx ** 2)
            if n > 1:
                A[k, k - (Nx - 2)] = 4 /(6*dx ** 2)
            if n < Ny - 2: #good
                A[k, k + Nx - 2] = 4 /(6*dx ** 2)
                
        
            if m > 1 and n > 1:
                A[k , k - 1 - (Nx - 2)] = 1/(6*dx**2)
            if m > 1 and n < Ny - 2:
                A[k , k - 1 + Nx - 2] = 1/(6*dx**2)
            if m < Nx - 2 and n > 1:
                A[k, k + 1 - (Nx - 2)] = 1/(6*dx**2)
            if m < Nx - 2 and n < Ny - 2:
                A[k, k + 1 + Nx - 2] = 1/(6*dx**2)
            
            
            #the boundaries are 0, so it doesn't matter about subtracting them 
            #to the right hand side. All of these are 0.
#             if n == 1:
#                 b[k] = b[k] - bot(x[m]) / dy ** 2
#             if n == Ny-2:
#                 b[k] = b[k] - top(x[m]) / dy ** 2
#             if m == 1:
#                 b[k] = b[k] - left(y[n]) / dx ** 2
#             if m == Nx-2:
#                 b[k] = b[k] - right(y[n]) / dx ** 2
            b[k] = b[k] + f(x[m], y[n])


    u = scipy.sparse.linalg.spsolve(A, b)

    U_int = u.reshape((Ny-2, Nx-2))
    U = np.zeros((Ny, Nx))
    U[1:(Ny-1), 1:(Nx-1)] = U_int
    U[0, :] = bot(x) #lower boundary
    U[-1,:] = top(x)
    U[:,0] = left(y)
    U[:,-1] = right(y)

    X, Y = np.meshgrid(x, y)
    U_true = true_sol(X,Y)

    error = np.max(np.abs(U-U_true))
    
    error_list9 = np.append(error_list9, error)

C:\Users\samip\anaconda3\lib\site-packages\scipy\sparse\linalg\_dsolve\linsolve.py:229: SparseEfficiencyWarning: spsolve requires A be CSC or CSR matrix format
  warn('spsolve requires A be CSC or CSR matrix format',

import pandas as pd

df = pd.DataFrame({
    'Δx': dx_list,
    'Error': error_list9
})

# Round to 3 significant digits
formatted_df = df.applymap(lambda x: f'{x:.3e}')

# Print the DataFrame without the index
print(formatted_df.to_string(index=False))

# Plotting
plt.loglog(dx_list, error_list9, 'ro')
plt.title(r'Error vs $\Delta x$ Plot using 9-Point Laplacian')
plt.xlabel(r'$\Delta x$')
plt.ylabel('Max Error')
plt.show()

       Δx     Error
3.125e-02 8.143e-05
1.562e-02 2.035e-05
7.812e-03 5.086e-06
3.906e-03 1.272e-06
1.953e-03 3.179e-07

As seen as above, the 9-point Laplacian for our given problem produces an error that is $O(\Delta x^2)$. This is actually what we expect for our Poisson Equation BVP. This error can be found by expanding the 9-Point stencil equation Taylor's series, revealing that the $O(\Delta x^2)$ comes from the right hand side of our BVP. If the right hand side was $0$, as in the Laplace equation, then this method would produce an error that is $O(\Delta x^6)$.

Below, I plot both the errors from part a) and part b) to compare them side by side. Clearly, both methods are $O(\Delta x^2)$.

plt.loglog(dx_list, error_list9, 'ro', label= '9-Point Laplacian')
plt.loglog(dx_list, error_list5, 'bo', label= '5-Point Laplacian')
plt.title('Comparison Between 5-Point Laplacian and 9-Point Laplacian')
plt.xlabel(r'$\Delta x$')
plt.ylabel('Max Error')
plt.legend()
plt.show()

(c) Calculate $\nabla^2_9 u(x,y) - f(x,y)$ at each grid point, for $f(x,y) = \sin(\pi x) \sin(\pi y)$. This is the discretization error of the 9-point Laplacian. Plot this as a 3-D surface plot. On a separate graph, plot the function $f_{xx} + f_{yy}$.

import numpy as np
import matplotlib.pyplot as plt

dx = 0.01
dy = dx

def sol(x, y):
    return 1/(2*np.pi**2)*np.sin(np.pi*x)*np.sin(np.pi*y)

def f(x,y):
    return -np.sin(np.pi*x)*np.sin(np.pi*y)

def dis_error(x,y):
    lap = (1/6*dx**2) * (4* sol(x - dx, y) + 4*sol(x+dx, y) + 4*sol(x, y-dx) + 4*sol(x, y+dx)
                        + sol(x - dx, y - dx) + sol(x - dx, y+ dx) + sol(x + dx, y - dx) +
                         sol(x + dx, y + dx) - 20* sol(x,y))
    return lap - f(x,y)

x_list = np.arange(x0, xN + dx/2, dx)
y_list = np.arange(y0, yN + dy/2, dy)

N = len(x_list)

dis_error_mat = np.zeros((N,N))

for i in range(0, len(x_list)):
    for j in range(0, len(y_list)):
        x = x_list[i]
        y = y_list[j]
        
        dis_error_mat[j,i] = dis_error(x,y)

X, Y = np.meshgrid(x_list, y_list)

ax = plt.axes(projection='3d')
# Plot the solution
ax.plot_surface(X, Y, dis_error_mat, alpha = 0.7, cmap='viridis')
ax.set_xlabel('x')
ax.set_ylabel('y')

plt.title('Discretization Error of 9-Point Laplacian')

plt.show()

Now I will plot $f_{xx} + f_{yy}$. This should look very similar to the plot of the discretization error.

import numpy as np
import matplotlib.pyplot as plt

def fxxfyy(x,y):
    return 2 * np.pi**2*np.sin(np.pi*x)*np.sin(np.pi*y)

fxxfyy_mat = np.zeros((N,N))

for i in range(0, len(x_list)):
    for j in range(0, len(y_list)):
        x = x_list[i]
        y = y_list[j]
        
        fxxfyy_mat[j,i] = fxxfyy(x,y)
        
X, Y = np.meshgrid(x_list, y_list)

ax = plt.axes(projection='3d')
# Plot the solution
ax.plot_surface(X, Y, fxxfyy_mat, alpha = 0.7, cmap='viridis')
ax.set_xlabel('x')
ax.set_ylabel('y')

plt.title(r'$f_{xx} + f_{yy}$')

plt.show()

It is known that the leading term in the discretization error is porportional to $f_{xx} + f_{yy}$, which is shown by the plots above. For this given problem, the 9-Point Laplacian and 5-Point Laplacian both had an error that was $O(\Delta x ^2)$, which may suggest that we have gained nothing from using a more complicated method.

However, for a different problem where $f(x,y)$ was chosen such that $f_{xx} + f_{yy} = 0$, the leading $O(\Delta x^2)$ error term for the 9-Point method would be eliminated, resulting in $O( \Delta x^6)$ error overall.