# High Running Time and Suboptimal Accuracy of 2D Wave Equation Solver with Finite Differences

Im trying to solve the following 2D wave equation: $$u_{tt} = u_{xx} + u_{yy}, \hspace{3mm} u(x,y,0) = \cos(4 \pi x) \sin(4 \pi y), \hspace{3mm} u_t(x,y,0) = 0$$ with the periodic boundary condition in both spatial directions $$u(x+1,y,t) = u(x,y,t), \hspace{3mm} u(x,y+1,t) = u(x,y,t)$$ which has analytical solution $$\cos(4 \pi x) \sin(4 \pi y) \cos(\sqrt{32}\pi t)$$

I used finite differences to solve the problem, with second order central differences in time, $$x$$ and $$y$$ from which I obtain an explicit iteration scheme. Details of the implementation can be viewed here under the "implementation" section: https://hplgit.github.io/fdm-book/doc/pub/book/sphinx/._book008.html

My issue is that Im not happy with the numerical accuracy of the method, and the running time of the program is also very high, which means I can't choose a large number of points in $$x$$- and $$y$$-directions to get higher accuracy as the code would take too long to process. Below is my code in Python and a convergence plot:


class wave_eq_2D:
def __init__(self,init_m,init_n,init_t):
self.M = init_m
self.N = init_n
self.T = init_t

self.x_grid = np.linspace(0, 1, self.M[0] + 2)
self.y_grid = np.linspace(0,1,self.N[0] + 2)
self.t_grid = np.linspace(0, 1, self.T + 2)

self.solution = []
self.rel_error = []

def wave_solver2(self,init_c):

M = self.M
N = self.N
T = self.T

h_x = self.x_grid[1] - self.x_grid[0]
h_y = self.y_grid[1] - self.y_grid[0]
k = self.t_grid[1] - self.t_grid[0]

r_x = k / h_x
r_y = k / h_y

u = np.zeros(shape=(M + 2, N + 2, T + 2)) #initialize numerical solution u

u[:, :, 0] = init_c(self.x_grid, self.y_grid) #apply initial condition g(x,y) = 4cos(4 * pi * x) * 4*sin(4* pi * x)

for i in range(T+1): #iterate in time

u_xx = u[2:, 1:-1, i] - 2 * u[1:-1, 1:-1, i] + u[:-2, 1:-1, i]
u_yy = u[1:-1,2:,i] -2*u[1:-1,1:-1,i] + u[1:-1,:-2,i]
if (i == 0):  # using Neumann condition at first time step

u[1:-1,1:-1, i + 1] = u[1:-1,1:-1, i] + 0.5*(r_x ** 2) * u_xx + 0.5*(r_y**2)*u_yy #apply time scheme for n = 0

#Apply periodic BCS using time scheme in both spatial dimensions: u(x+1,y,t) = u(x,y,t) = u(x,y+1,t)
u[0, 1:-1, i + 1] = u[0,1:-1, i] + 0.5*(r_x ** 2) * (u[1,1:-1, i] - 2 * u[0,1:-1, i] + u[-2,1:-1, i]) + (r_y**2)/2*(u[0,2:,i] -2*u[0,1:-1,i] + u[0,:-2,i])
u[-1,:,i+1] = u[0,:,i+1]
u[1:-1, 0, i + 1] = u[1:-1, 0, i] + 0.5 * (r_x ** 2)/2 * (u[2:, 0, i] - 2 * u[1:-1, 0, i] + u[:-2, 0, i]) + (r_y ** 2)/2 * (u[1:-1, 1, i] - 2 * u[1:-1, 0, i] + u[1:-1, -2, i])
u[:, -1, i + 1] = u[:, 0, i + 1]

else:

u[1:-1,1:-1, i + 1] = 2*u[1:-1,1:-1,i] -u[1:-1,1:-1,i-1] + (r_x ** 2)/2 * u_xx + (r_y**2)/2*u_yy #apply time scheme for general n

# Apply periodic BCS using time scheme in both spatial dimensions: u(x+1,y,t) = u(x,y,t) = u(x,y+1,t)
u[0, 1:-1, i + 1] = 2*u[0,1:-1,i] -u[0,1:-1,i-1] + 0.5*(r_x ** 2) * (u[1,1:-1, i] - 2 * u[0,1:-1, i] + u[-2,1:-1, i]) +(0.5*r_y**2)*(u[0,2:,i] -2*u[0,1:-1,i] + u[0,:-2,i])
u[-1, :, i + 1] = u[0, :, i + 1]  # on the right side
u[1:-1, 0, i + 1] = 2*u[1:-1,0,i] -u[1:-1,0,i-1] + (r_x ** 2) * (u[2:,0, i] - 2 * u[1:-1,0, i] + u[:-2,0, i]) +(r_y**2)*(u[1:-1,1,i] -2*u[1:-1,0,i] + u[1:-1,-2,i])
u[:,-1,i+1] = u[:,0,i+1]

self.solution = u[:,:,-1]



The convergence plot also looks weird as the convergence rate is not smooth. This makes me think there is an error somewhere. Convergence rate is approximately of second order which is as expected.

Would appreciate any feedback on how I could improve my wave equation solver implementation.

Edit: This is my full code and how I test convergence by running several iterations and doubling the number of points in $$x$$- and $$y$$ each time:


import numpy as np
import matplotlib.pyplot as plt
from matplotlib import cm

from plotting_tools import set_ax

class wave_eq_2D:
def __init__(self,init_m,init_n,init_t):
self.M = [init_m]
self.N = [init_n]
self.T = init_t

self.x_grid = [np.linspace(0, 1, self.M[0] + 2)]
self.y_grid = [np.linspace(0,1,self.N[0] + 2)]
self.t_grid = np.linspace(0, 1, self.T + 2)

self.solution = []
self.rel_error = []

def wave_solver2(self,init_c,iteration):

M = self.M[iteration]
N = self.N[iteration]
T = self.T

h_x = self.x_grid[iteration][1] - self.x_grid[iteration][0]
h_y = self.y_grid[iteration][1] - self.y_grid[iteration][0]
k = self.t_grid[1] - self.t_grid[0]

r_x = k / h_x
r_y = k / h_y

u = np.zeros(shape=(M + 2, N + 2, T + 2))

u[:, :, 0] = init_c(self.x_grid[iteration], self.y_grid[iteration])

for i in range(T+1):

u_xx = u[2:, 1:-1, i] - 2 * u[1:-1, 1:-1, i] + u[:-2, 1:-1, i]
u_yy = u[1:-1,2:,i] -2*u[1:-1,1:-1,i] + u[1:-1,:-2,i]
if (i == 0):  # using Neumann condition at first time step

u[1:-1,1:-1, i + 1] = u[1:-1,1:-1, i] + 0.5*(r_x ** 2) * u_xx + 0.5*(r_y**2)*u_yy

#Apply periodic BCS in both spatial dimensions: u(x+1,y,t) = u(x,y,t) = u(x,y+1,t)
u[0, 1:-1, i + 1] = u[0,1:-1, i] + 0.5*(r_x ** 2) * (u[1,1:-1, i] - 2 * u[0,1:-1, i] + u[-2,1:-1, i]) + (r_y**2)/2*(u[0,2:,i] -2*u[0,1:-1,i] + u[0,:-2,i])
u[-1,:,i+1] = u[0,:,i+1]
u[1:-1, 0, i + 1] = u[1:-1, 0, i] + 0.5 * (r_x ** 2)/2 * (u[2:, 0, i] - 2 * u[1:-1, 0, i] + u[:-2, 0, i]) + (r_y ** 2)/2 * (u[1:-1, 1, i] - 2 * u[1:-1, 0, i] + u[1:-1, -2, i])
u[:, -1, i + 1] = u[:, 0, i + 1]

else:

u[1:-1,1:-1, i + 1] = 2*u[1:-1,1:-1,i] -u[1:-1,1:-1,i-1] + (r_x ** 2)/2 * u_xx + (r_y**2)/2*u_yy

# Apply periodic BCS in both spatial dimensions: u(x+1,y,t) = u(x,y,t) = u(x,y+1,t)
u[0, 1:-1, i + 1] = 2*u[0,1:-1,i] -u[0,1:-1,i-1] + 0.5*(r_x ** 2) * (u[1,1:-1, i] - 2 * u[0,1:-1, i] + u[-2,1:-1, i]) +(0.5*r_y**2)*(u[0,2:,i] -2*u[0,1:-1,i] + u[0,:-2,i])
u[-1, :, i + 1] = u[0, :, i + 1]  # on the right side
u[1:-1, 0, i + 1] = 2*u[1:-1,0,i] -u[1:-1,0,i-1] + (r_x ** 2) * (u[2:,0, i] - 2 * u[1:-1,0, i] + u[:-2,0, i]) +(r_y**2)*(u[1:-1,1,i] -2*u[1:-1,0,i] + u[1:-1,-2,i])
u[:,-1,i+1] = u[:,0,i+1]

self.solution.append(u[:,:,-1])

print("Solution", iteration + 1, "found.")

def get_error(self, u, iteration):
# Defining the discrete l2 norm.
def disc_norm(v):
return np.sqrt(np.sum(v ** 2) / len(v))

# Defining a meshgrid from our x- and y-grids.
x_disc = self.x_grid[iteration]
y_disc = self.y_grid[iteration]
X_disc, Y_disc = np.meshgrid(x_disc, y_disc)

# Defining the numerical and analytical solution arrays
U_disc = self.solution[iteration]
u_disc = u(x_disc, y_disc,1)

# Finding the relative error using the l2 norm.
error_disc = disc_norm(u_disc - U_disc) / disc_norm(u_disc)
# Saving the error data in the relative error list.
self.rel_error.append(error_disc)

def check_convergence(self, iterations, u, init_c):

# Finding the l2 error of the initial system.
self.wave_solver2(init_c, 0)
self.get_error(u, 0)

# Calculating the l2 errors of a few more iterations of the system.
for i in range(1, iterations):
# Doubling the size of the x-grid and y-grid.
self.M.append(2 * self.M[i - 1])
M = self.M[i]
self.x_grid.append(np.linspace(0, 1, M + 2))
self.N.append(2*self.N[i-1])
N = self.N[i]
self.y_grid.append(np.linspace(0,1,N + 2))

# Solving for the new system.
self.wave_solver2(init_c, i)
self.get_error(u, i)

def plot_convergence(self, line_order=2, space='x'):
disc_error = self.rel_error
print('M-values:', self.M)

if space == 'x':
gridsize = self.M
x_label = r'$$M$$'
line_label = r'$$O(h^{%s})$$'
elif space == 'y':
gridsize = self.N
x_label = r'$$N$$'
line_label = r'$$O(k^{%s})$$'
else:
return

x_array = np.linspace(gridsize[0], gridsize[-1], 1000)

def line(x, order, init_error):
return init_error * x ** (-order) / (x[0]) ** (-order)

print('M-values:', gridsize)
print('disc error:', disc_error)

fig = plt.figure()
set_ax(ax,
title='Convergence testing for 2D Wave Equation',
xscale='log', yscale='log',
x_label=x_label, y_label='Relative error')

ax.plot(gridsize, disc_error, marker='o', label=r'$$l_2$$ norm error', lw=3, c='#1F77B4')
ax.plot(x_array, line(x_array, line_order, disc_error[0]), label=line_label % line_order,
lw=3, ls='--', c='#1F77B4')

ax.legend(fontsize=15)
plt.show()
plt.close()

#Boundary condition g(x,y)
def g(x,y):
return np.cos(4* np.pi * x) * np.sin(4 * np.pi * y)

#Analytical solution
def u(x,y,t):
return np.cos(4*np.pi*x)*np.sin(4*np.pi*y)*np.cos(np.sqrt(32)*np.pi*t)

System = wave_eq_2D(70,70,1550)

System.check_convergence(4,u,g)
System.plot_convergence()

$$$$

• I do not understand why you have computations in the boundary conditions. Shouldn't you just copy the already existing data to have a periodic border for the next Laplace computation, u[0,1:-1,i+1] = u[-2,1:-1,i+1]; u[-1,1:-1,i+1]=u[1,1:-1,i+1]; u[:,0,i+1]=u[:,-2,i+1]; u[:,-1,i+1]=u[:,1,i+1]? – Lutz Lehmann Mar 22 at 20:38
• You seem to mix 2 ideas for the periodic continuation. You use the identification of index 0 with index -1. For a grid with $M\times N$ rectangles, this would require that the subdivisions be constructed using length $M+1$ and $N+1$ in linspace. But you use $M+2$ and $N+2$, which would make sense if, as I assumed in my previous comment, you identify the first two indices [0,1] with the last two [-2,-1], so that $x_N=1$, $x_{N+1}=1+h_x$ etc. Then you also do not need the extra computations for row and column 0. – Lutz Lehmann Mar 23 at 8:26
• @LutzLehmann My idea was to have $x_0 = 0, x_{M+1} = 1$, but wouldn't then $x_0$ and $x_{M+1}$ be equal? With $x_0$ = $u[0]$ and $x_{M+1}$ being $u[-1]$ i.e point number $M+2$. – Pame Mar 23 at 10:26
• There is nothing wrong in the code in that regard, I just wonder about the semantics of $M$ and $N$. So if $M=N=10$, you get a grid of $11\times 11$ squares. Is that your intention, or is the $M+2$ a left-over from a Matlab implementation where the array indices start at $1$? – Lutz Lehmann Mar 23 at 10:50

Regarding performance, Python is definitely the bottleneck. I have experienced the same issue with a 2D Euler code I had developed, even with vectorised operations everywhere possible. It was actually even worse, as I was using solve_ivp time schemes which reallocated memory at every step... You can try and profile your code to see where the bottlenecks are. A profiler is included in the Spyder IDE (shortcut F10). Can you share your code (with the convergence analysis) ?
Also, I see that there are sometimes (r_x ** 2)/2, sometimes only (r_x ** 2)`in your code. Have you double-checked that this is correct? I remember that there was an issue with that in your previous question for the 1D case.
• I edited my post now to supply the full code where I calculate error and convergence plots. Regarding $r_x^2$ vs $r_x^2$/2 it is because the formula changes for $n = 0$. It is explained in the link I posted. Since the code runs slowly when choosing a high number of points in $x$ and $y$ means very high running time which is why I have to choose the low for the convergence plot. Im also wondering if you are able to achieve better results with your own revised version of this code? (without using temporal integration). – Pame Mar 23 at 9:43