 Theta scheme for 1D wave equation

Hello Community,

we are currently working on a complicated 1D PDE that involves a second-order time derivative. As a toy model, we are looking at the 1D wave equation to come up with a good discretization scheme for our actual problem. To be precise, we are looking at the 1D wave equation
with homogeneous Dirichlet boundary conditions
on the unit interval x\in [0,1] with initial values given by u(x,0)=u_0(x) and u_t(x,0)=v_0(x).

Multiplying the wave equation with the test function v and integrating over the domain of x, we can derive a weak formulation of the semi-discretized problem

\int_{0}^{1}u_{tt}v\,dx+c^{2}\int_{0}^{1}u_{x}v_{x}\,dx-\left.u_{x}v\right|_{0}^{1}=0

Note that the boundary term does not vanish since u_x is not zero at the boundaries. To arrive at a fully discretized problem, we approximate the second time derivative of u at time t_n=n \cdot \delta t by the second-order difference

u_{tt}=\frac{u^{n+1}-u^{n}+u^{n-1}}{\delta t^{2}}+\mathcal{O}(\delta t^{2}),

where u^n=u(n \cdot \delta t). For the second-order difference, the leading term of the truncation error scales quadratic with the time step \delta t. Substituting the approximation for u_{tt} into the weak formulation of the wave equation yields

\int_{0}^{1}\frac{u^{n+1}-u^{n}+u^{n-1}}{\delta t^{2}}v+c^{2}\int_{0}^{1}u_{x}^{n}v_{x}-\left.u_{x}^{n}v\right|_{0}^{1}=0,

which can be solved for u^{n+1} if u^n and u^{n-1} are known. Using the initial values u_0 and v_0, we can set u^0=u_0 and approximate u^1=u^0 + v_0\cdot \delta t which then allows us to compute u^2 which allows us to compute u^3 and so forth.

In the paper , the authors recap the three-point theta scheme for the wave equation while also introducing a new fourth-order scheme. The idea of the three-point theta scheme is that we do not use u^n in the second integral of the weak formulation, but instead use a weighted sum of the previous, the current, and the upcoming time steps given by

u_{\theta}^{n}=\theta u^{n+1}+(1-2\theta)u_{h}^{n}+\theta u_{h}^{n-1}.

For the three-point theta scheme, it can be shown that the leading term of the truncation error introduced by approximating u_{tt} by the second-order difference and u^n by u_{\theta}^{n} scales as

\delta t^{2}\left(\theta-\frac{1}{12}\right)+\mathcal{O}(\delta t^{4}),

as shown in equation (16) in . This suggests that we effectively end up with a fourth-order accurate scheme for \theta=1/12.

To verify this finding, we implemented the theta scheme for the 1D wave equation in FEniCS. See the MWE below. (Please skip the implementation of the analytic solution since it should not be relevant to my question)

Results

The top panel in the figure below compares the total energy E, the kinetic energy T, and the potential energy V for \theta=1/12. The bottom panel shows the relative error of the total energy compared to the analytic prediction which is simply a constant. The plots show a good agreement between numerical results and analytic prediction which gives me confidence in the solution obtained from the scheme.

To compute an error for the displacement u(x,t), we compare the numerical solution to the analytic solution which is given by a superposition of normal modes. Here, we are using the maximum absolute error for each time step given by

e(t_n)=\max_{x\in[0,1]}\left|u^{n}-u_{\text{theory}}^{n}\right|

The second figure shows that the error e(t) becomes smaller for larger \theta and does not show a minimum for \theta = 1/12 as I would naively expect. We would be thankful for any ideas on why e(t) shows the displayed behavior.

We also found that including the boundary term \left.u_{x}v\right|_{0}^{1}=0 does not make any difference. Is this to be expected?

Also, if somebody can hint us in the direction of other reliable discretization schemes for wave-like equations. This would be very much appreciated.

Thank you so much for taking the time to read this!

Third-party imports

import numpy as np
from fenics import *
import matplotlib.pyplot as plt

FEniCS

def compute_energy_fenics(u, u_t, c):

T = 0.5 * assemble(u_t**2 * dx)
V = 0.5 * assemble(c**2 * grad(u)**2 * dx)
E = T + V

return E, T, V

def solve_wave_equation(N, dt, T, u0, v0, theta, c = 1.0):

n = int(T/dt) + 1

# init mesh
mesh = UnitIntervalMesh(N - 1)

# init function spaces
FS = FunctionSpace(mesh, 'P', 1)

u = TrialFunction(FS)
v = TestFunction(FS)

# Matrix to store simulation displacement vector u for all time steps
u_mat = np.zeros((n, N))

# Init previous and current time step
u_mat[0, :] = u0
u_mat[1, :] = u0 + v0*dt

u_n = Function(FS)
u_m = Function(FS)

u_n.vector()[:] = u_mat[0, :]
u_m.vector()[:] = u_mat[1, :]

# Dirichlet boundary condition
bc = DirichletBC(FS, 0, lambda x, o: o)

# Weak formulation
F = (u - 2 * u_m + u_n) * v * dx + dt**2 * c**2 * (theta * grad(u) + (1 - 2*theta) * grad(u_m) + theta * grad(u_n)) * grad(v) * dx

#- (theta * grad(u) + (1 - 2*theta) * grad(u_m) + theta * grad(u_n)) * v * ds # Boundary term makes no difference

F_op, L = lhs(F), rhs(F)

# Arrays to store energy for all time steps
E_arr = np.zeros(n)
V_arr = np.zeros_like(E_arr)
T_arr = np.zeros_like(E_arr)

# Iterate through time
# Simulation starts at second time step
for i in range(2, n):

u = Function(FS)

solve(F_op == L, u, bcs = bc)

u_mat[i,:] = u.vector().get_local()

if i == 2:
# Compute Energy for first and second time step
u_t = 0.5*(-3*u_n + 4*u_m - u)/dt
E, T, V = compute_energy_fenics(u_n, u_t, c)
E_arr = E
T_arr = T
V_arr = V

u_t = 0.5*(u - u_n)/dt
E, T, V = compute_energy_fenics(u_m, u_t, c)
E_arr = E
T_arr = T
V_arr = V

# Compute Energy
u_t =  0.5*(3*u - 4*u_m + u_n)/dt

E, T, V = compute_energy_fenics(u, u_t, c)
E_arr[i] = E
T_arr[i] = T
V_arr[i] = V

# Update u
u_n.assign(u_m)
u_m.assign(u)

return u_mat, E_arr, T_arr, V_arr

Theory:

def init_single_mode(k, a, b, x_arr, c):

u0 = a * np.sin(k*pi*x_arr)
v0 = b * k * np.pi * c * np.sin(k*pi*x_arr)

return u0, v0

def init_super_position_of_modes(k_arr, a_arr, b_arr, x_arr, c):

u = np.zeros(len(x_arr))
v = np.zeros_like(u)

for i, k in enumerate(k_arr):

u_k, v_k = init_single_mode(k, a_arr[i], b_arr[i], x_arr, c)
u += u_k
v += v_k

return u, v

def compute_time_evolution_single_mode(k, a, b, t_arr, x_arr, c):

u_k = np.zeros((len(t_arr), len(x_arr)))
v_k = np.zeros((len(t_arr), len(x_arr)))

for i, t in enumerate(t_arr):

u_k[i, :] += (a * np.cos(k*np.pi*c*t) + b * np.sin(k*np.pi*c*t))*np.sin(k*pi*x_arr)
v_k[i, :] += (- a * k * np.pi * c * np.sin(k*np.pi*c*t) + b * k *np.pi * c *  np.cos(k*np.pi*c*t))*np.sin(k*pi*x_arr)

return u_k, v_k

def compute_time_evolution_superposition_of_modes(k_arr, a_arr, b_arr, x_arr, t_arr, c):

u_t = np.zeros((len(t_arr), len(x_arr)))
v_t = np.zeros_like(u_t)

for i, k in enumerate(k_arr):

u_k, v_k = compute_time_evolution_single_mode(k, a_arr[i], b_arr[i], t_arr, x_arr, c)
u_t += u_k
v_t += v_k

return u_t, v_t

def compute_energy_for_singe_mode(k, a, b, t_arr, c):

gamma = np.sqrt(a**2 + b**2)
phi = np.arctan(a/b)
E = 0.25 * c**2 * (k * gamma * pi)**2

T = E*np.cos(k * np.pi * t_arr + phi)**2
V = E*np.sin(k * np.pi * t_arr + phi)**2

return E, T, V

def compute_energy_for_superposition_of_modes(k_arr, a_arr, b_arr, t_arr, c):

E = np.zeros(len(t_arr))
T = np.zeros(len(t_arr))
V = np.zeros(len(t_arr))

for i, k in enumerate(k_arr):

E_k, T_k, V_k = compute_energy_for_singe_mode(k, a_arr[i], b_arr[i], t_arr, c)
E += E_k
T += T_k
V += V_k

return E, T, V

def compute_mode_coefficients(E, k_arr, c):

# Distribute energy evenly among the first ten modes
E_k = E/len(k_arr)

a_arr = np.zeros(len(k_arr))
b_arr = np.zeros(len(k_arr))

for i, k in enumerate(k_arr):

gamma = np.sqrt(4*E_k)/(c*np.pi*k)
a_arr[i] = gamma/np.sqrt(2)
b_arr[i] = a_arr[i]

return a_arr, b_arr

#===============================================================================

def plot_error_and_energy():

T = 5.0 # Simulation time
N  = 101 # Number of vertices
dt = 0.001 # time step
c = 1.0 # propagation speed
E = 100 # Total energy of the wave
theta = 1.0/12 # Theta

x_arr = np.linspace(0, 1, N)
t_arr = np.arange(0, T + 0.1*dt, dt)

# Initialize displacement vector u0 and velocity vector v0

# Set the first six normal modes active
k_arr = np.arange(1,6,1)
# Compute mode coeffients by distributing energy evenly among all modes
a_arr, b_arr = compute_mode_coefficients(E, k_arr, c)

# Initialize u0 and v0
u0, v0 = init_super_position_of_modes(k_arr, a_arr, b_arr, x_arr, c)

# Theory
u_theory, _ = compute_time_evolution_superposition_of_modes(k_arr, a_arr, b_arr, x_arr, t_arr, c)
E_theory, T_theory, V_theory = compute_energy_for_superposition_of_modes(k_arr, a_arr, b_arr, t_arr, c)

# FEniCS
u, E, T, V = solve_wave_equation(N, dt, T, u0, v0, theta, c)

# Compute error
max_err = np.max(np.abs((u - u_theory)), axis = 1)

# Plotting
lfz = 18

gs = plt.GridSpec(3, 1)
ax0 = plt.subplot(gs)
ax1 = plt.subplot(gs)
ax2 = plt.subplot(gs)

ax0.plot(t_arr, max_err)
ax0.set_ylabel(r'$e(t)$', fontsize = lfz)
ax1.plot(t_arr, E_theory, c = 'k', ls = '-', label = r'$E_{theory}$')
ax1.plot(t_arr, T_theory, c = 'b', ls = '-', label = r'$T_{theory}$')
ax1.plot(t_arr, V_theory, c = 'r', ls = '-', label = r'$V_{theory}$')
ax1.plot(t_arr, E, c = 'k', ls = '--', label = r'$E_{fenics}$')
ax1.plot(t_arr, T, c = 'b', ls = '--', label = r'$T_{fenics}$')
ax1.plot(t_arr, V, c = 'r', ls = '--', label = '$V_{fenics}$')
ax1.legend(fontsize = 12)
ax1.set_ylabel(r'Energy', fontsize = lfz)
ax2.plot(t_arr, (E - E_theory)/E_theory, c = 'k', ls = '-')
ax2.set_ylabel(r'Rel. error total energy', fontsize = lfz)
ax2.set_xlabel(r'$t$', fontsize = lfz)

plt.show()

return

def compare_error_for_different_theta():

T = 5.0 # Simulation time
N  = 101 # Number of vertices
dt = 0.001 # time step
c = 1.0 # propagation speed
E = 100 # Total energy of the wave

x_arr = np.linspace(0, 1, N)
t_arr = np.arange(0, T + 0.1*dt, dt)

# Initialize displacement vector u0 and velocity vector v0

# Set the first six normal modes active
k_arr = np.arange(1,6,1)
# Compute mode coeffients by distributing energy evenly among all modes
a_arr, b_arr = compute_mode_coefficients(E, k_arr, c)

# Initialize u0 and v0
u0, v0 = init_super_position_of_modes(k_arr, a_arr, b_arr, x_arr, c)

# theta_arr = np.arange(0, 2+0.01, 0.1)

theta_arr = np.array([0, 0.05, 1.0/12, 0.1, 0.2,  0.3, 0.4, 0.5, 0.75, 1.0])

col_arr = plt.cm.cool(theta_arr/theta_arr[-1])

# Theory
u_theory, _ = compute_time_evolution_superposition_of_modes(k_arr, a_arr, b_arr, x_arr, t_arr, c)

# FEniCS
for i, theta in enumerate(theta_arr):
u, _, _, _ = solve_wave_equation(N, dt, T, u0, v0, theta, c)

max_err = np.max(np.abs((u - u_theory)), axis = 1)
plt.plot(t_arr, max_err, c = col_arr[i], label = r'$\theta=%.2f$' % theta)

lfz = 18

plt.ylabel(r'$e(t)$', fontsize = lfz)
plt.xlabel(r'$t$', fontsize = lfz)
plt.legend()
plt.show()

return

if name == ‘main’:

#plot_error_and_energy()
compare_error_for_different_theta()
1. Chabassier, J., & Imperiale, S. (2013). Introduction and study of fourth order theta schemes for linear wave equations. Journal of Computational and Applied Mathematics

Hi, did you also consider Newmark methods ? These are classically used for dynamic analysis in structural mechanics.

Hey bleyerj, thank you for the suggestion, I will have a look. Interestingly, for \theta = 1/4 the stability of the scheme is improved for larger time steps as predicted by the theory. However, as mentioned in my post, the error is not minimal for \theta=1/12 as predicted. Since the theta-scheme is well known, I am prone to believe that it may have to do with FEniCS. Kind Regards, Lukas

The comparison of different \theta values for a fixed spatial mesh and time step, using an analytical solution to the PDE to compute error, is potentially misleading for multiple reasons:

• First, reducing only the time stepping error in the semi-discrete (viz., discretized in space only) problem does not necessarily bring you closer to the solution of the continuous (in both space and time) problem. You can end up better resolving the unsteady behavior of the spatial discretization error. This is the motivation behind the generalized-\alpha family of time integrators, which allow control over numerical dissipation in just the high-frequency modes of the problem, since that is where spatial discretization error is concentrated.
• Even for solving just an ODE system, the error at a fixed \Delta t depends on both the stability and the truncation error. The sense in which \theta = 1/12 is “most accurate” is that its truncation error converges to zero at a higher rate as \Delta t\to 0. That does not necessarily mean that, at a fixed \Delta t, its error is smaller than a lower-order method using the same value of \Delta t, since the constants in the error bound may be different.
1 Like