Consider the intial-boundary value problem \begin{equation} \begin{aligned} u_t &= \kappa (u_{xx} + u_{yy}) \quad \text{for all } 0 < x, y < 1 \text{ and } t > 0 \\ u(t,0,y) &= u(t,1,y) = 0 \quad \text{for all } 0 \le y \le 1 \text{ and } t > 0 \\ u(t,x,0) &= u(t,x,1) = 0 \quad \text{for all } 0 \le x \le 1 \text{ and } t > 0 \\ u(0,x,y) &= \begin{cases} e^{-10((x-0.5)^2 + (y-0.5)^2)} & 0 < x, y < 1 \\ 0 & x = 0, y = 0, x = 1 \text{ or } y = 1 \end{cases} \end{aligned} \end{equation}

where $\kappa = 0.1$.

In this problem, you wil solve the 2-D PDE using the method of lines. For all parts, use a second order central difference scheme for both $u_{xx}$ and $u_{yy}$ with $\Delta x = \Delta y = 0.1$

At a given time step $t$, I construct a vector of the solution at every grid point in the following way.

$$\vec u(t) = \begin{bmatrix} u_{0,0}(t) \\ u_{0,1}(t) \\ ...\\ u_{N-1,0}(t) \\ u_{0,1}(t) \\ u_{1,1}(t) \\ ...\\ u_{N-1,1}(t)\\ ...\\ u_{N-1,N-1}(t) \\ ... \end{bmatrix}$$

where $u(t_k, x_m, y_n)$ = $u_{m,n}(t_k)$

By using the 2nd order central difference scheme for both $u_{xx}$ and $u_{yy}$, each internal grid point produces a equation:

$$\partial_t u_{m,n} =\frac{\kappa}{\Delta x^2} \left( u_{m+1,n} + u_{m-1,n} + u_{m,n+1} + u_{m,n-1} - 4u_{m,n} \right)$$

We obtain a system equations:

$$\vec u_t = A \vec u(t) = f(t, \vec u(t))$$

In this form, time integration is done with any typical numerical IVP solvers.

(a) Use the forward Euler to solve the PDE with $\Delta t = 1/3$. Make a surface plot of the solution at $t = 0,\space t = 1/3, \space t = 2/3, \text{ and } t = 1$.

Forward Euler time integration is given by: $$\vec{u}(t_{k+1}) = \vec{u}(t_k) + \Delta t f(t_k,\vec{u}(t_k))$$

import numpy as np
import matplotlib.pyplot as plt
import scipy

x0 = 0
xf = 1
y0 = 0
yf = 1

dx = 0.1
dy = 0.1

x= np.arange(x0, xf + dx/2, dx)
y= np.arange(x0, yf + dy/2, dy)

N = len(x)

t0=0
tf = 1
dt = 1/3
t= np.arange(t0, tf + dt/2, dt)

Nt = len(t)

def point2ind(m, n):# Each row in the matrix has all points in the grid. So this function allows you to input the grid coordinate and output the row index for that coordinate
    return n * N + m


#fill in A
A = np.zeros((N**2,N**2))
U = np.zeros((N**2, len(t))) #every column is the flattened lattice
for n in range(N): #u_(m,n)
    for m in range(N): #looping over x y grid
        k = point2ind(m, n) #get row index for some point in the grid
        if n == 0: #if y coordinate is 0, bottem boundary
            A[k, k] = 1
        elif m == 0 or m == N-1 or n == N-1: #all other boundaries
            A[k, k] = 1
        else: #all other interior points
            A[k, k] = -4 / dx ** 2
            A[k, k + 1] = 1 / dx ** 2
            A[k, k - 1] = 1 / dx ** 2
            A[k, k + N] = 1 / dx ** 2
            A[k, k - N] = 1 / dx ** 2

A = 0.1*A

def f(t, u):
    return A @ u

#set up initial condition for U, eveyrthing but the boundaries
for m in range(1,len(x)-1):
    for n in range(1, len(y)-1):
        k = point2ind(m,n)
        U[k,0] = np.exp(-10* ((x[m] - 0.5)**2 + (y[n] - 0.5)**2))

#loop over forward euler
for k in range(Nt - 1):
    U[:, (k + 1):(k + 2)] = U[:, k:(k + 1)] + dt * f(t[k], U[:, k:(k + 1)])


U1 = U[:,0].reshape((N,N))
U2 = U[:,1].reshape((N,N))
U3 = U[:,2].reshape((N,N))
U4 = U[:,3].reshape((N,N))

from mpl_toolkits.mplot3d import Axes3D 

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

fig = plt.figure(figsize=(12, 8))

# Plot for U1
ax1 = fig.add_subplot(2, 2, 1, projection='3d')
ax1.plot_surface(X, Y, U1, cmap='viridis')
ax1.set_title(r'$t = 0$')
#ax1.set_zlim([-10, 10])

# Plot for U2
ax2 = fig.add_subplot(2, 2, 2, projection='3d')
ax2.plot_surface(X, Y, U2, cmap='viridis')
ax2.set_title(r'$t = 1/3$')
#ax2.set_zlim([-10, 10])

# Plot for U3
ax3 = fig.add_subplot(2, 2, 3, projection='3d')
ax3.plot_surface(X, Y, U3, cmap='viridis')
ax3.set_title(r'$t = 2/3$')
#ax3.set_zlim([-10, 10])

# Plot for U4
ax4 = fig.add_subplot(2, 2, 4, projection='3d')
ax4.plot_surface(X, Y, U4, cmap='viridis')
ax4.set_title(r'$t = 1$')
#ax4.set_zlim([-10, 10])

fig.suptitle(r'Forward Euler with $\Delta t = 1/3$', fontsize=16)

# Adjust layout
plt.tight_layout()

# Show the plot
plt.show()
    

(b) Use the forward Euler method to solve the PDE with $\Delta t = 0.01$ Make a surface plot of the solution at times $t=0, \space t=0.33, \space t = 0.67, \text{ and } t= 1$.

t0=0
tf = 1
dt = 0.01
t= np.arange(t0, tf + dt/2, dt)

Nt = len(t)

U = np.zeros((N**2, len(t))) #every column is the flattened lattice

#set up initial condition for U, eveyrthing but the boundaries
for m in range(1,len(x)-1):
    for n in range(1, len(y)-1):
        k = point2ind(m,n)
        U[k,0] = np.exp(-10* ((x[m] - 0.5)**2 + (y[n] - 0.5)**2))
        
#loop over forward euler
for k in range(Nt - 1):
    U[:, (k + 1):(k + 2)] = U[:, k:(k + 1)] + dt * f(t[k], U[:, k:(k + 1)])



U1 = U[:,0].reshape((N,N))
time = np.where(t == 0.33)[0]
U2 = U[:,time].reshape((N,N))
time = np.where(t == 0.67)[0]
U3 = U[:,time].reshape((N,N))
time = np.where(t == 1)[0]
U4 = U[:,time].reshape((N,N))

from mpl_toolkits.mplot3d import Axes3D 

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

fig = plt.figure(figsize=(12, 8))

# Plot for U1
ax1 = fig.add_subplot(2, 2, 1, projection='3d')
ax1.plot_surface(X, Y, U1, cmap='viridis')
ax1.set_title(r'$t = 0$')
ax1.set_zlim([0, 1])

# Plot for U2
ax2 = fig.add_subplot(2, 2, 2, projection='3d')
ax2.plot_surface(X, Y, U2, cmap='viridis')
ax2.set_title(r'$t = 0.33$')
ax2.set_zlim([0, 1])

# Plot for U3
ax3 = fig.add_subplot(2, 2, 3, projection='3d')
ax3.plot_surface(X, Y, U3, cmap='viridis')
ax3.set_title(r'$t = 0.67$')
ax3.set_zlim([0, 1])

# Plot for U4
ax4 = fig.add_subplot(2, 2, 4, projection='3d')
ax4.plot_surface(X, Y, U4, cmap='viridis')
ax4.set_title(r'$t = 1$')
ax4.set_zlim([0, 1])

fig.suptitle(r'Forward Euler with $\Delta t = 0.01$', fontsize=16)

# Adjust layout
plt.tight_layout()

# Show the plot
plt.show()

(c) Use the forward trapezoidal method to solve the PDE with $\Delta t = 1/3$. Make a surface plot of the solution at $t = 0,\space t = 1/3, \space t = 2/3, \text{ and } t = 1$.

The trapezoid method is given by: $$\vec u_k = \vec u_k + \frac{\Delta t}{2} \left(f(t_k, \vec u_k) + f(t_{k+1}, \vec u _{k+1}) \right)$$ $$\vec u_k = \vec u_k + \frac{\Delta t}{2} \left( A\vec u_k + A\vec u _{k+1} \right)$$

Rearranging, we have

$$(I - \frac{\Delta t}{2}A) \vec u_{k+1} = (I + \frac{\Delta t}{2}A) \vec u_k$$

t0=0
tf = 1
dt = 1/3
t= np.arange(t0, tf + dt/2, dt)

Nt = len(t)

U = np.zeros((N**2, len(t))) #every column is the flattened lattice

I = np.eye(*A.shape)

#set up initial condition for U, eveyrthing but the boundaries
for m in range(1,len(x)-1):
    for n in range(1, len(y)-1):
        k = point2ind(m,n)
        U[k,0] = np.exp(-10* ((x[m] - 0.5)**2 + (y[n] - 0.5)**2))
        
#loop over trapazoid
for k in range(Nt - 1):
    U[:, (k + 1):(k + 2)] = np.linalg.solve(I - dt*A/2, (I + dt*A/2) @ U[:, k:(k + 1)])



U1 = U[:,0].reshape((N,N))
#time = np.where(t == 0.33)[0]
U2 = U[:,1].reshape((N,N))
#time = np.where(t == 0.67)[0]
U3 = U[:,2].reshape((N,N))
#time = np.where(t == 1)[0]
U4 = U[:,3].reshape((N,N))

from mpl_toolkits.mplot3d import Axes3D 

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

fig = plt.figure(figsize=(12, 8))

# Plot for U1
ax1 = fig.add_subplot(2, 2, 1, projection='3d')
ax1.plot_surface(X, Y, U1, cmap='viridis')
ax1.set_title(r'$t = 0$')
ax1.set_zlim([0, 1])

# Plot for U2
ax2 = fig.add_subplot(2, 2, 2, projection='3d')
ax2.plot_surface(X, Y, U2, cmap='viridis')
ax2.set_title(r'$t = 1/3$')
ax2.set_zlim([0, 1])

# Plot for U3
ax3 = fig.add_subplot(2, 2, 3, projection='3d')
ax3.plot_surface(X, Y, U3, cmap='viridis')
ax3.set_title(r'$t = 2/3$')
ax3.set_zlim([0, 1])

# Plot for U4
ax4 = fig.add_subplot(2, 2, 4, projection='3d')
ax4.plot_surface(X, Y, U4, cmap='viridis')
ax4.set_title(r'$t = 1$')
ax4.set_zlim([0, 1])

fig.suptitle(r'Trapezoid with $\Delta t = 1/3$', fontsize=16)

# Adjust layout
plt.tight_layout()

# Show the plot
plt.show()

(d) Use the forward trapezoidal method to solve the PDE with $\Delta t = 0.01$. Make a surface plot of the solution at $t = 0,\space t = 0.33, \space t = 0.67, \text{ and } t = 1$.

t0=0
tf = 1
dt = 0.01
t= np.arange(t0, tf + dt/2, dt)

Nt = len(t)

U = np.zeros((N**2, len(t))) #every column is the flattened lattice

I = np.eye(*A.shape)

#set up initial condition for U, eveyrthing but the boundaries
for m in range(1,len(x)-1):
    for n in range(1, len(y)-1):
        k = point2ind(m,n)
        U[k,0] = np.exp(-10* ((x[m] - 0.5)**2 + (y[n] - 0.5)**2))
        
#loop over trapazoid
for k in range(Nt - 1):
    U[:, (k + 1):(k + 2)] = np.linalg.solve(I - dt*A/2, (I + dt*A/2) @ U[:, k:(k + 1)])



U1 = U[:,0].reshape((N,N))
time = np.where(t == 0.33)[0]
U2 = U[:,time].reshape((N,N))
time = np.where(t == 0.67)[0]
U3 = U[:,time].reshape((N,N))
time = np.where(t == 1)[0]
U4 = U[:,time].reshape((N,N))

from mpl_toolkits.mplot3d import Axes3D 

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

fig = plt.figure(figsize=(12, 8))

# Plot for U1
ax1 = fig.add_subplot(2, 2, 1, projection='3d')
ax1.plot_surface(X, Y, U1, cmap='viridis')
ax1.set_title(r'$t = 0$')
ax1.set_zlim([0, 1])

# Plot for U2
ax2 = fig.add_subplot(2, 2, 2, projection='3d')
ax2.plot_surface(X, Y, U2, cmap='viridis')
ax2.set_title(r'$t = 1/3$')
ax2.set_zlim([0, 1])

# Plot for U3
ax3 = fig.add_subplot(2, 2, 3, projection='3d')
ax3.plot_surface(X, Y, U3, cmap='viridis')
ax3.set_title(r'$t = 2/3$')
ax3.set_zlim([0, 1])

# Plot for U4
ax4 = fig.add_subplot(2, 2, 4, projection='3d')
ax4.plot_surface(X, Y, U4, cmap='viridis')
ax4.set_title(r'$t = 1$')
ax4.set_zlim([0, 1])

fig.suptitle(r'Trapezoid with $\Delta t = 0.01$', fontsize=16)

# Adjust layout
plt.tight_layout()

# Show the plot
plt.show()

Which of the solutions exhibit stability issues? Looking at the plots above, Forward Euler had instability issues when $\Delta t = 1/3$, which is what we expect. With a smaller $\Delta t$, Forward Euler can become stable. Being an implicit method, Trapezoidal method is always stable.

Which method would be the best to get a rough picture of the solution at $\\\textbf{ t = 0, t = 1/3, t = 2/3, t = 1?}$

Trapezoid method with $\Delta t = 1/3$ did fine job at getting a rough picture of the PDE at different times. This is the clear choice, as we only need to iterate 4 times, as opposed to the methods with $\Delta t = 0.01$, which would take 30 times more iterations.