## 5 The Non-Linear Seismic Wave Equation  
In the second concrete application with our method, we seek solutions of the non-linear seismic wave equation. Its general form, which describes the propagation of waves in elastic media, is [22]  
\[-\frac{1}{\left[c(x,y,z)\right]^{2}}u_{tt}+u_{xx}+u_{yy}+u_{zz}=0, \tag{22}\]  
where the function \(c(x,y,z)\) is the wave velocity reflecting spatially variable elastic properties of the medium. It should be noted that the non-linearity of the seismic wave Equation (22) enters through the square of the spatially dependent velocity \(c(x,y,z)\). Equation (22) has two types of solutions that correspond to compressional and shear waves, which share many common characteristics with the solutions of the \(1\)D wave Equation (15). Moreover, an expression equivalent to (22) can be obtained by starting from the Hamiltonian system (1) where [14]  
\[M= \begin{bmatrix}0&-1/c(x)&0\\ 1/c(x)&0&0\\ 0&0&0\end{bmatrix} \qquad\text{and}\qquad K= \begin{bmatrix}0&0&1\\ 0&0&0\\ -1&0&0\end{bmatrix} , \tag{23}\]  
for \(z=[u\;\;v\;\;w]^{T}\) and \(S(z)=(v^{2}-w^{2})/2\). The pre-symplectic forms \(\omega\) and \(\kappa\) given by (2) can also be defined so that the Hamiltonian system satisfies the multisymplectic conservation law (3). In addition, Equation (22) can also be derived from the Lagrangian  
\[L\big{(}u,u_{t},u_{x},u_{y},u_{z}\big{)}=\frac{1}{2c^{2}}u_{t}^{2}-\frac{1}{2} \Big{(}u_{x}^{2}+u_{y}^{2}+u_{z}^{2}\Big{)}, \tag{24}\]  
Figure 3: Evolution of the discrete energy error for the linear wave equation using the standard multisymplectic scheme of [11] (blue) versus the proposed one (red). Top: temporal evolution of the discrete energy error (20). Bottom: spatial evolution of the discrete energy error (21). inserted into the Euler-Lagrange field Equations (7). For simplicity, we consider a plane wave propagating in \(x\)-direction; i.e., we seek solutions of (22) for compressional waves travelling in \(x\)-direction. The resulting equation is written as  
\[-\frac{1}{\left[c(x)\right]^{2}}u_{tt}+u_{xx}=0. \tag{25}\]  
By discretising the corresponding Lagrangian of (24) and using a simple quadrature rule for the discrete version of the wave velocity, we obtain the integrator  
\[-\frac{1}{\left[c\left(\frac{x_{i}+x_{i+1}}{2}\right)\right]^{2}}\frac{u_{i}^{ j+1}-2u_{i}^{j}+u_{i}^{j-1}}{\left[\phi(\Delta t)\right]^{2}}+\frac{u_{i+1}^{j}-2 u_{i}^{j}+u_{i-1}^{j}}{\left[\psi(\Delta x)\right]^{2}}=0. \tag{26}\]  
Furthermore, in order to compare the latter integrator with previously presented schemes, we follow the setup of [27]. The wave velocity is assumed to be given by  
\[[c(x)]^{2}=1+0.1e^{-x^{2}}. \tag{27}\]  
The solution is sought in the domain \(-10<x<10\), with the initial conditions [27]  
\[u(x,0) = sechx,\] \[u_{t}(x,0) = 0,\]  
and with the periodic boundary condition  
\[u(-10,t)=u(10,t). \tag{28}\]  
Figure 4 shows the wave-forms of the seismic wave obtained via the discrete Equation (26) with \(\Delta x=0.1\) and \(\Delta t=0.01\). In addition, Figure 5 illustrates the wave-forms at four distinct time instances: \(t=0\), \(t=4\), \(t=10\) and \(t=15\). These have been selected following [27] to visualise the split of energy of the initial wave to two smaller but equal waves travelling in different (opposite) directions, which, due to periodic boundary conditions, reform the initial wave. Figure 4: Calculated waveforms of the non-linear seismic wave Equation (22) for \(\Delta x=0.1\) and \(\Delta t=0.01\). In a similar way to that followed in the previous Section, we calculated the evolution of the relative energy error in the time interval \(t\in[0,400]\), as the difference between the current and the initial energy. The results are plotted in Figure 6, where the blue dashed line corresponds to the error obtained by the standard multisymplectic scheme of [11], which is identical to the one presented in [27], and the red line represents the error evolution in the new scheme.