## 7 Numerical Examples  
In this section, we present numerical examples justifying the theoretical results of the previous two sections. For convergence tests the domain decomposition algorithm (see Algorithm 2), implemented in the software package WavES [48], was used. We note that because of using explicit FE and FD schemes in \(\Omega_{FEM}\) and \(\Omega_{FDM}\), correspondingly, we need to choose a time step \(\tau\) according to the CFL stability condition (9) derived in [44] so that the whole hybrid scheme remains stable. Numerical tests are performed in time interval \((0,T)=(0,0.25)\) and in the spatial dimensionless computational domain  
\[\Omega=\{(x,y):x\in[0,1],y\in[0,1]\}, \tag{74}\]  
which is split into the finite element domain  
\[\Omega_{FEM}=\{(x,y):x\in[0.25,0.75],y\in[0.25,0.75]\} \tag{75}\]  
and the finite difference domain \(\Omega_{FDM}\), thus \(\Omega=\Omega_{FEM}\cup\Omega_{FDM}\), see Figure 1.  
The model problem in all our tests that is stated for the electric field \(E=(E_{1},E_{2})\) is as follows:  
\[ \begin{array}{ll}\varepsilon\partial_{tt}E+\
abla(\
abla\cdot E)-\triangle E -\
abla\
abla\cdot(\varepsilon E)=F&\text{in }\Omega\times(0,T),\\ E(\cdot,0)=0\text{ and }\partial_{t}E(\cdot,0)=0&\text{in }\Omega,\\ E=0&\text{on }\partial\Omega\times(0,T).\end{array}  \tag{76}\]  
We have the functions  
\[ \begin{array}{l}E_{1}=\frac{1}{\varepsilon}2\pi\sin^{2}\pi x\cos\pi y\sin \pi y\;\frac{t^{2}}{2},\\ E_{2}=-\frac{1}{\varepsilon}2\pi\sin^{2}\pi y\cos\pi x\sin\pi x\;\frac{t^{2}}{ 2},\end{array}  \tag{77}\]  
as the exact solution \(E=(E_{1},E_{2})\) of the model problem (76) with the source data \(F=(F_{1},F_{2})\) which corresponds to this exact solution. The function \(\varepsilon\) in (76) is defined as  
\[\varepsilon(x,y)=\left\{ \begin{array}{ll}1+\sin^{m}\pi(2x-0.5)\cdot\sin^{m} \pi(2y-0.5)&\text{in }\Omega_{FEM},\\ 1&\text{in }\Omega_{FDM}.\end{array} \right. \tag{78}\]  
We choose \(m=2,4,6,8\) in our numerical examples, see Figure 2, for these functions in the domain \(\Omega_{FEM}\). We note that the exact solution (77) satisfies the divergence-free condition \(\
abla\cdot(\varepsilon E)=0\) for \(\varepsilon\) defined by (78), the homogeneous initial conditions, as well as the homogeneous Dirichlet conditions for all times. The computational domain \(\Omega_{FEM}\times(0,T)\) was discretized into triangular elements with mesh sizes \(h_{l}=2^{-l},l=3,\ 4,\ 5,\ 6\), and the mesh in \(\Omega_{FDM}\times(0,T)\) was decomposed into squares of the same mesh sizes as described in Section 3, see Figure 1. The time step was chosen corresponding to the stability criterion (9) as \(\tau_{l}=0.025\cdot 2^{-l}\) for \(l=3,\ 4,\ 5,\ 6\). Convergence results of the proposed finite element scheme computed in \(L_{2}\) and \(H^{1}\) norms are presented in Tables 1-4 for \(m=2,4,6,8\) in (78). Relative norms in these tables were computed as  
\[ \begin{split} e_{l}^{1}&=\frac{\max_{1\leq k\leq N }\|E^{k}-E_{h}^{k}\|}{\max_{1\leq k\leq N}\|E^{k}\|},\\ e_{l}^{2}&=\frac{\max_{1\leq k\leq N}\|\
abla(E^{k} -E_{h}^{k})\|}{\max_{1\leq k\leq N}\|\
abla E^{k}\|}.\end{split}  \tag{79}\]  
\(E\) and \(E_{h}\) are the exact and computed FE solutions in \(\Omega_{FEM}\times(0,T)\), respectively, and \(N=T/\tau_{l}\). Logarithmic convergence rates \(r_{1},r_{2}\) in these tables are computed, viz. \[r_{1}=\frac{\left|\log\left(\frac{el_{h}^{1}}{el_{2h}^{1}}\right)\right|}{| \log(0.5)|},\quad\text{and}\quad r_{2}=\frac{\left|\log\left(\frac{el_{h}^{2} }{el_{2h}^{2}}\right)\right|}{|\log(0.5)|}, \tag{80}\]  
where \(el_{h}^{1,2},el_{2h}^{1,2}\) are relative norms computed via (79) on the mesh \(\mathcal{T}_{h}\) with the mesh size \(h\) and \(2h\), respectively. Figure 3 shows convergence of the relative \(L_{2}\)- and \(H^{1}\)-norms computed via (79) and compared with exact behavior of \(h\) and \(h^{2}\). Figure 2: Function \(\varepsilon(x,y)\) in the domain \(\Omega_{FEM}\) for different values of \(m\) in (78).