# Error too large in leapfrog method for solving the wave equation of a vibrating string

I have been trying to figure out what I did wrong for the last two days. I do not know if I actually did something wrong or if the error is supposed to be this large in usual leapfrog problems. I tried checking my numerical solution with an exact one. Specifically, where the initial condition of the string is $$y(x,0) = x(1-x)$$ and $$\frac{\partial y}{\partial t} \vert_{t=0} = 0$$ which has the exact solution $$y(x,t) = \sum_{k=1\\ \textrm{k odd}}^{\infty} \frac{8}{k^3 \pi^3} \sin{}k \pi x \cos{v k \pi t}$$ The pde was $$v^2 y_{xx} = y_{tt}$$

The first initial condition gives me $$y_{i,0}$$ for all $$i$$ and the second initial condition gives me $$y_{i,1} = y_{i,0} + \frac{C^2}{2}(y_{i+1, 0} + y_{i-1, 0} - 2y_{i,0})$$ where C is the Courant number $$v\Delta t / \Delta x$$ Further time steps are found using $$y_{i,j+1} = 2y_{i,j} - y_{i,j-1} + C^2(y_{i+1,j} + y_{i-1,j} - 2y_{i,j})$$

I made sure the Courant condition was satisfied. The amplitude of the wave is 0.25 units, while there is an error of 0.1 units after just 316 time steps. Which is about 0.316s. I don't think leapfrog method has such weak predictability (does it?). If I increase time stepsize, then the number of time steps until error reaches 0.1 just decreases by the same ratio. I tried matching this with another exact solution today and the errors were even larger. I am convinced there is something wrong with my approach but I can't seem to figure out what. I am pretty confident my recurrence relations are correct. If the error really is supposed to be this big, then is there any analytic formula for the error for this sort of problems which I could use to choose my parameters properly?

Here is my code in python (edited with a faster algorithm, since I only need to know the values at the last two time steps to calculate the values at the next time step) :

import numpy as np
from math import *

# Exact Soln
def y_exact(x,t):
if t==0:
return x*(1-x)
else:
sum = 0
for i in range(1,n):
if i%2 == 1:
sum += 8*sin(i*pi*x)*cos(v*i*pi*t)/((i**3)*(pi**3))
return sum

n = 200
v = 1
tolerance = 0.1
L = 1
T = 1
x = np.linspace(0,L,10)
t = np.linspace(0,T,1000)
dx = x[1] - x[0]
dt = t[1] - t[0]
c = dx/dt
C = v/c
C2 = C*C
nx = len(x)
y = np.zeros((nx,3))
max_error = 0

# Initial condition 1
for i in range(nx):
y[i,0] = y_exact(x[i],0)

# Initial condition 2
for i in range(nx):
if i == 0 or i == (nx - 1) : # Boundary Condition
y[i,1] = 0
else:
y[i,1] = y[i,0] + 0.5*C2*(y[i+1,0] + y[i-1,0] - 2*y[i,0])

for j in range(2,len(t)):
if max_error >= tolerance :
print(f"Error exceeded {tolerance} at time step number {j-1}.")
break
for i in range(nx):
if i==0  or i == (nx -1):
y[i,2] = 0
else:
y[i,2] = 2*y[i,1] - y[i,0] + C2*(y[i+1,1] + y[i-1,1] - 2*y[i,1])
error = abs(y[i,2] - y_exact(x[i],dt*j))
if error >= max_error :
max_error = error
for k in range(nx):
y[i,0] = y[i,1]
y[i,1] = y[i,2]

print(f"Maximum value of error was {max_error}")
$$$$

• Did you check that your x sequence ends at 1 (it does not)? You can force this with x = np.arange(0,1+0.5*dx,dx) so that the last node is at least close to 1, but better would be x=np.linspace(0,1,N+1) where N=round(1/h) and then h=1/N or h=x[1]-x[0]. There is also an option to return the stepsize in linspace. – Lutz Lehmann Jan 12 at 15:37
• @LutzLehmann I don't think that matters. If it doesn't end at 1, it just means the length of the string is (1-dx) units instead of 1 unit. I did what you said btw, replacing arange with linspace and testing the result. The error is the same. – Brain Stroke Patient Jan 12 at 16:16
• Do you only change dt? Then the exact solution to compare against is the solution of the ODE system corresponding to the $x$-discretization, which is $O(dx^2)$ distant from the fully exact solution . You only get convergence to the fully exact solution if you keep C about constant. – Lutz Lehmann Jan 12 at 16:44
• No I tried changing just dx, just dt and both keeping C the same or increasing/decreasing it. If I decreased the time stepsize from say 0.001 to 0.0001 for example, then the number of time steps till error > 0.1 changes from 316 to 3162 but the result is the same. Error exceeds 0.1 after 0.316 s. Changing dx to keep C constant or not changing dx only changes the error from 0.104... to 0.100... , but the error exceeds 0.1 after 0.316 s nonetheless. – Brain Stroke Patient Jan 12 at 17:51
• I tried bumping up the number of terms summed from 200 to 1000, but the error was still greater than 0.1. The only change to the error was in the 6'th or 7th significant figure which means it's not an issue of my exact solution not being exact enough, at least. – Brain Stroke Patient Jan 12 at 17:57

I used another variant for the exact solution of the wave equation (which is the same as the series sum), and shortened the computation by transforming loops into numpy array operations. With this changed code I do not observe the numerical errors you got

v = 1
def f(x):
'''2-periodic odd-symmetric continuation of y(x,0)'''
x = (x%2+1)%2-1;
return x*(1-abs(x))
def y_exact(x,t): return 0.5*(f(x+v*t)+f(x-v*t))

tolerance = 0.1
L = 1
T = 10
x = np.linspace(0,L,10)
t = np.linspace(0,T,1000)
dx = x[1] - x[0]
dt = t[1] - t[0]
c = dx/dt
C = v/c
C2 = C*C
nx = len(x)
y = np.zeros((nx,3))
max_error = 0

# Initial condition 1
y[:,0] = y_exact(x,0)

# Initial condition 2
y[[0,-1],1] = 0
y[1:-1,1] = y[1:-1,0] + 0.5*C2*(y[2:,0] + y[:-2,0] - 2*y[1:-1,0])

for j in range(2,len(t)):
if max_error >= tolerance :
print(f"Error exceeded {tolerance} at time step number {j-1}.")
break
y[[1,-1],2] = 0
y[1:-1,2] = 2*y[1:-1,1] - y[1:-1,0] + C2*(y[2:,1] + y[:-2,1] - 2*y[1:-1,1])
error = max(abs( y[:,2] - y_exact(x,dt*j) ))
max_error = max(max_error,error)

y[:,0] = y[:,1]
y[:,1] = y[:,2]

print(f"Maximum value of error was {max_error}")


This runs to the end, with a maximal reported error of 0.005687, well below the tolerance. One would have to check more carefully where the replacements erased a coding error.

Found it. In the last step, the wrap-around, you have

    for k in range(nx):
y[i,0] = y[i,1]
y[i,1] = y[i,2]


the iteration variable is k, the array index used is i`.

• Wow that was very sloppy of me. Thanks! – Brain Stroke Patient Jan 12 at 19:53