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

Question

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}")
```

Answer 1

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.