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
u_{tt}-c^{2}u_{xx}=0,\quad c>0,
with homogeneous Dirichlet boundary conditions
u(0,t)=u(1,t)=0,\quad\forall\,t\in\mathbb{R}
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 [1], 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 [1]. 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 pltFEniCS
def grad(function): return Dx(function, 0)
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, Vdef 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[0] = E T_arr[0] = T V_arr[0] = V u_t = 0.5*(u - u_n)/dt E, T, V = compute_energy_fenics(u_m, u_t, c) E_arr[1] = E T_arr[1] = T V_arr[1] = 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_arrTheory:
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, v0def 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, vdef 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_kdef 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_tdef 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, Vdef 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, Vdef 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[0]) ax1 = plt.subplot(gs[1]) ax2 = plt.subplot(gs[2]) 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() returndef 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() returnif name == ‘main’:
#plot_error_and_energy() compare_error_for_different_theta()
- Chabassier, J., & Imperiale, S. (2013). Introduction and study of fourth order theta schemes for linear wave equations. Journal of Computational and Applied Mathematics