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}")
```
x
sequence ends at1
(it does not)? You can force this withx = np.arange(0,1+0.5*dx,dx)
so that the last node is at least close to 1, but better would bex=np.linspace(0,1,N+1)
whereN=round(1/h)
and thenh=1/N
orh=x[1]-x[0]
. There is also an option to return the stepsize inlinspace
. $\endgroup$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 keepC
about constant. $\endgroup$