# Trouble Implementing 1d Wave Equation Finite Difference Solver

Im trying to solve the 1d Wave Equation on $$x \in \mathbb{R}, t > 0$$: $$u_{tt} = c^2u_{xx}, \hspace{5mm} u(x,0) = \cos(4 \pi x), \hspace{5mm} u_t(x,0) = 0$$ with $$c = 1$$ and a periodic boundary condition $$u(x,t) = u(x+1,t)$$ Im using central difference in time and space:

$$\frac{u_j^{n+1} - 2u_j^n + u_j^{n-1}}{k^2} = c^2 \frac{u_{j+1}^n - 2u_j^n + u_{j-1}^n}{h^2}$$ and the notation $$r = c\frac{k}{h}$$ which gives the iteration $$u_{j}^{n+1} = 2u_j^n - u_j^{n-1} + r^2(u_{j+1}^n - 2u_j^n + u_{j-1}^n)$$ And then for $$n = 0$$, to find $$u_j^{n-1}$$ I use the Neumann condition $$u_t(x,0) = \frac{u_j^{1} - u_j^{-1}}{2k} = 0$$ which gives $$u_j^{-1} = u_j^1$$. Inserting this into the original discretisation then gives for $$n = 0$$ $$u_j^1 = u_j^0 + \frac{1}{2}r^2(u_{j+1}^0 - 2u_j^0 + u_{j-1}^0)$$ Below is my implementation in Python:


class wave_eq:
def __init__(self,init_m,init_n):

self.M = init_m
self.N = init_n

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

self.solution = []

def forward_euler_solver(self,g):
M = self.M
N = self.N

h = self.x_grid[1] - self.x_grid[0]

k = self.t_grid[1] - self.t_grid[0]

r = k/h
u = np.zeros(shape = (M+2,N+1))

u[:,0] = g(self.x_grid) #apply initial condition u(x,0) = g(x)

for i in range(N+1):
new_u  = np.zeros(shape = (M+2,N+1))
for j in range(M+1):

if (i == 0): #using Neumann condition at first time step
new_u[j][i] = u[j][i] + ((r**2)/2)*(u[j+1][i] - 2*u[j][i] + u[j-1][I])
else:
new_u[j][i] = 2*u[j][i] - u[j][i-1] + (r**2)*(u[j+1][i] -2*u[j][i] + u[j-1][i])

new_u[-1][i] = new_u[0][i] #Apply periodic boundary condition

u = new_u

self.solution.append(u)

def init_c_g(x):
return np.cos(4*np.pi*x)



The relative $$\ell_2$$ error is very small, and comes out as $$7 \cdot 10^{-7}$$ but the solution does not converge and in fact grows slightly when adding more points in the $$x$$-space, as can be seen from the figure below where I plotted the relative $$\ell_2$$ error versus number of points in the $$x$$-space.

The analytical solution I used as reference to compute the plot is $$\frac{1}{2} \big (\cos(4 \pi (x-t) + \cos(4 \pi (x+t)) \big)$$I belive there is only a small error somewhere (perhaps when applying the boundary condition?) but I can't seem to figure it out. Would greatly appreciate any help.

• I believe your application of the perdiodic BC is wrong. You should not enforce the values of the solution itself. What you need to do is compute $\partial_{xx} u$ at all inner nodes (not the outer ones) first, using your standard centered-scheme. Then handle the outer nodes separately: for instance at point $x=0$, you conceptually introduce a "ghost point" at $x=-dx$, set its value to u[-1] (i.e. it will mimic the opposite boundary node), and apply your scheme: dudxx[0] = (u[1] - 2*u[0] + u[-1])/dx**2. Do a similar thing on the opposite side and that should work. – Laurent90 Mar 20 at 9:28
• btw, you should avoid using for loops in Python, and instead vectorize your code. There a lot of examples online on how to do that ;) – Laurent90 Mar 20 at 9:30
• @Laurent90 How could I solve for the ghost point at $x = -dx$? Since when solving for $u[0]$ I have not yet found $u[-1]$. Im new to finite differences so the method is still somewhat confusing to me. – Pame Mar 20 at 12:07
• just like when solving the heat equation I think: your vector $u$ has the values of the discrete solution at all nodes (i.e. u[0]and u[-1] are part of your solution and do not need to be solved for). The RHS at $c^2 \partial_{xx} u$ can be evaluated in a straightforward manner for the nodes 1 to N-2 as you've done. The "trick" to compute this RHS for the first point (and can be done for the alst one in a similar fashion) is to introduce an additional point befort that first point, i.e. ass if you had extended the domain by $dx$ to its left, and to use it as I have written before. – Laurent90 Mar 20 at 12:14