Long time error analysis of finite difference time domain methods for the nonlinear Klein-Gordon equation with weak nonlinearity

TL;DR

This paper derives long-time error bounds for finite difference time domain methods applied to the nonlinear Klein-Gordon equation with weak nonlinearity, providing guidelines for mesh size and time step in simulations.

Contribution

It establishes rigorous error bounds for FDTD methods over long times in the weak nonlinearity regime, with explicit dependence on the small parameter and validation through numerical tests.

Findings

Error bounds valid up to time O(1/ε^β) for 0 ≤ β ≤ 2

Error bounds depend explicitly on mesh size h, time step τ, and ε

Numerical results confirm the sharpness of the bounds

Abstract

We establish error bounds of the finite difference time domain (FDTD) methods for the long time dynamics of the nonlinear Klein-Gordon equation (NKGE) with a cubic nonlinearity, while the nonlinearity strength is characterized by $\varepsilon^2$ with $0 <\varepsilon \leq 1$ a dimensionless parameter. When $0 < \varepsilon \ll 1$, it is in the weak nonlinearity regime and the problem is equivalent to the NKGE with small initial data, while the amplitude of the initial data (and the solution) is at $O(\varepsilon)$. Four different FDTD methods are adapted to discretize the problem and rigorous error bounds of the FDTD methods are established for the long time dynamics, i.e. error bounds are valid up to the time at $O(1/\varepsilon^{\beta})$ with $0 \le \beta \leq 2$, by using the energy method and the techniques of either the cut-off of the nonlinearity or the mathematical induction to bound the numerical approximate solutions. In the error bounds, we pay particular attention to how error bounds depend explicitly on the mesh size $h$ and time step $\tau$ as well as the small parameter $\varepsilon\in (0,1]$, especially in the weak nonlinearity regime when $0 < \varepsilon \ll 1$. Our error bounds indicate that, in order to get ``correct'' numerical solutions up to the time at $O(1/\varepsilon^{\beta})$, the $\varepsilon$-scalability (or meshing strategy) of the FDTD methods should be taken as: $h = O(\varepsilon^{\beta/2})$ and $\tau = O(\varepsilon^{\beta/2})$. Extensive numerical results are reported to confirm our error bounds and to demonstrate that they are sharp.

Tables

Table 1. Table 1: Spatial errors of the CNFD ( 7 ) for the NKGE ( 5 ) with a = 0 𝑎 0 a=0 , b = 2 ​ π 𝑏 2 𝜋 b=2\pi , β = 0 𝛽 0 \beta=0 and ( 69 )

$e_{h,{\tau }_e}\left(t=1\right)$	$h_0=\pi /16$	$h_0/2$	$h_0/2^2$	$h_0/2^3$	$h_0/2^4$	$h_0/2^5$
${\epsilon }_0=1$	3.77E-2	9.65E-3	2.43E-3	6.09E-4	1.52E-4	3.84E-5
order	-	1.97	1.99	2.00	2.00	1.98
${\epsilon }_0/2$	3.33E-2	8.35E-3	2.09E-3	5.22E-4	1.31E-4	3.34E-5
order	-	2.00	2.00	2.00	1.99	1.97
${\epsilon }_0/2^2$	3.48E-2	8.74E-3	2.19E-3	5.47E-4	1.37E-4	3.50E-5
order	-	1.99	2.00	2.00	2.00	1.97
${\epsilon }_0/2^3$	3.55E-2	8.92E-3	2.23E-3	5.58E-4	1.40E-4	3.57E-5
order	-	1.99	2.00	2.00	1.99	1.97
${\epsilon }_0/2^4$	3.57E-2	8.97E-3	2.24E-3	5.61E-4	1.40E-4	3.59E-5
order	-	1.99	2.00	2.00	2.00	1.96

Table 2. Table 2: Temporal errors of the CNFD ( 7 ) for the NKGE ( 5 ) with a = 0 𝑎 0 a=0 , b = 2 ​ π 𝑏 2 𝜋 b=2\pi , β = 0 𝛽 0 \beta=0 and ( 69 )

$e_{h_e,\tau }\left(t=1\right)$	${\tau }_0=0.05$	${\tau }_0/2$	${\tau }_0/2^2$	${\tau }_0/2^3$	${\tau }_0/2^4$	${\tau }_0/2^5$
${\epsilon }_0=1$	3.27E-2	8.57E-3	2.19E-3	5.53E-4	1.39E-4	3.48E-5
order	-	1.93	1.97	1.99	1.99	2.00
${\epsilon }_0/2$	2.10E-2	5.45E-3	1.39E-3	3.49E-4	8.76E-5	2.20E-5
order	-	1.96	1.97	1.99	1.99	1.99
${\epsilon }_0/2^2$	1.84E-2	4.75E-3	1.21E-3	3.04E-4	7.63E-5	1.91E-5
order	-	1.95	1.97	1.99	1.99	2.00
${\epsilon }_0/2^3$	1.78E-2	4.59E-3	1.17E-3	2.94E-4	7.37E-5	1.85E-5
order	-	1.96	1.97	1.99	2.00	1.99
${\epsilon }_0/2^4$	1.77E-2	4.56E-3	1.16E-3	2.91E-4	7.31E-5	1.83E-5
order	-	1.96	1.97	2.00	1.99	2.00

Table 3. Table 3: Spatial errors of the CNFD ( 7 ) for the NKGE ( 5 ) with a = 0 𝑎 0 a=0 , b = 2 ​ π 𝑏 2 𝜋 b=2\pi , β = 1 𝛽 1 \beta=1 and ( 69 )

$e_{h,{\tau }_e}\left(t=1/\epsilon \right)$	$h_0=\pi /16$	$h_0/2$	$h_0/2^2$	$h_0/2^3$	$h_0/2^4$	$h_0/2^5$
${\epsilon }_0=1$	3.77E-2	9.65E-3	2.43E-3	6.09E-4	1.52E-4	3.84E-5
order	-	1.97	1.99	2.00	2.00	1.98
${\epsilon }_0/4$	7.31E-2	1.77E-2	4.38E-3	1.09E-3	2.74E-4	7.02E-5
order	-	2.05	2.01	2.01	1.99	1.96
${\epsilon }_0/4^2$	6.60E-1	1.71E-1	4.31E-2	1.08E-2	2.70E-3	6.91E-4
order	-	1.95	1.99	2.00	2.00	1.97
${\epsilon }_0/4^3$	2.78E+0	7.25E-1	1.80E-1	4.50E-2	1.13E-2	2.88E-3
order	-	1.94	2.01	2.00	1.99	1.97
${\epsilon }_0/4^4$	5.67E+0	8.48E-1	3.96E-1	1.10E-1	2.81E-2	7.22E-3
order	-	2.74	1.10	1.85	1.97	1.96

Table 4. Table 4: Temporal errors of the CNFD ( 7 ) for the NKGE ( 5 ) with a = 0 𝑎 0 a=0 , b = 2 ​ π 𝑏 2 𝜋 b=2\pi , β = 1 𝛽 1 \beta=1 and ( 69 )

$e_{h_e,\tau }\left(t=1/\epsilon \right)$	${\tau }_0=0.05$	${\tau }_0/2$	${\tau }_0/2^2$	${\tau }_0/2^3$	${\tau }_0/2^4$	${\tau }_0/2^5$
${\epsilon }_0=1$	3.27E-2	8.57E-3	2.19E-3	5.53E-4	1.39E-4	3.48E-5
order	-	1.93	1.97	1.99	1.99	2.00
${\epsilon }_0/4$	4.01E-2	9.95E-3	2.49E-3	6.22E-4	1.56E-4	3.89E-5
order	-	2.01	2.00	2.00	2.00	2.00
${\epsilon }_0/4^2$	3.45E-1	8.79E-2	2.21E-2	5.53E-3	1.38E-3	3.46E-4
order	-	1.97	1.99	2.00	2.00	2.00
${\epsilon }_0/4^3$	1.47E+0	3.69E-1	9.19E-2	2.29E-2	5.74E-3	1.43E-3
order	-	1.99	2.01	2.00	2.00	2.01
${\epsilon }_0/4^4$	8.58E-1	7.05E-1	2.20E-1	5.75E-2	1.45E-2	3.64E-3
order	-	0.28	1.68	1.94	1.99	1.99

Table 5. Table 5: Spatial errors of the CNFD ( 7 ) for the NKGE ( 5 ) with a = 0 𝑎 0 a=0 , b = 2 ​ π 𝑏 2 𝜋 b=2\pi , β = 2 𝛽 2 \beta=2 and ( 69 )

$e_{h,{\tau }_e}\left(t=1/{\epsilon }^2\right)$	$h_0=\pi /16$	$h_0/2$	$h_0/2^2$	$h_0/2^3$	$h_0/2^4$	$h_0/2^5$
${\epsilon }_0=1$	3.77E-2	9.65E-3	2.43E-3	6.09E-4	1.52E-4	3.84E-5
order	-	1.97	1.99	2.00	2.00	1.98
${\epsilon }_0/2$	3.98E-2	9.56E-3	2.39E-3	5.97E-4	1.49E-4	3.81E-5
order	-	2.06	2.00	2.00	2.00	1.97
${\epsilon }_0/2^2$	7.17E-1	1.82E-1	4.55E-2	1.14E-2	2.85E-3	7.27E-4
order	-	1.98	2.00	2.00	2.00	1.97
${\epsilon }_0/2^3$	2.78E+0	6.54E-1	1.58E-1	3.92E-2	9.78E-3	2.50E-3
order	-	2.09	2.05	2.01	2.00	1.97
${\epsilon }_0/2^4$	3.31E+0	1.78E+0	5.92E-1	1.55E-1	3.93E-2	1.01E-2
order	-	0.89	1.59	1.93	1.98	1.96

Table 6. Table 6: Temporal errors of the CNFD ( 7 ) for the NKGE ( 5 ) with a = 0 𝑎 0 a=0 , b = 2 ​ π 𝑏 2 𝜋 b=2\pi , β = 2 𝛽 2 \beta=2 and ( 69 )

$e_{h_e,\tau }\left(t=1/{\epsilon }^2\right)$	${\tau }_0=0.05$	${\tau }_0/2$	${\tau }_0/2^2$	${\tau }_0/2^3$	${\tau }_0/2^4$	${\tau }_0/2^5$
${\epsilon }_0=1$	3.27E-2	8.57E-3	2.19E-3	5.53E-4	1.39E-4	3.48E-5
order	-	1.93	1.97	1.99	1.99	2.00
${\epsilon }_0/2$	2.56E-2	6.32E-3	1.58E-3	3.94E-4	9.86E-5	2.47E-5
order	-	2.02	2.00	2.00	2.00	2.00
${\epsilon }_0/2^2$	3.91E-1	9.83E-2	2.46E-2	6.16E-3	1.54E-3	3.85E-4
order	-	1.99	2.00	2.00	2.00	2.00
${\epsilon }_0/2^3$	1.40E+0	3.32E-1	8.14E-2	2.03E-2	5.06E-3	1.26E-3
order	-	2.08	2.03	2.00	2.00	2.01
${\epsilon }_0/2^4$	1.81E+0	1.13E+0	3.16E-1	8.07E-2	2.03E-2	5.07E-3
order	-	0.68	1.84	1.97	1.99	2.00

Table 7. Table 7: Spatial errors of the CNFD ( 74 ) for the oscillatory NKGE ( 92 ) with d = 1 𝑑 1 d=1 , β = 1 𝛽 1 \beta=1 and ( 93 )

$e_{h,k_e}\left(s=1\right)$	$h_0=1/8$	$h_0/2$	$h_0/2^2$	$h_0/2^3$	$h_0/2^4$	$h_0/2^5$
${\epsilon }_0=1$	1.68E-2	4.26E-3	1.07E-3	2.68E-4	6.72E-5	1.76E-5
order	-	1.98	1.99	2.00	2.00	1.93
${\epsilon }_0/4$	5.60E-2	1.44E-2	3.63E-3	9.08E-4	2.27E-4	5.69E-5
order	-	1.96	1.99	2.00	2.00	2.00
${\epsilon }_0/4^2$	2.00E-1	5.68E-2	1.45E-2	3.63E-3	9.07E-4	2.27E-4
order	-	1.82	1.97	2.00	2.00	2.00
${\epsilon }_0/4^3$	4.83E-1	2.02E-1	5.70E-2	1.45E-2	3.63E-3	9.12E-4
order	-	1.26	1.83	1.97	2.00	1.99
${\epsilon }_0/4^4$	6.21E-1	4.86E-1	2.03E-1	5.74E-2	1.48E-2	3.97E-3
order	-	0.35	1.26	1.82	1.96	1.90

Table 8. Table 8: Temporal errors of the CNFD ( 74 ) for the oscillatory NKGE ( 92 ) with d = 1 𝑑 1 d=1 , β = 1 𝛽 1 \beta=1 and ( 93 )

$e_{h_e,k}\left(s=1\right)$	$k_0=0.025$	$k_0/4$	$k_0/4^2$	$k_0/4^3$	$k_0/4^4$	$k_0/4^5$
${\epsilon }_0=1$	4.11E-3	2.64E-4	1.66E-5	1.05E-6	7.82E-8	1E-8
order	-	1.98	2.00	1.99	1.87	-
${\epsilon }_0/4^{2/3}$	4.88E-2	3.24E-3	2.04E-4	1.28E-5	8.29E-7	6.48E-8
order	-	1.96	1.99	2.00	1.97	1.84
${\epsilon }_0/4^{4/3}$	4.98E-1	5.06E-2	3.23E-3	2.02E-4	1.28E-5	8.73E-7
order	-	1.65	1.98	2.00	1.99	1.94
${\epsilon }_0/4^{6/3}$	1.75E+0	5.18E-1	5.13E-2	3.23E-3	2.02E-4	1.28E-5
order	-	0.88	1.67	1.99	2.00	1.99
${\epsilon }_0/4^{8/3}$	1.93E+0	1.71E+0	5.27E-1	5.18E-2	3.24E-3	2.02E-4
order	-	0.09	0.85	1.67	2.00	2.00

Table 9. Table 9: Spatial errors of the CNFD ( 74 ) for the oscillatory NKGE ( 92 ) with d = 1 𝑑 1 d=1 , β = 2 𝛽 2 \beta=2 and ( 93 )

$e_{h,k_e}\left(s=1\right)$	$h_0=1/8$	$h_0/2$	$h_0/2^2$	$h_0/2^3$	$h_0/2^4$	$h_0/2^5$
${\epsilon }_0=1$	1.68E-2	4.26E-3	1.07E-3	2.68E-4	6.72E-5	1.76E-5
order	-	1.98	1.99	2.00	2.00	1.93
${\epsilon }_0/2$	5.64E-2	1.46E-2	3.66E-3	9.16E-4	2.30E-4	5.74E-5
order	-	1.95	2.00	2.00	2.00	2.00
${\epsilon }_0/2^2$	2.01E-1	5.71E-2	1.46E-2	3.65E-3	9.12E-4	2.28E-4
order	-	1.82	1.97	2.00	2.00	2.00
${\epsilon }_0/2^3$	4.83E-1	2.03E-1	5.71E-2	1.45E-2	3.64E-3	9.14E-4
order	-	1.25	1.83	1.98	1.99	1.99
${\epsilon }_0/2^4$	6.22E-1	4.86E-1	2.03E-1	5.74E-2	1.48E-2	3.97E-3
order	-	0.36	1.26	1.82	1.96	1.90

Table 10. Table 10: Temporal errors of the CNFD ( 74 ) for the oscillatory NKGE ( 92 ) with d = 1 𝑑 1 d=1 , β = 2 𝛽 2 \beta=2 and ( 93 )

$e_{h_e,k}\left(s=1\right)$	$k_0=0.025$	$k_0/4$	$k_0/4^2$	$k_0/4^3$	$k_0/4^4$	$k_0/4^5$
${\epsilon }_0=1$	4.11E-3	2.64E-4	1.66E-5	1.05E-6	7.82E-8	1E-8
order	-	1.98	2.00	1.99	1.87	-
${\epsilon }_0/4^{1/3}$	4.99E-2	3.31E-3	2.08E-4	1.31E-5	8.48E-7	9.37E-8
order	-	1.96	2.00	1.99	1.97	1.59
${\epsilon }_0/4^{2/3}$	5.03E-1	5.13E-2	3.28E-3	2.05E-4	1.29E-5	8.85E-7
order	-	1.65	1.98	2.00	2.00	1.93
${\epsilon }_0/4^{3/3}$	1.77E+0	5.21E-1	5.17E-2	3.26E-3	2.04E-4	1.29E-5
order	-	0.88	1.67	1.99	2.00	1.99
${\epsilon }_0/4^{4/3}$	1.93E+0	1.72E+0	5.28E-1	5.19E-2	3.25E-3	2.03E-4
order	-	0.08	0.85	1.67	2.00	2.00

Equations

$$\begin{array}{rl}& {\partial }_{tt}u(\mathbf{x},t)-\Delta u(\mathbf{x},t)+u(\mathbf{x},t)+{\epsilon }^2{u}^3(\mathbf{x},t)=0,\mathbf{x}\in {\mathbb{T}}^d,t>0,\\ & u(\mathbf{x},0)=\varphi (\mathbf{x}),{\partial }_{t}u(\mathbf{x},0)=\gamma (\mathbf{x}),\mathbf{x}\in {\mathbb{T}}^d.\end{array}$$

$$\begin{array}{rl}E(t):=& {\int }_{{\mathbb{T}}^d}\left[|{\partial }_{t}u(\mathbf{x},t){|}^2+|\nabla u(\mathbf{x},t){|}^2+|u(\mathbf{x},t){|}^2+\frac{{\epsilon }^2}{2}|u(\mathbf{x},t){|}^4\right]d\mathbf{x}\\ \equiv & {\int }_{{\mathbb{T}}^d}\left[|\gamma (\mathbf{x}){|}^2+|\nabla \varphi (\mathbf{x}){|}^2+|\varphi (\mathbf{x}){|}^2+\frac{{\epsilon }^2}{2}|\varphi (\mathbf{x}){|}^4\right]d\mathbf{x}:=E(0)=O(1),t\ge 0.\end{array}$$

$$\begin{array}{rl}& {\partial }_{tt}w(\mathbf{x},t)-\Delta w(\mathbf{x},t)+w(\mathbf{x},t)+w^3(\mathbf{x},t)=0,\mathbf{x}\in {\mathbb{T}}^d,t>0,\\ & w(\mathbf{x},0)=\epsilon \varphi (\mathbf{x}),{\partial }_{t}w(\mathbf{x},0)=\epsilon \gamma (\mathbf{x}),\mathbf{x}\in {\mathbb{T}}^d.\end{array}$$

$$\begin{array}{rl}\bar{E}(t):=& {\int }_{{\mathbb{T}}^d}\left[|{\partial }_{t}w(\mathbf{x},t){|}^2+|\nabla w(\mathbf{x},t){|}^2+|w(\mathbf{x},t){|}^2+\frac{1}{2}|w(\mathbf{x},t){|}^4\right]d\mathbf{x}={\epsilon }^{2}E(t)\\ \equiv & {\int }_{{\mathbb{T}}^d}\left[{\epsilon }^2|\gamma (\mathbf{x}){|}^2+{\epsilon }^2|\nabla \varphi (\mathbf{x}){|}^2+{\epsilon }^2|\varphi (\mathbf{x}){|}^2+\frac{{\epsilon }^4}{2}|\varphi (\mathbf{x}){|}^4\right]d\mathbf{x}:=\bar{E}(0)=O({\epsilon }^2).\end{array}$$

$$h=O({\epsilon }^{\beta /2})\text{mbox}and\tau =O({\epsilon }^{\beta /2}).$$

$$\begin{array}{rl}& {\partial }_{tt}u(x,t)-{\partial }_{xx}u(x,t)+u(x,t)+{\epsilon }^2{u}^3(x,t)=0,x\in \Omega =(a,b),t>0,\\ & u(x,0)=\varphi (x),{\partial }_{t}u(x,0)=\gamma (x),x\in \overline{\Omega }=[a,b],\end{array}$$

$$x_j:=a+jh,j=0,1,\dots ,M;t_n:=n\tau ,n=0,1,2,\dots .$$

$$\Vert u{\Vert }_{l^2}^2=h\sum _{j=0}^{M-1}|u_j{|}^2,\Vert {\delta }_x^{+}u{\Vert }_{l^2}^2=h\sum _{j=0}^{M-1}|{\delta }_x^{+}{u}_j{|}^2,\Vert u{\Vert }_{l^{\infty }}=\underset{0\le j\le M-1}{\max }|u_j|,(u,v)=h\sum _{j=0}^{M-1}{u}_j{v}_j,$$

$${\delta }_t^{+}{u}_j^n=\frac{u_j^{n+1}-u_j^n}{\tau },{\delta }_t^{-}{u}_j^n=\frac{u_j^n-u_j^{n-1}}{\tau },{\delta }_t^2{u}_j^n=\frac{u_j^{n+1}-2u_j^n+u_j^{n-1}}{{\tau }^2},$$

$${\delta }_x^{+}{u}_j^n=\frac{u_{j+1}^n-u_j^n}{h},{\delta }_x^{-}{u}_j^n=\frac{u_j^n-u_{j-1}^n}{h},{\delta }_x^2{u}_j^n=\frac{u_{j+1}^n-2u_j^n+u_{j-1}^n}{h^2}.$$

$${\delta }_t^2{u}_j^n-\frac{1}{2}{\delta }_x^2\left(u_j^{n+1}+u_j^{n-1}\right)+\frac{1}{2}\left(u_j^{n+1}+u_j^{n-1}\right)+{\epsilon }^{2}G\left(u_j^{n+1},u_j^{n-1}\right)=0,n\ge 1;$$

$${\delta }_t^2{u}_j^n-{\delta }_x^2{u}_j^n+\frac{1}{2}\left(u_j^{n+1}+u_j^{n-1}\right)+{\epsilon }^{2}G\left(u_j^{n+1},u_j^{n-1}\right)=0,n\ge 1;$$

$${\delta }_t^2{u}_j^n-\frac{1}{2}{\delta }_x^2\left(u_j^{n+1}+u_j^{n-1}\right)+\frac{1}{2}\left(u_j^{n+1}+u_j^{n-1}\right)+{\epsilon }^2{\left(u_j^n\right)}^3=0,n\ge 1;$$

$${\delta }_t^2{u}_j^n-{\delta }_x^2{u}_j^n+u_j^n+{\epsilon }^2{\left(u_j^n\right)}^3=0,j=0,1,\dots ,M-1,n\ge 1.$$

$$G(v,w)=\frac{F(v)-F(w)}{v-w},\forall v,w\in \mathbb{R},F(v)={\int }_0^v{s}^{3}ds=\frac{v^4}{4},v\in \mathbb{R}.$$

$$u_0^{n+1}=u_M^{n+1},u_{-1}^{n+1}=u_{M-1}^{n+1},n\ge 0;u_j^0=\varphi (x_j),j=0,1,\dots ,M,$$

$$u_j^1=\varphi (x_j)+\tau \gamma (x_j)+\frac{{\tau }^2}{2}\left[{\delta }_x^2\varphi (x_j)-\varphi (x_j)-{\epsilon }^2{\left(\varphi (x_j)\right)}^3\right],j=0,1,\dots ,M.$$

$${\sigma }_{\max }:=\underset{0\le n\le T_0{\epsilon }^{-\beta }/\tau }{\max }\Vert u^n{\Vert }_{l^{\infty }}^2.$$

$$0<\tau <2h/\sqrt{4-h^2},h>0,0<\epsilon \le 1.$$

$$0<\tau <2/\sqrt{{\epsilon }^2{\sigma }_{\max }-1},h>0,0<\epsilon \le 1.$$

$$0<\tau <2h/\sqrt{4+h^2(1+{\epsilon }^2{\sigma }_{\max })},h>0,0<\epsilon \le 1.$$

$$\begin{array}{rl}{E}^n& :=\Vert {\delta }_t^{+}{u}^n{\Vert }_{l^2}^2+\frac{1}{2}\sum _{k=n}^{n+1}\Vert {\delta }_x^{+}{u}^k{\Vert }_{l^2}^2+\frac{1}{2}\sum _{k=n}^{n+1}\Vert u^k{\Vert }_{l^2}^2+\frac{{\epsilon }^{2}h}{4}\sum _{j=0}^{M-1}\left[(u_j^n{)}^4+(u_j^{n+1}{)}^4\right]\equiv E^0.\end{array}$$

$$\begin{array}{rl}{\tilde{E}}^n& :=\Vert {\delta }_t^{+}{u}^n{\Vert }_{l^2}^2+h\sum _{j=0}^{M-1}({\delta }_x^{+}{u}_j^n)({\delta }_x^{+}{u}_j^{n+1})+\frac{1}{2}\sum _{k=n}^{n+1}\Vert u^k{\Vert }_{l^2}^2+\frac{{\epsilon }^{2}h}{4}\sum _{j=0}^{M-1}\left[(u_j^n{)}^4+(u_j^{n+1}{)}^4\right]\\ & \equiv {\tilde{E}}^0,n\ge 0.\end{array}$$

$$[[u]{]}_j^n=\frac{u_j^{n+1}+u_j^{n-1}}{2},j=0,1,\dots ,M,n\ge 1.$$

$$[[u]{]}^n=u^n+\frac{{\tau }^2}{2}{F}^n([[u]{]}^n),n\ge 1,$$

$$F_j^n(v)={\delta }_x^2{v}_j-\left[1+\frac{{\epsilon }^2}{2}(|u_j^{n-1}{|}^2+|2v_j-u_j^{n-1}{|}^2)\right]v_j,j=0,1,\dots ,M,n\ge 1.$$

$$K^n(v)=v-u^n-\frac{{\tau }^2}{2}{F}^n(v),v\in X_M,n\ge 1.$$

$$\begin{array}{rl}\left(K^n(v),v\right)& =\Vert v{\Vert }_{l^2}^2-(u^n,v)+\frac{{\tau }^2}{2}\left[\Vert {\delta }_x^{+}v{\Vert }_{l^2}^2+\Vert v{\Vert }_{l^2}^2+\frac{{\epsilon }^2}{2}\left(|u^{n-1}{|}^2+|2v-u^{n-1}{|}^2,v^2\right)\right]\\ & \ge \left(\Vert v{\Vert }_{l^2}-\Vert u^n{\Vert }_{l^2}\right)\Vert v{\Vert }_{l^2},n\ge 1,\end{array}$$

$$\underset{\Vert v{\Vert }_{l^2}\to \infty }{\lim }\frac{\left(K^n(v),v\right)}{\Vert v{\Vert }_{l^2}}=\infty ,n\ge 1.$$

$$\Vert u^n{\Vert }_{l^2}^2+\Vert {\delta }_x^{+}{u}^n{\Vert }_{l^2}^2\le 2E^n=2E^0,n\ge 0.$$

Full text

\ams

35L70, 65M06, 65M12, 65M15, 81-08

Dedicated to Professor Jie Shen on the occasion of his 60th birthday

1 Introduction

Consider the nonlinear Klein-Gordon equation (NKGE) with a cubic nonlinearity on a torus $\mathbb{T}^{d}\;(d=1,2,3)$ [23, 27, 36, 37] as

[TABLE]

Here $t$ is time, $\bf{x}\in{\mathbb{R}}^{d}$ is the spatial coordinates, $u:=u({\bf{x}},t)$ is a real-valued scalar field, $0<\varepsilon\leq 1$ is a dimensionless parameter, and $\phi({\bf{x}})$ and $\gamma({\bf{x}})$ are two given real-valued functions which are independent of $\varepsilon$. The NKGE is a relativistic (and nonlinear) version of the Schrödinger equation and it is widely used in quantum electrodynamics, particle and/or plasma physics to describe the dynamics of a spinless particle in some extra potential [36, 13, 4, 7, 22, 33, 34]. Provided that $u(\cdot,t)\in H^{1}(\mathbb{T}^{d})$ and $\partial_{t}u(\cdot,t)\in L^{2}(\mathbb{T}^{d})$, the NKGE (1) is time symmetric or time reversible and conserves the energy [5, 19], i.e.,

[TABLE]

We remark here that, when $0<\varepsilon\ll 1$, rescaling the amplitude of the wave function $u$ by introducing $w({\bf{x}},t)=\varepsilon u({\bf{x}},t)$, then the NKGE (1) with weak nonlinearity can be reformulated as the following NKGE with small initial data, while the amplitude of the initial data (and the solution) is at $O(\varepsilon)$:

[TABLE]

Again, the above NKGE (3) is time symmetric or time reversible and conserves the energy [5, 19], i.e.,

[TABLE]

In other words, the NKGE with weak nonlinearity and $O(1)$ initial data, i.e. (1), is equivalent to it with small initial data and $O(1)$ nonlinearity, i.e. (3). In the following, we only present numerical methods and their error bounds for the NKGE with weak nonlinearity. Extensions of the numerical methods and their error bounds to the NKGE with small initial data are straightforward.

There are extensive analytical results in the literature for the NKGE (1) (or (3)). For the existence of global classical solutions and almost periodic solutions as well as asymptotic behavior of solutions, we refer to [11, 12, 15, 42, 10, 40, 41] and references therein. For the Cauchy problem with small initial data (or weak nonlinearity), the global existence and asymptotic behavior of solutions were studied in different space dimensions and with different nonlinear terms [26, 25, 31, 35, 38]. Recently, more attentions have been devoted to analyzing the life-span of the solutions of the NKGE (3) [25, 32]. The results indicate that the life-span of a smooth solution to the NKGE (3) (or (1)) is at least up to time at $O(\varepsilon^{-2})$ [18, 16]. For more details related to this topic, we refer to [17, 21] and references therein.

For the numerical aspects of the NKGE (1) (or (3)), different numerical methods have been proposed and analyzed in the literatures [5, 14, 20, 44], including the finite difference time domain (FDTD) methods [5, 14, 20, 44], exponential wave integrator Fourier pseudospectral (EWI-FP) method [5, 6, 9], multiscale time integrator Fourier pseudospectral (MTI-FP) method [4], etc. In these results, the error bounds are normally valid up to the time at $O(1)$. Since the life-span of the solution of the NKGE (1) can be up to the time at $O(\varepsilon^{-2})$, it is a natural question to ask how the performance of a numerical method for (1) up to the time at $O(\varepsilon^{-2})$, i.e. long time error analysis. In other words, one has to establish error bounds of the numerical method for (1) up to the time at $O(\varepsilon^{-2})$ instead of the classical error bounds which are only valid up to the time at $O(1)$. The purpose of this paper is to carry out rigorous error analysis of four widely used FDTD methods for the NKGE (1) in the long time regime. In our error bounds, we pay particular attention to how the error bounds depend explicitly on the mesh size $h$ and time step $\tau$ as well as the small parameter $\varepsilon\in(0,1]$. In our numerical analysis, besides the standard technique of the energy method and the inverse inequality, we adapt the cut-off of the nonlinearity for the conservative methods, and resp., the mathematical induction for nonconservative methods, to obtain a priori bound of the numerical solution in the $l^{\infty}$ norm. Based on our rigorous error bounds, in order to obtain “correct” numerical approximations of the NKGE (1) (or (3)) up to the long time at $(\varepsilon^{-\beta})$ with $0\leq\beta\leq 2$ a fixed constant, the $\varepsilon$-scalability (or meshing strategy requirement) of the FDTD methods when $0<\varepsilon\ll 1$ is:

[TABLE]

As a by-product, by rescaling the time as $t\to t/\varepsilon^{\beta}$ with $0\leq\beta\leq 2$ in (1), then the problem (1) can be re-formulated as an oscillatory NKGE whose solution propagates waves with wavelength at $O(1)$ in space and $O(\varepsilon^{\beta})$ in time. The FDTD methods to (1) and their error bounds over long time can be extended straightforwardly to the oscillatory NKGE up to the time at $O(1)$. With the error bounds, the $\varepsilon$-scalability (or meshing strategy) of the FDTD methods for the oscillatory NKGE can be drawn.

The rest of the paper is organized as follows. In Section 2, different explicit/semi-implicit/implicit and conservative/nonconservative FDTD discretizations are presented for the NKGE (1) and their properties of the stability, conservation and solvability are analyzed. In Section 3, we establish rigorous error estimates of the FDTD methods for the NKGE (1) over long time dynamics. Extensive numerical results are reported in Section 4 to confirm our error bounds. In Section 5, we extend the FDTD methods and their error bounds to an oscillatory NKGE. Finally, some conclusions are drawn in Section 6. Throughout this paper, we adopt the notation $p\lesssim q$ to represent that there exists a generic constant $C>0$, which is independent of the mesh size $h$ and time step $\tau$ as well as $\varepsilon$ such that $|p|\leq Cq$.

2 FDTD methods and their analysis

In this section, we adapt four different FDTD methods to discretize the NKGE (1) and analyze their properties, such as stability, energy conservation and solvability. For simplicity of notations, we shall only present the numerical methods and their analysis for the NKGE (1) in one space dimension (1D). Thanks to tensor grids, generalizations to higher dimensions are straightforward and results remain valid with minor modifications. In 1D, consider the following NKGE

[TABLE]

with periodic boundary conditions.

2.1 FDTD methods

Choose the temporal step size $\tau:=\Delta t>0$ and the spatial mesh size $h:=\Delta x>0$, and denote $M=(b-a)/h$ being a positive integer and the grid points and time steps as:

[TABLE]

Denote $X_{M}=\{u=(u_{0},u_{1},\ldots,u_{M})^{T}|u_{j}\in\mathbb{R},j=0,1,2,\ldots,M,u_{0}=u_{M}\}$ and we always use $u_{-1}=u_{M-1}$ and $u_{M+1}=u_{1}$ if they are involved. The standard discrete $l^{2}$, semi-$H^{1}$ and $l^{\infty}$ norms and inner product in $X_{M}$ are defined as

[TABLE]

with $\delta^{+}_{x}u\in X_{M}$ defined as $\delta^{+}_{x}u_{j}=(u_{j+1}-u_{j})/h$ for $j=0,1,\ldots,M-1$.

Let $u^{n}_{j}$ be the numerical approximation of $u(x_{j},t_{n})$ for $j=0,1,\ldots,M$, $n\geq 0$ and denote the numerical solution at time $t=t_{n}$ as $u^{n}=(u_{0}^{n},u_{1}^{n},\ldots,u_{M}^{n})^{T}\in X_{M}$. We introduce the finite difference operators as

[TABLE]

[TABLE]

Here we consider four frequently used FDTD methods to discretize the NKGE (5):

I. The Crank-Nicolson finite difference (CNFD) method

[TABLE]

II. A semi-implicit energy conservative finite difference (SIFD1) method

[TABLE]

III. Another semi-implicit finite difference (SIFD2) method

[TABLE]

IV. The leap-frog finite difference (LFFD) method

[TABLE]

Here,

[TABLE]

The initial and boundary conditions in (5) are discretized as

[TABLE]

where the initial velocity $\gamma(x)$ is employed to update the first step $u^{1}$ by the Taylor expansion and the NKGE (5) as

[TABLE]

It is easy to check that the above FDTD methods are all time symmetric or time reversible, i.e. they are unchanged if interchanging $n+1\leftrightarrow n-1$ and $\tau\leftrightarrow-\tau$. In addition, the LFFD (10) is explicit and might be the simplest and the most efficient discretization for the NKGE (5) with the computational cost per time step at $O(M)$. The others are implicit schemes. Nevertheless, the CNFD (7) and SIFD1 (8) can be solved via either a direct solver or an iterative solver with the computational cost per time step depending on the solver, which is usually larger than $O(M)$, especially in two dimensions (2D) and three dimensions (3D). Meanwhile, the solution of the SIFD2 (9) can be explicitly updated in the Fourier space with $O(M{\rm ln}M)$ computational cost per time step, and such approach is valid in higher dimensions.

2.2 Stability, energy conservation and solvability

Let $T_{0}>0$ be a fixed constant and $0\leq\beta\leq 2$, and denote

[TABLE]

Following the von Neumann linear stability analysis of the classical FDTD methods for the NKGE in the nonrelativistic limit regime [5, 29], we can conclude the linear stability of the above FDTD methods for the NKGE (5) in the following lemma.

Lemma 2.1**.**

(linear stability) For the above FDTD methods applied to the NKGE (5) up to the time $t=T_{0}\varepsilon^{-\beta}$, we have:

(i) The CNFD (7) is unconditionally stable for any $h>0,\tau>0$ and $0<\varepsilon\leq 1$.

(ii) When $h\geq 2$, the SIFD1 (8) is unconditionally stable for any $h>0$ and $\tau>0$; and when $0<h<2$, this scheme is conditionally stable under the stability condition

[TABLE]

(iii) When $\sigma_{\rm max}\leq\varepsilon^{-2}$, the SIFD2 (9) is unconditionally stable for any $h>0$ and $\tau>0$; and when $\sigma_{\rm max}>\varepsilon^{-2}$, this scheme is conditionally stable under the stability condition

[TABLE]

(iv) The LFFD (10) is conditionally stable under the stability condition

[TABLE]

Remark 2.2**.**

The stability of schemes (9) - (10) is related to $\sigma_{\max}$, dependent on the boundedness of the $l^{\infty}$ norm of the numerical solution $u^{n}$ at the previous time step. The convergence estimates up to the previous time step could ensure such a bound in the $l^{\infty}$ norm, by making use of the inverse inequality, and such an error estimate could be recovered at the next time step, as given by the Theorems presented in Section 3.

For the CNFD (7) and SIFD1 (8), we can show that they conserve the energy in the discretized level with the proofs proceeding in the analogous lines as those in [5, 30, 37] and we omit the details here for brevity.

Lemma 2.3**.**

(energy conservation) For $n\geq 0$, the CNFD (7) conserves the discrete energy as

[TABLE]

Similarly, the SIFD1 (8) conserves the discrete energy as

[TABLE]

Based on Lemma 19, we can show the unique solvability of the CNFD (7) at each time step as follows.

Lemma 2.4**.**

(solvability of CNFD) For any given $u^{n},u^{n-1}$ ($n\geq 1$), the solution $u^{n+1}$ of the CNFD (7) is unique at each time step.

Proof 2.5**.**

Firstly, we prove the existence of the solution for the CNFD (7). To simplify the notations, we denote the grid function $[\![u]\!]^{n}\in X_{M}$ with

[TABLE]

For any $u^{n-1},u^{n},u^{n+1}\in X_{M}$, we rewrite the CNFD (7) as

[TABLE]

where $F^{n}:X_{M}\to X_{M}$ with

[TABLE]

Define a map $K^{n}:X_{M}\to X_{M}$ as

[TABLE]

It is obvious that $K^{n}$ ($n\geq 1$) is continuous from $X_{M}$ to $X_{M}$. Moreover, the fact

[TABLE]

implies

[TABLE]

Then, we can conclude that there exists a solution $v^{\ast}$ such that $K^{n}(v^{\ast})=0$ by applying the Brouwer fixed point theorem [2, 8, 28]. In other words, the CNFD (7) is solvable.

Now, we proceed to verify the uniqueness. From (18), we can get

[TABLE]

Hence, by employing the discrete Sobolev inequality [2, 39], we can obtain

[TABLE]

For any $v\in X_{M}$, we define a functional $S(v):X_{M}\to\mathbb{R}$ as

[TABLE]

It is easy to check that $S(v)$ is strictly convex with the gradient of it denoted as $\nabla S(v):=[\partial_{v_{0}}S(v),\ldots,$ $\partial_{v_{M}}S_{M}(v)]^{T}$ turning out to be

[TABLE]

By the strict convexity of $S(v)$, we can get the uniqueness of $\nabla S(v)=0$, which yields the uniqueness of $u^{n+1}\in X_{M}$ immediately. Thus, the proof is completed.

Remark 2.6**.**

The solvability of the SIFD1 (8) can be obtain similarly to the CNFD (7) in Lemma 2.3. There exists a unique solution for the SIFD2 due to the fact that it solves a linear system with a strictly diagonally dominant matrix. The solvability and uniqueness for (10) are straightforward since it is explicit.

3 Error estimates

In this section, we will establish error bounds of the FDTD methods.

3.1 Main results

Motivated by the analytical results in [26, 25, 31, 35, 38, 18, 16] and references therein, we make the following assumptions on the exact solution $u$ of the NKGE (5) up to the time $t=T_{0}/\varepsilon^{2}$:

[TABLE]

here $L^{\infty}=L^{\infty}([0,T_{0}/\varepsilon^{2}];L^{\infty})$ and $W^{m,\infty}_{p}=\{u\in W^{m,\infty}|\frac{\partial^{l}}{\partial x^{l}}u(a)=\frac{\partial^{l}}{\partial x^{l}}u(b),\quad 0\leq l<m\}$ for $m\geq 1$.

Denote $M_{0}=\sup_{\varepsilon\in(0,1]}\|u(x,t)\|_{L^{\infty}}$ and the grid ‘error’ function ${e}^{n}\in X_{M}(n\geq 0)$ as

[TABLE]

where $u^{n}\in X_{M}$ is the numerical approximation of the NKGE (5).

For the CNFD (7), we can establish the following error estimates (see its detailed proof in Section 3.2):

Theorem 3.1**.**

Under the assumption ($A$), there exist constants $h_{0}>0$ and $\tau_{0}>0$ sufficiently small and independent of $\varepsilon$, such that, for any $0<\varepsilon\leq 1$, when $0<h\leq h_{0}\varepsilon^{\beta/2}$ and $0<\tau\leq\tau_{0}\varepsilon^{\beta/2}$, we have the following error estimates for the CNFD (7) with (12) and (13)

[TABLE]

For the LFFD (10), the error estimates can be established as follows (see its detailed proof in Section 3.3):

Theorem 3.2**.**

Assume $\tau\leq\frac{1}{2}\min\{1,h\}$ and under the assumption ($A$), there exist constants $h_{0}>0$ and $\tau_{0}>0$ sufficiently small and independent of $\varepsilon$, such that for any $0<\varepsilon\leq 1$, when $0<h\leq h_{0}\varepsilon^{\beta/2}$ and $0<\tau\leq\tau_{0}\varepsilon^{\beta/2}$ and under the stability condition (17), we have the error estimates for the LFFD (10) with (12) and (13) as

[TABLE]

Similarly, for the SIFD1 (8) and SIFD2 (9), we have the following error estimates (their proofs are quite similar and thus they are omitted for brevity):

Theorem 3.3**.**

Assume $\tau\lesssim h$ and under the assumption ($A$), there exist constants $h_{0}>0$ and $\tau_{0}>0$ sufficiently small and independent of $\varepsilon$, such that for any $0<\varepsilon\leq 1$, when $0<h\leq h_{0}\varepsilon^{\beta/2}$, $0<\tau\leq\tau_{0}\varepsilon^{\beta/2}$ and under the stability condition (15), we have the following error estimates for the SIFD1 (8) with (12) and (13)

[TABLE]

Theorem 3.4**.**

Assume $\tau\lesssim h$ and under the assumption ($A$), there exist constants $h_{0}>0$ and $\tau_{0}>0$ sufficiently small and independent of $\varepsilon$, such that for any $0<\varepsilon\leq 1$, when $0<h\leq h_{0}\varepsilon^{\beta/2}$, $0<\tau\leq\tau_{0}\varepsilon^{\beta/2}$ and under the stability condition (16), we have the following error estimates for the SIFD2 (9) with (12) and (13)

[TABLE]

Remark 3.5**.**

In 2D with $d=2$ and 3D with $d=3$ cases, the above theorems are still valid under the technical conditions $0<h\lesssim\varepsilon^{\beta/2}\sqrt{C_{d}(h)}$ and $0<\tau\lesssim\varepsilon^{\beta/2}\sqrt{C_{d}(h)}$ where $C_{d}(h)=1/|\ln h|$ when $d=2$, and $C_{d}(h)=h^{1/2}$ when $d=3$.

Hence, the four FDTD methods studied here share the same spatial/temporal resolution capacity for the NKGE (5) up to the long time at $O(\varepsilon^{-\beta})$ with $0\leq\beta\leq 2$. In fact, given an accuracy bound $\delta_{0}>0$, the $\varepsilon$-scalability (or meshing strategy) of the FDTD methods should be taken as

[TABLE]

This implies that, in order to get “correct” numerical solution up to the time at $O(\varepsilon^{-1})$, one has to take the meshing strategy: $h=O(\varepsilon^{1/2})$ and $\tau=O(\varepsilon^{1/2})$; and resp., in order to get “correct” numerical solution up to the time at $O(\varepsilon^{-2})$, one has to take the meshing strategy: $h=O(\varepsilon)$ and $\tau=O(\varepsilon)$. These results are very useful for practical computations on how to select mesh size and time step such that the numerical results are trustable!

3.2 The proof of Theorem 3.1

For the CNFD (7), we establish the error estimates in Theorem 3.1. The key of the proof is to deal with the nonlinearity and overcome the main difficulty in uniformly bounding the numerical solution $u^{n}$, i.e., $\|u^{n}\|_{l^{\infty}}\lesssim 1$. Here, we adapt the cut-off technique which has been widely used in the literature[1, 2, 39], i.e., the nonlinearity is truncated to a global Lipschitz function with compact support.

Denote $B=(1+M_{0})^{2}$, choose a smooth function $\rho(\theta)\in C_{0}^{\infty}(\mathbb{R}^{+})$ and define

[TABLE]

then $F_{B}(\theta)$ has compact support and is smooth and global Lipschitz, i.e., there exists $C_{B}$ independent of $h$, $\tau$ and $\varepsilon$, such that

[TABLE]

Set $\hat{u}^{0}=u^{0}$, $\hat{u}^{1}=u^{1}$ and determine $\hat{u}^{n+1}\in X_{M}$ for $n\geq 1$ as follows

[TABLE]

In fact, $\hat{u}^{n}_{j}$ can be viewed as another approximation of $u(x_{j},t_{n})$ for $j=0,1,\ldots,M$ and $n\geq 0$. It is easy to verify that the scheme (38) is uniquely solvable for sufficiently small $\tau$ by using the properties of $\rho$ and standard techniques in Section 2. Define the corresponding ‘error’ function $\hat{e}^{n}\in X_{M}$ as

[TABLE]

and we can establish the following estimates:

Theorem 3.6**.**

Under the assumption (A), there exist constants $h_{0}>0$ and $\tau_{0}>0$ sufficiently small and independent of $\varepsilon$, such that for any $0<\varepsilon\leq 1$, when $0<h\leq h_{0}\varepsilon^{\beta/2}$ and $0<\tau\leq\tau_{0}\varepsilon^{\beta/2}$, we have the following error estimates

[TABLE]

We begin with the local truncation error $\hat{\xi}^{n}\in X_{M}$ of the scheme (38) given as

[TABLE]

The following estimates hold for $\hat{\xi}^{n}$.

Lemma 3.7**.**

Under the assumption (A), we have

[TABLE]

Proof 3.8**.**

Under the assumption (A), by applying the Taylor expansion to (41), it leads to

[TABLE]

Similarly, we have $|\delta^{+}_{x}\hat{\xi}^{0}_{j}|\lesssim h^{2}+\tau^{2}$ for $0\leq j\leq M-1$. These immediately imply (42).

Next, we control the nonlinear term as follows.

Lemma 3.9**.**

For $j=0,1,\ldots,M$ and $1\leq n\leq{T_{0}\varepsilon^{-\beta}}/{\tau}-1$, denote the error of the nonlinear term

[TABLE]

under the assumption (A), we have

[TABLE]

Proof 3.10**.**

Noticing (37) and (43), direct calculation for $j=0,1,\ldots,M$ and $1\leq n\leq{T_{0}\varepsilon^{-\beta}}/{\tau}-1$ leads to

[TABLE]

where the constant $C$ is independent of $h,\tau$ and $\varepsilon$. Under the assumption (A) and the properties of $F_{B}$, we have

[TABLE]

*which completes the proof. *

Now, we proceed to study the growth of the errors and verify Theorem 40. Subtracting (38) from (41), the error $\hat{e}^{n}\in X_{M}$ satisfies

[TABLE]

Define the ‘energy’ for the error vector $\hat{e}^{n}$ as

[TABLE]

It is easy to see that

[TABLE]

Proof 3.11**.**

(Proof of Theorem 40) When $n=0$, the estimates in (40) are obvious and the $n=1$ case is already verified in Lemma 3.1 for sufficiently small $0<\tau<\tau_{1}$ and $0<h<h_{1}$. Thus, we only need to prove (40) for $2\leq n\leq{T_{0}\varepsilon^{-\beta}}/{\tau}$.

Multiplying both sides of (47) by $h\left(\hat{e}^{n+1}_{j}-\hat{e}^{n-1}_{j}\right)$, summing up for $j$, noticing the fact $0\leq\beta\leq 2$ and making use of the Young’s inequality and Lemmas 42 &44, we derive

[TABLE]

Summing the above inequalities for time steps from 1 to $n$, there exists a constant $C>0$ such that

[TABLE]

Hence, the discrete Gronwall’s inequality suggests that there exists a constant $\tau_{2}>0$ sufficiently small, such that when $0<\tau\leq\tau_{2}$, the following holds

[TABLE]

Recalling $\|\hat{e}^{n+1}\|_{l^{2}}^{2}+\|\delta^{+}_{x}\hat{e}^{n+1}\|_{l^{2}}^{2}\leq 2\hat{S}^{n}$ when $0<\varepsilon\leq 1$, we can obtain the error estimate

[TABLE]

Finally, we estimate $\|\hat{u}^{n+1}\|_{l^{\infty}}$ for $1\leq n\leq{T_{0}\varepsilon^{-\beta}}/{\tau}-1$. The discrete Sobolev inequality implies

[TABLE]

Thus, there exist $h_{2}>0$ and $\tau_{3}>0$ sufficiently small, when $0<h\leq h_{2}\varepsilon^{\beta/2}$ and $0<\tau\leq\tau_{3}\varepsilon^{\beta/2}$, we obtain

[TABLE]

The proof is completed by choosing $h_{0}=\min\{h_{1},h_{2}\}$ and $\tau_{0}=\min\{\tau_{1},\tau_{2},\tau_{3}\}$.

Proof 3.12**.**

(Proof of Theorem 3.1) In view of the definition of $\rho$, Theorem 40 implies that (38) collapses to (7). By the unique solvability of the CNFD, $\hat{u}^{n}$ is identical to $u^{n}$. Thus, Theorem 3.1 is a direct consequence of Theorem 40.

3.3 The proof of Theorem 3.2.

For the LFFD (10), we establish the error estimates in Theorem 3.2. Throughout this section, the stability condition (17) is assumed. Here, we sketch the proof and omit those parts similar to the proof of Theorem 3.1 in Section 3.2.

Proof 3.13**.**

Denote the local truncation error as $\tilde{\xi}^{n}\in X_{M}$

[TABLE]

and the error of the nonlinear term as $\tilde{\eta}^{n}\in X_{M}$

[TABLE]

Similar to Lemma 42, under the assumption (A), we have

[TABLE]

The error equation for the LFFD (10) can be derived as

[TABLE]

We adapt the mathematical induction to prove Theorem 3.2, i.e. we want to demonstrate that there exist $h_{0}>0$ and $\tau_{0}>0$, such that, when $0<h<h_{0}$ and $0<\tau<\tau_{0}$, under the stability condition (17), the error bounds hold

[TABLE]

for all $0\leq n\leq{T_{0}\varepsilon^{-\beta}}/{\tau}$ and $0\leq\beta\leq 2$, where $C_{1}$, $\tau_{0}$ and $h_{0}$ will be classified later. For $n=0$, (60) is trivial. For $n=1$, the error equation (59) and the estimate (58) imply

[TABLE]

In view of the triangle inequality, discrete Sobolev inequality and the assumption (A), there exist $h_{1}>0$ and $\tau_{1}>0$ sufficiently small, when $0<h\leq h_{1}$ and $0<\tau\leq\tau_{1}$, we have

[TABLE]

In other words, (60) hold for $n=1$.

Now we assume that (60) is valid for all $0\leq n\leq m-1\leq{T_{0}\varepsilon^{-\beta}}/{\tau}-1$, then we need to show shat it is still valid when $n=m$. From (57), the error of the nonlinear term can be controlled as

[TABLE]

Define the ‘energy’ for the error vector $e^{n}(n=0,1,\ldots)$ as

[TABLE]

where

[TABLE]

Under the assumption $\tau\leq\frac{1}{2}\min\{1,h\}$, we have $1-\tau^{2}/2-\tau^{2}/h^{2}\geq\frac{1}{4}>0$. Since

[TABLE]

we can conclude that

[TABLE]

Similar to the proof in Section 3.2, there exists $\tau_{2}>0$ sufficiently small, when $0<\tau\leq\tau_{2}$,

[TABLE]

where $C_{5}$ depends on $T_{0}$ and the exact solution $u(x,t)$. Letting $n=m$, we have

[TABLE]

where $C_{6}$ depends on $T_{0}$ and the exact solution $u(x,t)$.

It remains to estimate $\|{u}^{m}\|_{l^{\infty}}$ for $n=m$. In fact, the discrete Sobolev inequality implies

[TABLE]

Thus, there exist $h_{2}>0$ and $\tau_{3}>0$ sufficiently small, when $0<h\leq h_{2}\varepsilon^{\beta/2}$ and $0<\tau\leq\tau_{3}\varepsilon^{\beta/2}$, we obtain

[TABLE]

Under the stability condition (17) and the choices of $h_{0}=\min\{h_{1},h_{2}\}$, $\tau_{0}=\min\{\tau_{1},\tau_{2},\tau_{3}\}$ and $C_{1}=\max\{C_{2},C_{6}\}$, the estimates in (60) are valid when $n=m$. Hence, the mathematical induction process is done and the proof of Theorem 3.2 is completed.

4 Numerical results

In this section, we present numerical results of the FDTD methods for the NKGE (5) up to the long time at $O(\varepsilon^{-\beta})$ with $0\leq\beta\leq 2$ to verify our error bounds. We only show numerical results for the CNFD (7) and the results for other FDTD methods are quite similar which are omitted for brevity. In the numerical experiments, we take $a=0$, $b=2\pi$ and choose the initial data as

[TABLE]

The ‘exact’ solution is obtained numerically by the exponential-wave integrator Fourier pseudospectral method [5, 19] with a very fine mesh size and a very small time step, e.g. $h_{e}=\pi/2^{15}$ and $\tau_{e}=10^{-5}$. Denote $u^{n}_{h,\tau}$ as the numerical solution at time $t=t_{n}$ obtained by a numerical method with mesh size $h$ and time step $\tau$. In order to quantify the numerical results, we define the error function as follows:

[TABLE]

Here we study the following three cases with respect to different $0\leq\beta\leq 2$:

Case I. Fixed time dynamics up to the time at $O(1)$, i.e., $\beta=0$;

Case II. Intermediate long time dynamics up to the time at $O(\varepsilon^{-1})$, i.e., $\beta=1$;

Case III. Long time dynamics up to the time at $O(\varepsilon^{-2})$, i.e., $\beta=2$.

We first test the spatial discretization errors at $t_{\varepsilon}=1/\varepsilon^{\beta}$ for different $0<\varepsilon\leq 1$. In order to do this, we fix the time step as $\tau_{e}=10^{-5}$ such that the temporal error can be ignored, and solve the NKGE (5) under different mesh size $h$. Tables 1, 3 and 5 depict the spatial errors for $\beta=0$, $\beta=1$ and $\beta=2$, respectively. Then we check the temporal errors at $t_{\varepsilon}=1/\varepsilon^{\beta}$ for different $0<\varepsilon\leq 1$ with different time step $\tau$ and a fine mesh size $h_{e}=\pi/2^{11}$ such that the spatial errors can be neglected. Tables 2, 4 and 6 show the temporal errors for $\beta=0$, $\beta=1$ and $\beta=2$, respectively.

From Tables 1-6 for the CNFD and additional similar numerical results for other FDTD methods not shown here for brevity, we can draw the following observations:

(i) For any fixed $\varepsilon=\varepsilon_{0}>0$ or $\beta=0$, the FDTD methods are uniformly second-order accurate in both spatial and temporal discretizations (cf. Tables 1 & 2 and the first rows in Tables 3-6), which agree with those results in the literature. (ii) In the intermediate long time regime, i.e. $\beta=1$, the second order convergence in space and time of the FDTD methods can be observed only when $0<h\lesssim\varepsilon^{1/2}$ and $0<\tau\lesssim\varepsilon^{1/2}$ (cf. upper triangles above the diagonals (corresponding to $h\sim\varepsilon^{1/2}$ and $\tau\sim\varepsilon^{1/2}$, and being labelled in bold letters) in Tables 3-4), which confirm our error bounds. (iii) In the long time regime, i.e. $\beta=2$, the second order convergence in space and time of the FDTD methods can be observed only when $0<h\lesssim\varepsilon$ and $0<\tau\lesssim\varepsilon$ (cf. upper triangles above the diagonals (corresponding to $h\sim\varepsilon$ and $\tau\sim\varepsilon$, and being labelled in bold letters) in Tables 5-6), which again confirm our error bounds. In summary, our numerical results confirm our rigorous error bounds and show that they are sharp.

5 Extension to an oscillatory NKGE

Introducing a rescaling in time by $s={\varepsilon}^{\beta}t$ with $0\leq\beta\leq 2$ and denoting $v({\bf x},s):=u({\bf x},s/\varepsilon^{\beta})=u({\bf x},t)$, we can reformulate the NKGE (1) into the following oscillatory NKGE

[TABLE]

Again, the oscillatory NKGE (71) is time symmetric or time reversible and conserves the energy [5, 19], i.e.,

[TABLE]

In fact, the long time dynamics of the NKGE (1) up to the time at $t=O(\varepsilon^{-\beta})$ is equivalent to the dynamics of the oscillatory NKGE (71) up to the fixed time at $s=O(1)$. Of course, the solution of of the NKGE (1) propagates waves with wavelength at $O(1)$ in both space and time, and wave speed in space at $O(1)$ too. On the contrary, the solution of the oscillatory NKGE (71) propagates waves with wavelength at $O(1)$ in space and $O(\varepsilon^{\beta})$ in time, and wave speed in space at $O(\varepsilon^{-\beta})$. To illustrate this, Figures 1 & 2 show the solutions $v(0,s)$ and $v(x,1)$, respectively, of the oscillatory NKGE (71) with $d=1$, $\mathbb{T}=(0,2\pi)$ and initial data (69) for different $0<\varepsilon\leq 1$ and $\beta$. We remark here that the oscillatory nature of the oscillatory NKGE (71) is quite different with that of the NKGE in the nonrelativistic limit regime. In fact, in the nonrelativistic limit regime of the NKGE [3, 4, 5, 7], the solution propagates waves with wavelength at $O(1)$ in space and $O(\varepsilon^{2})$ in time, and wave speed in space at $O(1)$!

In the following, we extend the FDTD methods and their error bounds for the NKGE (1) in previous sections to the oscillatory NKGE (71). Again, for simplicity of notations, the FDTD methods and their error bounds are only presented in 1D, and the results can be easily generalized to high dimensions with minor modifications. In addition, the proofs for the error bounds are quite similar to those in Sections 2&3, and thus they are omitted for brevity. We adopt similar notations as those used in Sections 2&3 except stated otherwise. In 1D, consider the following oscillatory NKGE

[TABLE]

with periodic boundary conditions.

5.1 FDTD methods

Choose the temporal step size $k:=\Delta s>0$ and denote time steps as $s_{n}:=nk$ for $n\geq 0$. Let $v^{n}_{j}$ be the numerical approximation of $v(x_{j},s_{n})$ for $j=0,1,\ldots,M$ and $n\geq 0$, and denote the numerical solution at time $s=s_{n}$ as $v^{n}$. Introduce the temporal finite difference operators as

[TABLE]

We consider the following four FDTD methods:

I. The Crank-Nicolson finite difference (CNFD) method

[TABLE]

II. A semi-implicit energy conservative finite difference (SIFD1) method

[TABLE]

III. Another semi-implicit finite difference (SIFD2) method

[TABLE]

IV. The Leap-frog finite difference (LFFD) method

[TABLE]

The initial and boundary conditions are discretized as

[TABLE]

Using the Taylor expansion and noticing (73), the first step $v^{1}\in X_{M}$ can be computed as

[TABLE]

In fact, if we take $k=\tau\varepsilon^{\beta}$ in the FDTD methods in this section, then they are consistent with those FDTD methods presented in Section 2. Thus they have the same solutions.

We remark here that, in practical computations, in order to uniformly bound the first step value $v^{1}\in X_{M}$ for $\varepsilon\in(0,1]$, in the above approximation (79), $k{\varepsilon^{-\beta}}$ and $k^{2}{\varepsilon^{-2\beta}}$ are replaced by $\sin({k}{\varepsilon^{-\beta}})$ and $k\sin({k}{\varepsilon^{-2\beta}})$, respectively [5, 8].

5.2 Stability and energy conservation

Denote

[TABLE]

Similar to Section 2, following the von Neumann linear stability analysis of the classical FDTD methods for the NKGE in the nonrelativistic limit regime [5, 29], we can conclude the linear stability of the above FDTD methods for oscillatory NKGE (73) up to the fixed time $s=T_{0}$ in the following lemma.

Lemma 5.1**.**

For the above FDTD methods applied to the oscillatory NKGE (73) up to the fixed time $s=T_{0}$, we have:

(i) The CNFD (74) is unconditionally stable for any $h>0,k>0$ and $0<\varepsilon\leq 1$.

(ii) When $h\geq 2$, the SIFD1 (75) is unconditionally stable for any $h>0$ and $k>0$; and when $0<h<2$, this scheme is conditionally stable under the stability condition

[TABLE]

(iii) When $\tilde{\sigma}_{\rm max}\leq\varepsilon^{-2}$, the SIFD2 (76) is unconditionally stable for any $h>0$ and $k>0$; and when $\tilde{\sigma}_{\rm max}>\varepsilon^{-2}$, this scheme is conditionally stable under the stability condition

[TABLE]

(iv) The LFFD (77) is conditionally stable under the stability condition

[TABLE]

For the CNFD (74) and SIFD1 (75), we have the following energy conservation properties:

Lemma 5.2**.**

The CNFD (74) conserves the discrete energy as

[TABLE]

Similarly, the SIFD1 (75) conserves the discrete energy as

[TABLE]

5.3 Main results

Again, motivated by the analytical results and the assumptions on the NKGE (5), we assume that the exact solution $v$ of the oscillatory NKGE (73) satisfies

[TABLE]

Define the grid ‘error’ function $\tilde{e}^{n}\in X_{M}(n\geq 0)$ as

[TABLE]

where $v^{n}\in X_{M}$ is the numerical approximation of the oscillatory NKGE (73) obtained by one of the FDTD methods.

By taking $k=\tau\varepsilon^{\beta}$ in the above FDTD methods and noting the error bounds in Section 3, we can immediately obtain error bounds of the above FDTD methods for the oscillatory NKGE (73).

Theorem 5.3**.**

Under the assumption (B), there exist constants $h_{0}>0$ and $k_{0}>0$ sufficiently small and independent of $\varepsilon$, such that for any $0<\varepsilon\leq 1$, when $0<h\leq h_{0}\varepsilon^{\beta/2}$ and $0<k\leq k_{0}\varepsilon^{3\beta/2}$, we have the following error estimates for the CNFD (74) with (78) and (79)

[TABLE]

Theorem 5.4**.**

Assume $k\lesssim h\varepsilon^{\beta}$ and under the assumption (B), there exist constants $h_{0}>0$ and $k_{0}>0$ sufficiently small and independent of $\varepsilon$, such that for any $0<\varepsilon\leq 1$, when $0<h\leq h_{0}\varepsilon^{\beta/2}$, $0<k\leq k_{0}\varepsilon^{3\beta/2}$ and under the stability condition (81), we have the following error estimates for the SIFD1 (75) with (78) and (79)

[TABLE]

Theorem 5.5**.**

Assume $k\lesssim h\varepsilon^{\beta}$ and under the assumption (B), there exist constants $h_{0}>0$ and $k_{0}>0$ sufficiently small and independent of $\varepsilon$, such that for any $0<\varepsilon\leq 1$, when $0<h\leq h_{0}\varepsilon^{\beta/2}$, $0<k\leq k_{0}\varepsilon^{3\beta/2}$ and under the stability condition (82), we have the following error estimates for the SIFD2 (76) with (78) and (79)

[TABLE]

Theorem 5.6**.**

Assume $k\lesssim h\varepsilon^{\beta}$ and under the assumption (B), there exist constants $h_{0}>0$ and $k_{0}>0$ sufficiently small and independent of $\varepsilon$, such that for any $0<\varepsilon\leq 1$, when $0<h\leq h_{0}\varepsilon^{\beta/2}$, $0<k\leq k_{0}\varepsilon^{3\beta/2}$ and under the stability condition (83), we have the following error estimates for the LFFD (77) with (78) and (79)

[TABLE]

The above four FDTD methods share the same spatial/temporal resolution capacity for the oscillatory NKGE (73) up to the fixed time at $O(1)$. In fact, given an accuracy bound $\delta_{0}>0$, the $\varepsilon$-scalability (or meshing strategy) of the FDTD methods for the oscillatory NKGE (73) should be taken as

[TABLE]

Again, these results are very useful for practical computations on how to select mesh size and time step such that the numerical results are trustable!

5.4 Numerical results of the oscillatory NKGE in the whole space

Consider the following oscillatory NKGE in $d$-dimensional ($d=1,2,3$) whole space

[TABLE]

Similar to the oscillatory NKGE (71), the solution of of the oscillatory NKGE (92) propagates waves with wavelength at $O(1)$ in space and $O(\varepsilon^{\beta})$ in time, and wave speed in space at $O(\varepsilon^{-\beta})$. To illustrate the rapid wave propagation in space at $O(\varepsilon^{-\beta})$, Figure 3 shows the solution $v(x,1)$ of the oscillatory NKGE (92) with $d=1$ and initial data

[TABLE]

Similar to those in the literature, by using the fast decay of the solution of the oscillatory NKGE (92) at the far field (see [5, 20, 37] and references therein), in practical computation, we usually truncate the originally whole space problem onto a bounded domain $\Omega$ with periodic boundary conditions, provided that $\Omega$ is large enough such that the truncation error is negligible. Then the truncated problem can be solved by the FDTD methods. Of course, due to the rapid wave propagation in space of the oscillatory NKGE (92) (cf. Fig. 3), in order to compute numerical solution up to the time at $O(1)$, in general, the size of the bounded domain $\Omega$ has to be taken as $O(\varepsilon^{-\beta})$.

In the following, we report numerical results of the oscillatory NKGE (92) with $d=1$. The initial data is chosen as (93) and the bounded computational domain is taken as $\Omega_{\varepsilon}=[-4-{\varepsilon^{-\beta}},4+{\varepsilon^{-\beta}}]$. The ‘exact’ solution is obtained numerically by the exponential-wave integrator Fourier pseudospectral method with a very fine mesh size and a very small time step, e.g. $h_{e}=1/2^{13}$ and $k_{e}=2\times 10^{-6}$. Denote $v^{n}_{h,k}$ as the numerical solution at $s=s_{n}$ obtained by a numerical method with mesh size $h$ and time step $k$. In order to quantify the numerical results, we define the error function as follows:

[TABLE]

Tables 7 and 8 show the spatial and temporal errors, respectively, of the CNFD method with $\beta=1$, and Tables 9 and 10 show similar results for $\beta=2$. The results for other FDTD methods are quite similar and they are omitted here for brevity.

From Tables 7-10 for the CNFD and additional similar numerical results for other FDTD methods not shown here for brevity, we can draw the following observations on the FDTD methods for the oscillatory NKGE (71) (or (92)):

(i) For any fixed $\varepsilon=\varepsilon_{0}>0$ or $\beta=0$, the FDTD methods are uniformly second-order accurate in both spatial and temporal discretizations (cf. the first rows in Tables 7-10), which agree with those results in the literature. (ii) In the intermediate oscillatory case, i.e. $\beta=1$, the second order convergence in space and time of the FDTD methods can be observed only when $0<h\lesssim\varepsilon^{1/2}$ and $0<k\lesssim\varepsilon^{3/2}$ (cf. upper triangles above the diagonals (corresponding to $h\sim\varepsilon^{1/2}$ and $k\sim\varepsilon^{3/2}$, and being labelled in bold letters) in Tables 7-8), which confirm our error bounds. (iii) In the highly oscillatory case, i.e. $\beta=2$, the second order convergence in space and time of the FDTD methods can be observed only when $0<h\lesssim\varepsilon$ and $0<k\lesssim\varepsilon^{3}$ (cf. upper triangles above the diagonals (corresponding to $h\sim\varepsilon$ and $k\sim\varepsilon^{3}$, and being labelled in bold letters) in Tables 9-10), which again confirm our error bounds. In summary, our numerical results confirm our rigorous error bounds and show that they are sharp.

6 Conclusion

Four different finite difference time domain FDTD methods were adapted to discretize the nonlinear Klein-Gordon equation (NKGE) with a weak cubic nonlinearity, while the nonlinearity strength is characterized by $\varepsilon^{2}$ with $0<\varepsilon\leq 1$ a dimensionless parameter. Rigorous error estimates were established for the long time dynamics of the NKGE up to the time at $O(\varepsilon^{-\beta})$ with $0\leq\beta\leq 2$. The error bounds depend explicitly on the mesh size $h$ and time step $\tau$ as well as the small parameter $\varepsilon\in(0,1]$, which indicate the temporal and spatial resolution capacities of the FDTD methods for the long time dynamics of the NKGE. Based on the error bounds, in order to get “correct” numerical solution of the NKGE up to the long time at $O(\varepsilon^{-\beta})$ with $0<\beta\leq 2$, the $\varepsilon$-scalability (or meshing strategy) of the FDTD methods has to be taken as: $h=O(\varepsilon^{\beta/2})$ and $\tau=O(\varepsilon^{\beta/2})$. In addition, the FDTD methods were also applied to solve an oscillatory NKGE and their error bounds were also obtained. Extensive numerical results were reported to confirm our error bounds and to demonstrate that they are sharp.

Acknowledgments

We thank fruitful discussion with Dr Chunmei Su. This work was partially supported by the Ministry of Education of Singapore grant R-146-000-223-112.