Exact finite difference schemes for three-dimensional linear systems with constant coefficients

TL;DR

This paper develops implicit and explicit exact difference schemes for three-dimensional linear systems with constant coefficients, demonstrating their efficiency and accuracy through numerical simulations, and extends the approach to more general systems.

Contribution

It introduces new exact difference schemes for 3D linear systems and shows their effectiveness compared to high-order methods, with potential extensions to nonlinear systems.

Findings

Numerical simulations confirm efficiency for stiff problems.

Exact difference schemes outperform high-order methods.

Potential extension to nonlinear systems with stability preservation.

Abstract

In this paper implicit and explicit exact difference schemes (EDS) for system $\textbf{x}' = A\textbf{x}$ of three linear differential equations with constant coefficients are constructed. Numerical simulations for stiff problem and for problems with periodic solutions on very large time interval demonstrate the efficiency and exactness of the EDS compared with high-order numerical methods. This result can be extended for constructing EDS for general systems of $n$ linear differential equations with constant coefficients and nonstandard finite difference (NSFD) schemes preserving stability properties for quasi-linear system of equations $\textbf{x}' = A\textbf{x }+ f(\textbf{x})$.

Tables

Table 1. Table 1 : 3 × 3 3 3 3\times 3 Jordan form matrices

$\sigma \left(A\right)$	${\chi }_A\left(t\right)$	$m_A\left(t\right)$	$J$
$\{{\lambda }_1,{\lambda }_2,{\lambda }_3\}$	$\left(t-{\lambda }_1\right)\left(t-{\lambda }_2\right)\left(t-{\lambda }_3\right)$	$\left(t-{\lambda }_1\right)\left(t-{\lambda }_2\right)\left(t-{\lambda }_3\right)$
$\{{\lambda }_1,{\lambda }_1,{\lambda }_2\}$	${\left(t-{\lambda }_1\right)}^2\left(t-{\lambda }_2\right)$	$\left(t-{\lambda }_1\right)\left(t-{\lambda }_2\right)$
$\{{\lambda }_1,{\lambda }_1,{\lambda }_2\}$	${\left(t-{\lambda }_1\right)}^2\left(t-{\lambda }_2\right)$	${\left(t-{\lambda }_1\right)}^2\left(t-{\lambda }_2\right)$
$\{\lambda ,\lambda ,\lambda \}$	${\left(t-\lambda \right)}^3$	$\left(t-\lambda \right)$
$\{\lambda ,\lambda ,\lambda \}$	${\left(t-\lambda \right)}^3$	${\left(t-\lambda \right)}^2$
$\{\lambda ,\lambda ,\lambda \}$	${\left(t-\lambda \right)}^3$	${\left(t-\lambda \right)}^3$

Table 2. Table 2: The exact difference schemes for x ′ = A ​ x superscript x ′ 𝐴 x \textbf{x}^{\prime}=A\textbf{x} in dependence on the eigenvalues of A 𝐴 A

Case	$\sigma \left(A\right)$	$\psi ,\varphi ,\theta$ are determined by
$1$	$\{{\lambda }_1,{\lambda }_2,{\lambda }_3\}$	Theorem 3.1
Special case $1$	$\{0,{\lambda }_2,{\lambda }_3\}$	Theorem 3.2
Special case $2$	$\{\alpha +\beta i,\alpha -\beta i,\lambda \}$	Corollary 1
$2$	$\{{\lambda }_1,{\lambda }_1,{\lambda }_2\}$	Theorem 3.4
$3$	$\{\lambda ,\lambda ,\lambda \}$	Theorem 3.7

Table 3. Table 3: Error of the methods

$T,\lambda$	$h$	IEDS error	EEDS error	RK4 error	Taylor error	Trapezoidal error
$1,1$	${10}^{-5}$	3.2618e-011	3.8608e-011	3.1419e-014	3.4646e-011	1.0599e-010
	${10}^{-4}$	1.5561e-012	1.2415e-012	1.2990e-014	1.5561e-012	3.4218e-009
	${10}^{-3}$	4.9460e-013	8.1424e-013	3.3751e-014	1.2468e-013	3.4167e-007
	${10}^{-2}$	4.5852e-014	4.6851e-014	3.4000e-010	2.9421e-014	3.4167e-005
	${10}^{-1}$	3.2196e-015	7.7716e-015	3.2526e-006	7.7251e-010	0.0034
	$1$	7.7716e-016	1.1102e-016	0.0195	4.4648e-004	0.3829
$10,{10}^{-1}$	$h={10}^{-5}$	1.0909e-010	4.9326e-011	1.3178e-013	1.0909e-010	2.1743e-010
	$h={10}^{-4}$	3.0020e-011	4.3151e-011	8.3267e-015	3.7715e-011	1.1621e-008
	$h={10}^{-3}$	1.3045e-012	1.0316e-012	1.3023e-013	1.3045e-012	1.1548e-006
	$h={10}^{-2}$	7.4307e-013	1.6542e-013	1.1505e-009	1.3267e-013	1.1548e-004
	$h={10}^{-1}$	3.7637e-014	3.7970e-014	1.1280e-005	2.6824e-009	0.0115
	$h=1$	2.9976e-015	3.4417e-015	0.0995	0.0025	0.8470
	$h=10$	1.3323e-015	1.4433e-015	524.6383	1.6976e+003	1.2944
${10}^2,{10}^{-2}$	$h={10}^{-4}$	1.3930e-010	1.0862e-010	8.1490e-014	1.3930e-010	1.1415e-007
	$h={10}^{-3}$	3.5170e-011	4.5318e-011	1.1456e-012	3.5170e-011	1.1406e-005
	$h={10}^{-2}$	2.5135e-012	2.1225e-012	1.1430e-008	1.3967e-012	0.0011
	$h={10}^{-1}$	1.9151e-013	7.0943e-013	1.1612e-004	2.7668e-008	0.1150
	$h=1$	3.3640e-014	6.0507e-014	0.6364	0.0230	1.3021
	$h=10$	1.7208e-014	9.6589e-015	1.4626e+026	1.0099e+031	2.7849
	$h={10}^2$	1.1102e-016	3.3307e-016	4.3282e+006	1.4679e+009	2.6896
${10}^3,{10}^{-3}$	$h={10}^{-4}$	1.1436e-010	6.5578e-009	1.7082	4.8609e-009	1.1592e-006
	$h={10}^{-3}$	1.1436e-010	7.7965e-010	1.1645e-011	1.1436e-010	1.1577e-004
	$h={10}^{-2}$	3.9845e-011	5.0336e-011	1.1595e-007	3.6106e-011	0.0116
	$h={10}^{-1}$	2.0014e-012	7.6292e-012	0.0012	2.7916e-007	1.1135
	$h=1$	7.6816e-013	7.4385e-014	1.3873	0.2445	2.7559
	$h=10$	1.0358e-013	6.9056e-014	2.0312e+260	Inf	2.5752
	$h={10}^2$	3.2196e-015	4.4409e-015	2.0587e+066	3.6683e+091	1.5772
	$h={10}^3$	5.5511e-016	4.4409e-016	4.1833e+010	1.3972e+015	2.6670
${10}^4,{10}^{-4}$	$h={10}^{-3}$	2.1401e-009	4.7583e-009	1.6131	4.5953e-009	0.0010
	$h={10}^{-2}$	1.5582e-010	1.6500e-010	1.0436e-006	1.2046e-010	0.1024
	$h={10}^{-1}$	3.3033e-011	3.9898e-011	0.0100	2.3738e-006	2.4012
	$h=1$	6.3882e-012	1.0866e-011	1.2578	5.2064	2.0189
	$h=10$	7.4413e-013	1.0292e-013	NaN	NaN	2.6229
	$h={10}^2$	3.3529e-014	5.6566e-014	NaN	NaN	1.3602
	$h={10}^3$	6.0507e-015	3.4417e-015	1.6389e+106	2.8260e+151	2.2193
	$h={10}^4$	1.6653e-016	1.1102e-016	4.1683e+014	1.3897e+021	0.6356
${10}^5,{10}^{-5}$	$h={10}^{-2}$	8.3200e-009	2.8834e-009	8.6926e-006	4.6617e-009	1.0818
	$h={10}^{-1}$	2.3169e-010	9.1972e-011	0.0952	2.2129e-005	1.9557
	$h=1$	3.2853e-011	4.3130e-011	1.0351	3.8828e+006	1.0811
	$h=10$	2.0601e-011	7.6230e-012	NaN	NaN	1.4332
	$h={10}^2$	3.2153e-013	2.0207e-013	NaN	NaN	1.1208
	$h={10}^3$	5.5303e-014	5.1750e-014	NaN	NaN	2.3456
	$h={10}^4$	3.4431e-014	5.6760e-015	1.5835e+146	2.6870e+211	2.0414
	$h={10}^5$	4.9544e-015	1.1102e-016	4.1668e+018	1.3890e+027	0.3181

Table 4. Table 4: Error of the methods

$T$	$h$	IEDS error	EEDS error	RK4 error	Taylor error	Radau IIA error
${10}^{-3}$	${10}^{-6}$	6.6613e-016	6.6613e-016	2.4425e-015	2.7756e-014	2.4425e-015
	${10}^{-5}$	5.5511e-016	6.6613e-016	1.1102e-015	8.5487e-015	6.6613e-016
	${10}^{-4}$	5.5511e-016	3.3307e-016	7.6037e-012	5.5511e-016	1.5543e-015
	${10}^{-3}$	2.2204e-016	2.2204e-016	8.1964e-008	1.9596e-011	1.2359e-010
${10}^{-2}$	${10}^{-6}$	1.1102e-015	6.6613e-016	8.6597e-015	2.6584e-013	8.6597e-015
	${10}^{-5}$	1.0547e-015	6.6613e-016	6.6613e-015	3.8858e-014	3.3307e-015
	${10}^{-4}$	8.8818e-016	4.9960e-016	3.0913e-011	3.0531e-015	5.7732e-015
	${10}^{-3}$	6.1062e-016	4.9960e-016	3.3324e-007	7.9673e-011	5.0249e-010
	${10}^{-2}$	2.7756e-016	2.2204e-016	0.0071	1.7611e-004	4.5087e-005
${10}^{-1}$	${10}^{-6}$	8.5165e-015	2.9616e-015	3.4529e-014	2.2968e-012	3.4418e-014
	${10}^{-5}$	7.8753e-015	2.6691e-015	1.1999e-014	6.8204e-014	1.1991e-014
	${10}^{-4}$	7.2122e-015	2.7515e-015	3.0913e-011	9.2176e-015	5.7732e-015
	${10}^{-3}$	4.0069e-015	1.8644e-015	3.3324e-007	7.9673e-011	5.0249e-010
	${10}^{-2}$	3.6078e-015	1.2257e-015	0.0071	1.7611e-004	4.5087e-005
	${10}^{-1}$	3.6078e-015	4.3819e-016	291.0000	846.5555	0.0517
1	${10}^{-5}$	4.5214e-014	7.6050e-015	1.3323e-014	2.4566e-013	1.3212e-014
	${10}^{-4}$	4.1633e-014	7.3841e-015	3.0913e-011	2.3925e-014	5.7732e-015
	${10}^{-3}$	4.1633e-014	7.2164e-015	3.3324e-007	7.9673e-011	5.0249e-010
	${10}^{-2}$	2.3564e-014	4.7699e-015	0.0071	1.7611e-004	4.5087e-005
	${10}^{-1}$	1.6376e-014	3.7192e-015	4.3544e+024	1.8904e+029	0.0517
	$1$	5.2180e-015	1.1102e-016	4.0049e+006	1.3096e+009	0.0264

Equations

$$\textbf{x}^{\prime}(t)=A\textbf{x}(t),\quad\textbf{x}(t)=\big{(}x(t),y(t)\big{)}^{T},\qquad A\in M_{2\times 2}(\mathbb{R})$$

$$\dfrac{\textbf{x}_{k+1}-\textbf{x}_{k}}{\phi(h)}=A\big{[}\theta\textbf{x}_{k+1}+(1-\theta)\textbf{x}_{k}\big{]},$$

$$\frac{du}{dt}=au+bw,\frac{dw}{dt}=cu+dw,$$

$$\frac{u_{k+1}-\psi u_k}{\varphi }=a{u}_k+b{w}_k,\frac{w_{k+1}-\psi w_k}{\varphi }=c{u}_k+d{w}_k,$$

$$\psi =\frac{{\lambda }_1{e}^{{\lambda }_{2}h}-{\lambda }_2{e}^{{\lambda }_{1}h}}{{\lambda }_1-{\lambda }_2}=1+\mathcal{O}(h^2),\varphi =\frac{e^{{\lambda }_{1}h}-e^{{\lambda }_{2}h}}{{\lambda }_1-{\lambda }_2}=h+\mathcal{O}(h^2),$$

$$\det \left(\begin{array}{cc}a-\lambda & b\\ c& d-\lambda \end{array}\right)=0.$$

$$\textbf{x}^{\prime}(t)=A\textbf{x}(t),\quad\textbf{x}(t)=\big{(}x(t),y(t),z(t)\big{)}^{T},\quad t\in[0,T],\qquad A\in M_{3\times 3}(\mathbb{R}).$$

$$\dfrac{\textbf{x}_{k+1}-\psi(h)\textbf{x}_{k}}{\phi(h)}=A\big{[}\theta\textbf{x}_{k+1}+(1-\theta)\textbf{x}_{k}\big{]},$$

$$\frac{{\text{x}}_{k+1}-{\text{x}}_k}{h}=A{\text{x}}_k+{\alpha }_{2}h{A}^2{\text{x}}_k+{\alpha }_3{h}^2{A}^3{\text{x}}_k+\dots +{\alpha }_s{h}^{s-1}{A}^s{\text{x}}_k,$$

$$\frac{{\text{x}}_{k+1}-\psi (h){\text{x}}_k}{\varphi }=A{\text{x}}_k+{\alpha }_2\varphi A^2{\text{x}}_k+{\alpha }_3{\varphi }^2{A}^3{\text{x}}_k+\dots +{\alpha }_s{\varphi }^{s-1}{A}^s{\text{x}}_k.$$

$$\frac{{\text{x}}_{k+1}-\psi (h){\text{x}}_k}{\varphi (h)}=A{\text{x}}_{k+1}+{\alpha }_2\varphi A^2{\text{x}}_{k+1}+{\alpha }_3{\varphi }^2{A}^3{\text{x}}_{k+1}+\dots +{\alpha }_s{\varphi }^{s-1}{A}^s{\text{x}}_{k+1}.$$

$$\frac{{\text{x}}_{k+1}-\psi (h){\text{x}}_k}{\varphi (h)}=A[{\theta }_1{\text{x}}_k+(1-{\theta }_1){\text{x}}_{k+1}]+\sum _{m=2}^s{\alpha }_m{\varphi }^{m-1}{A}^m[{\theta }_m{\text{x}}_k+(1-{\theta }_m){\text{x}}_{k+1}].$$

$$\frac{{\text{x}}_{k+1}-\psi (h){\text{x}}_k}{\varphi }=A{\text{x}}_k+\theta \varphi A^2{\text{x}}_k.$$

$$\dfrac{\textbf{u}_{k+1}-\psi(h)\textbf{u}_{k}}{\phi(h)}=J\big{[}\theta\textbf{u}_{k+1}+(1-\theta)\textbf{u}_{k}\big{]}$$

$$J_1=\left(\begin{array}{ccc}{\lambda }_1& 0& 0\\ 0& {\lambda }_2& 0\\ 0& 0& {\lambda }_3\end{array}\right).$$

$$x(t)=c_1{e}^{{\lambda }_{1}t},y(t)=c_2{e}^{{\lambda }_{2}t},z(t)=c_3{e}^{{\lambda }_{3}t},$$

$$x_{k+1}=x_k{e}^{{\lambda }_{1}h},y_{k+1}=y_k{e}^{{\lambda }_{2}h},z_{k+1}=z_k{e}^{{\lambda }_{3}h}.$$

$$x_{k+1}=\frac{\psi +\varphi {\lambda }_1(1-\theta )}{1-\varphi {\lambda }_1\theta }{x}_k,y_{k+1}=\frac{\psi +\varphi {\lambda }_2(1-\theta )}{1-\varphi {\lambda }_2\theta }{y}_k,z_{k+1}=\frac{\psi +\varphi {\lambda }_3(1-\theta )}{1-\varphi {\lambda }_3\theta }{z}_k.$$

$$\begin{array}{r}\psi +\varphi {\lambda }_1(1-\theta )=e^{{\lambda }_{1}h}(1-\varphi {\lambda }_1\theta ),\\ \psi +\varphi {\lambda }_2(1-\theta )=e^{{\lambda }_{2}h}(1-\varphi {\lambda }_2\theta ),\\ \psi +\varphi {\lambda }_3(1-\theta )=e^{{\lambda }_{3}h}(1-\varphi {\lambda }_3\theta ).\end{array}$$

$$\begin{array}{r}\varphi ({\lambda }_1-{\lambda }_2)+\varphi \theta C_1=e^{{\lambda }_{1}h}-e^{{\lambda }_{2}h},\\ \varphi ({\lambda }_2-{\lambda }_3)+\varphi \theta C_2=e^{{\lambda }_{2}h}-e^{{\lambda }_{3}h},\end{array}$$

$$C_1={\lambda }_2-{\lambda }_1+{\lambda }_1{e}^{{\lambda }_{1}h}-{\lambda }_2{e}^{{\lambda }_{2}h},C_2={\lambda }_3-{\lambda }_2+{\lambda }_2{e}^{{\lambda }_{2}h}-{\lambda }_3{e}^{{\lambda }_{3}h}.$$

$$\frac{{\lambda }_1-{\lambda }_2+\theta C_1}{{\lambda }_2-{\lambda }_3+\theta C_2}=\frac{e^{{\lambda }_{1}h}-e^{{\lambda }_{2}h}}{e^{{\lambda }_{2}h}-e^{{\lambda }_{3}h}}.$$

$$\begin{array}{rl}& \theta =\frac{T_1}{T_2},T_1={\lambda }_1(e^{{\lambda }_{2}h}-e^{{\lambda }_{3}h})+{\lambda }_2(e^{{\lambda }_{3}h}-e^{{\lambda }_{1}h})+{\lambda }_3(e^{{\lambda }_{1}h}-e^{{\lambda }_{2}h}),\\ & T_2={\lambda }_1(1-e^{{\lambda }_{1}h})(e^{{\lambda }_{2}h}-e^{{\lambda }_{3}h})+{\lambda }_2(1-e^{{\lambda }_{2}h})(e^{{\lambda }_{3}h}-e^{{\lambda }_{1}h})+{\lambda }_3(1-e^{{\lambda }_{3}h})(e^{{\lambda }_{1}h}-e^{{\lambda }_{2}h}).\end{array}$$

$$\varphi =\frac{e^{{\lambda }_{1}h}-e^{{\lambda }_{2}h}}{{\lambda }_1-{\lambda }_2+\theta C_1},\psi =e^{{\lambda }_{3}h}-\varphi {\lambda }_3(e^{{\lambda }_{3}h}\theta +1-\theta ),$$

$$J_1^{\ast }=\left(\begin{array}{ccc}0& 0& 0\\ 0& {\lambda }_2& 0\\ 0& 0& {\lambda }_3\end{array}\right).$$

$$\psi =1,\varphi =\frac{({\lambda }_1-{\lambda }_2)(e^{{\lambda }_{1}h}-1)(e^{{\lambda }_{2}h}-1)}{{\lambda }_1{\lambda }_2(e^{{\lambda }_{1}h}-e^{{\lambda }_{2}h})},\theta =\frac{{\lambda }_2(e^{{\lambda }_{1}h}-1)-{\lambda }_1(e^{{\lambda }_{2}h}-1)}{({\lambda }_1-{\lambda }_2)(e^{{\lambda }_{1}h}-1)(e^{{\lambda }_{2}h}-1)}.$$

$${\lambda }_{1,2}=\alpha \pm \beta i,{\lambda }_3=\lambda ,\beta \ne 0,\alpha ,\beta ,\lambda \in \mathbb{R},$$

$$\begin{array}{rl}\theta & =\frac{T_1}{T_2},\varphi =\frac{2e^{\alpha h}\sin (\beta h)}{2\beta +T_3\theta },\psi =e^{\lambda h}-\varphi \lambda (\theta e^{\lambda h}+1-\theta ),\\ T_1& =2\beta (e^{\alpha h}\cos (\beta h)-e^{\lambda h})+2e^{\alpha h}(\lambda -\alpha )\sin (\beta h),\\ T_2& =\alpha (1-e^{\alpha h}\cos (\beta h))(-2e^{\alpha h}\sin (\beta h))+\alpha e^{\alpha h}\sin (\beta h)(2e^{\lambda h}-2e^{\alpha h}\cos (\beta h))\\ & +\beta (1-e^{\alpha h}\cos (\beta h))(2e^{\alpha h}\cos (\beta h)-2e^{\lambda h})+\beta e^{\alpha h}\sin (\beta h)(-2e^{\alpha h}\sin (\beta h))\\ & +\lambda (1-e^{\lambda h})(2e^{\alpha h}\sin (\beta h)),\\ T_3& =-2\beta +2\alpha e^{\alpha h}\sin (\beta h)+2\beta e^{\alpha h}\cos (\beta h).\end{array}$$

$$J_1^{\ast \ast }=\left(\begin{array}{ccc}\alpha +\beta i& 0& 0\\ 0& \alpha -\beta i& 0\\ 0& 0& \lambda \end{array}\right),$$

$$J_2=\left(\begin{array}{ccc}{\lambda }_1& 0& 0\\ 0& {\lambda }_1& 0\\ 0& 0& {\lambda }_2\end{array}\right),J_3=\left(\begin{array}{ccc}{\lambda }_1& 1& 0\\ 0& {\lambda }_1& 0\\ 0& 0& {\lambda }_2\end{array}\right).$$

Full text

1 Introduction

The concepts of nonstandard finite difference schemes (NSFD) and exact finite difference schemes for differential equations were introduced by R. E. Mickens in $1980$ (see mickens1 ; mickens2 ; mickens4 ). According to this the exact finite difference schemes are those NSFD, whose solution coincides with the exact solution of the differential equations at grid points. Nonstandard finite-difference (NSFD) schemes and exact finite difference schemes have become popular in recent years (see e.g. AL ; dk1 ; Mickens3 ; Roeger2 ; Roeger3 ; Roeger4 ; Roeger5 ; Roeger6 ; Roeger7 ; Roeger8 ; Roeger9 ) mainly because some methods are more efficient on preserving certain qualitative properties in the original differential equations or systems. A good review of NSFD methods can be found in mickens4 .

There are a lot of results of EDS for both ordinay and partial differential equations such as Mickens3 ; Mickensp2 ; Roegerp1 ; Roeger9 ; Roegerp3 ; Lapinska ; Zibaei . Among them EDS for linear differential equations or system of differential equations with constant coefficients have attracted a special interest Mickens3 ; Roeger3 ; Roeger4 ; Roegerp3 .

Recently in 2008, Roeger Roeger4 constructed exact difference schemes for the system of two differential equations with constant coefficients

[TABLE]

in the form

[TABLE]

where $\theta\in\mathbb{R}$ and $\phi(h)=h+\mathcal{O}(h^{2}),\quad h\to 0$.

The main idea of the construction is that instead of the system $\textbf{x}^{\prime}=A\textbf{x}$ with a general $2\times 2$ matrix $A$ the author considered the system $\textbf{u}^{\prime}=J\textbf{u}$, where $J$ is a $2\times 2$ Jordan canonical form. For each case of the Jordan forms $J$, the parameters $\theta$ and $\phi$ are determined so that (2) is an exact difference scheme for the system $\textbf{u}^{\prime}=J\textbf{u}$. Finally, due to the fact that $A$ is similar to $J$ the exact finite difference schemes for $\textbf{u}^{\prime}=J\textbf{u}$ also are exact finite difference schemes for $\textbf{x}^{\prime}=A\textbf{x}$. The obtained results show that $\theta$ and $\phi$ depend only on the eigenvalues of $A$.

Before, Mickens mickens1 ; mickens2 ; mickens4 constructed EDS for system of two linear differential equations with constant coefficients

[TABLE]

in the form

[TABLE]

where

[TABLE]

and $\lambda_{1}$, $\lambda_{2}$ are the roots of the characteristic equation

[TABLE]

EDS proposed by Mickens and Roeger contain two parameters. However, different from the Roeger’s scheme Mickens’ scheme contains additional parameter $\psi=1+\mathcal{O}(h^{2})$ in discretization of the first derivative while the right hand side is locally discretized.

It is possible to say that system of linear differential equations is simple, so if there are available eigenvalues and Jordan structure of the coefficient matrix then we have an explicit formula for solution. But it should be emphasized that for EDS we have to compute parameters once at beginning, after that compute solution recurrently, therefore the computational cost is cheaper than computing by formula of solution (including the exponential and trigonometric functions) with the guarantee of accuracy of the result.

Therefore, in this paper we extend the Roeger’s idea and Mickens’s idea for constructing exact finite difference schemes for the system of three equations with constant coefficients

[TABLE]

Differently from the Roeger’s approximation of the first derivative with the use of one parameter $\phi(h)$ in denominator, here we add a parameter $\psi$ in numerator so that the number of parameters is equal to the number of differential equations in the system. Namely, instead of the difference scheme (2) we use the difference scheme of the form

[TABLE]

where $\psi(h)=1+\mathcal{O}(h^{2})$.

Notice that the difference scheme (2) contains only two parameters $\phi$ and $\theta$, therefore it could be exact for two-dimensional system of equations. In some special cases it may be exact for three-dimensional system of equations but in general case of three-dimensional system two parameters difference scheme (2) will not be exact. It is the reason why we introduce an additional parameter $\psi$ into the scheme (4).

In general, a system of $n$ linear differential equations with constant coefficients has a fundamental system of solutions including $n$ functions, therefore an EDS must contain at least $n$ parameters. The Roeger’s exact scheme and Mickens’ scheme contain two parameters, therefore, they cannot ensure exactness for system of three linear equations. The addition of the parameter $\psi$ into the Roeger’s scheme is a simple extension of ours.

In general case, it is possible construct EDS for system of $n$ equations based on Runge-Kutta methods. Namely, applying a s-stage Runge-Kutta method Ascher ; Hairer1 ; Hairer2 with coefficient matrix A_{RK}=\big{(}a^{*}_{ij}\big{)}_{s\times s} and coefficients b_{RK}=\big{(}b^{*}_{1},b^{*}_{2},\ldots,b^{*}_{s}\big{)}^{T} and c=\big{(}c^{*}_{1},c^{*}_{2},\ldots,c^{*}_{s}\big{)}^{T} to the system $\textbf{x}^{\prime}=A\textbf{x}$ we obtain the scheme

[TABLE]

where the coefficients $\alpha_{m}=\alpha_{m}(a_{ij}^{*},b_{i}^{*}),m=\overline{2,s}$ depend on $A_{RK}$ and $b_{RK}$. In the scheme (5) replacing $h$ by the function $\phi(h)=h+\mathcal{O}(h^{2})$ and adding the parameter $\psi(h)=1+\mathcal{O}(h^{2})$ we obtain NSFD scheme for the system $\textbf{x}^{\prime}=A\textbf{x}$

[TABLE]

This scheme (6) contains $s+1$ parameters. Therefore, it is possible construct EDS for the system $\textbf{x}^{\prime}=A\textbf{x}$ with the dimension $n\leq s+1$ from the scheme (6). Analogously, consider the following implicit difference schemes

[TABLE]

[TABLE]

Depending on the dimension of the systems of equations we can choose the suitable number of parameters for the difference schemes to be exact. For example, in the case $n=2$ Roeger considers the scheme (8) with $\psi=1,\alpha_{m}=0,m=\overline{2,s}$. In the case $n=3$ under consideration we choose the scheme (4) which is a particular case of the scheme (8) with $\alpha_{m}=0,m=\overline{2,s}$. It may be considered as a natural extension of the results of Roeger and Mickens. Besides, we choose the explicit scheme (6) with $\alpha_{m}=0,m=\overline{3,s}$, that is, the scheme of the form

[TABLE]

In the case of $n\geq 3$ dimensions we can do in a similar way. In general, it is possible to construct EDSs based on the schemes of the form (6), (7), (8) combined with the use of Jordan forms of matrices.

In this work, we show that any three-dimensional linear system $\textbf{x}^{\prime}(t)=A\textbf{x}(t)$, \textbf{x}(t)=\big{(}x(t),y(t),z(t)\big{)}^{T},A\in M_{3\times 3}(\mathbb{R}), has an exact finite-difference method in the forms (4) and (9), where $\psi$, $\phi$ and $\theta$ can be found explicitly in terms of the step-size $h$ and the eigenvalues $\lambda_{1,2,3}$ of the coefficient matrix $A$. In Section $2$, we prove that if the parameters $\psi$, $\phi$ and $\theta$ are determined so that (4)/(9) is the exact difference scheme for $\textbf{u}^{\prime}=J\textbf{u}$, where $J$ is $3\times 3$ Jordan form matrix then (4)/(9) also will be exact for $\textbf{x}^{\prime}=A\textbf{x}$ if $A$ is similar to $J$. Based on this fact, in Section $3$ and Section $4$ we construct implicit and explicit exact difference schemes for the system $\textbf{x}^{\prime}=J\textbf{x}$. Next, in Section 5 we make a perturbation analysis for estimating the accuracy of EDS in the case of appearing of rounding errors due to the approximate computation of the parameters. In Section $6$ we report some numerical examples for stiff problems and problems with special properties on long time interval for demonstrating the efficiency and exactness of EDS in comparison with high-order numerical methods. Some concluding remarks will be given in the last section.

2 Why consider the system with Jordan form matrix?

From Linear Algebra it is well known that any $n\times n$ matrix $A$ is similar to a Jordan form matrix $J$. In the case $n=3$ it is easy to list all Jordan form matrices $J$ as stated in the following theorem (see e.g. (SW, , Chapter 1), (KW, , Chapter 6)).

Theorem 2.1

Let $A$ be any $3\times 3$ matrix. Then, $A$ is similar to one of the following Jordan form matrices $J$ depending on the set of its eigenvalues and the dimension of eigenspaces associated with the eigenvalues. Here, $\sigma(A)$ is the set of eigenvalues, $\chi_{A}(t)$ and $m_{A}(t)$ are characteristic and minimal polynomials of $A$, respectively.

Now, we consider the three-dimensional system of differential equations with constants coefficients (3) and NSDF schemes of the form (4). Denote by $\mathcal{J}$ the set of all Jordan form $3\times 3$ matrices.

Theorem 2.2

Suppose that the difference scheme

[TABLE]

is exact for the system $\textbf{u}^{\prime}=J\textbf{u},J\in\mathcal{J}$. Then the difference scheme (4) with the same parameters $\psi$, $\phi$, $\theta$ is exact for the system $\textbf{x}^{\prime}=A\textbf{x}$, $A\in M_{3\times 3}(\mathbb{R})$ if $A$ is similar to $J.$

Proof

Suppose that $A$ is similar to the Jordan form matrix $J\in\mathcal{J}$. Then there exists a invertable matrix $P$ such that $P^{-1}AP=J$. Making the transformations $\textbf{u}=P^{-1}\textbf{x}$ and $\textbf{u}_{k}=P^{-1}\textbf{x}_{k}$ we convert the system $\textbf{u}^{\prime}=J\textbf{u}$ to $\textbf{x}^{\prime}=A\textbf{x}$ and the difference scheme (10) to (4), respectively. Since (10) is exact for $\textbf{u}^{\prime}=J\textbf{u}$, the difference scheme (4) is exact for $\textbf{x}^{\prime}=A\textbf{x}$.

The following theorem is a generalization of Theorem 2.3 and is proved in a completely similar way.

Theorem 2.3

Suppose that the difference scheme (6)/(7)/(8) is exact for the system $\textbf{u}^{\prime}=J\textbf{u},J\in\mathcal{J}$. Then the difference scheme (6)/(7)/(8) is exact for the system $\textbf{x}^{\prime}=A\textbf{x}$, $A\in M_{n\times n}(\mathbb{R})$ if $A$ is similar to $J$.

From the above theorems we see that it suffices to construct exact difference schemes for three-dimensional systems with Jordan form matrices.

3 Implicit exact difference schemes (IEDS) for $\textbf{x}^{\prime}=J\textbf{x}$

In this section we construct implicit EDS for the system $\textbf{x}^{\prime}=A\textbf{x}$ in the form (4). To avoid the introduction of new variables we shall use x instead of u in the system of equations with Jordan form matrix. We construct exact difference schemes for this system and it suffices to replace $J$ by $A$ to obtain exact difference schemes for the system with the matrix $A$.

3.1 The case when $A$ has 3 distinct eigenvalues

In this case $A$ is similar to the Jordan form matrix

[TABLE]

First we consider the subcase $\lambda_{1}\lambda_{2}\lambda_{3}\neq 0$. Then the linear system $\textbf{x}^{\prime}=J_{1}\textbf{x}$ has the exact solution

[TABLE]

or equivalently

[TABLE]

Applying the difference scheme (4) for the system $\textbf{x}^{\prime}=J\textbf{x}$ we obtain

[TABLE]

Identifying (13) and (12) we come to the system for finding the parameters $\psi,\phi,\theta$:

[TABLE]

After some elementary transformations we can eliminate $\psi$ and obtain the system

[TABLE]

where for brevity we set

[TABLE]

Dividing the first equation by the second one in (15) we obtain the equation containing only one unknown $\theta$

[TABLE]

The solution of this equation is

[TABLE]

After $\theta$ is found, returning to (15) and (14) we obtain $\phi$ and $\psi$, namely

[TABLE]

where $C_{1}$ is given by (16). They are the parameters to be determined.

Theorem 3.1

The linear system (3) with the matrix of coefficients (11) has an exact difference scheme of the form (4), where the parameters $\psi,\phi,\theta$ are determined by (17) and (18).

Now we consider the special subcase $\lambda_{1}\lambda_{2}\lambda_{3}=0$. Without generality we suppose $\lambda_{1}=0$. Then, $A$ is similar to the Jordan form matrix

[TABLE]

In this subcase, analogously as in the previous subcase, we find the parameters in the exact difference scheme

[TABLE]

This result coincides with that of Roeger for the system of two linear equations with the coefficient matrix having two distinct eigenvalues $\lambda_{1}\neq\lambda_{2}$ and both are nonzeros.

Theorem 3.2

The linear system (3) with the coefficient matrix (19) has an exact difference scheme of the form (4), where the parameters $\psi,\phi,\theta$ are given by (20).

Remark that, when the matrix $A$ has a pair of complex conjugate eigenvalues

[TABLE]

following Theorem 3.1 we obtain the parameters $\psi,\phi,\theta$

[TABLE]

Corollary 1

The linear system (3) with the coefficient matrix

[TABLE]

has an exact difference scheme of the form (4), where the parameters $\psi,\phi,\theta$ are determined by (21).

3.2 The case $A$ has eigenvalues $\lambda_{1}=\lambda_{2}\neq\lambda_{3}$

In this case the matrix $A$ is similar to one of the following Jordan form matrices

[TABLE]

**(i). When $A$ is similar to $J_{2}$

**

In this subcase the linear system (3) with the coefficient matrix $J_{2}$ has exact solution

[TABLE]

or equivalently

[TABLE]

Applying the difference scheme (4) for $\textbf{x}^{\prime}=J_{2}\textbf{x}$ we obtain

[TABLE]

Identifying (24) and (23) we come to the system of conditions for determining $\psi,\phi,\theta$

[TABLE]

The above system (25) is of two equations but contains three parameters. Therefore, it has infinite number of solutions. It is easy to find a dependence of $\psi,\phi$ on $\theta$ as follows

[TABLE]

Theorem 3.3

The linear system (3) with the coefficient matrix

[TABLE]

has an exact difference scheme of the form (4), where the parameters $\psi,\phi,\theta$ satisfy the relations (26).

**(ii). When $A$ is similar to $J_{3}$

**In this case the system (3) with the coefficient matrix $J_{3}$ has an exact solution

[TABLE]

or equivalently

[TABLE]

Applying the difference scheme (4) to the system $\textbf{x}^{\prime}=J_{3}\textbf{x}$ we obtain

[TABLE]

Suppose the above system is not degenerate. It occurs if $(1-\phi\lambda_{1}\theta)(1-\phi\lambda_{2}\theta)\neq 0.$ Then the system has a unique solution

[TABLE]

Identifying (30) and (28) we come to the system of conditions for determining the parameters $\psi$, $\phi$, $\theta$:

[TABLE]

If $\lambda_{1}=0$ then from the above relations we get

[TABLE]

Otherwise, setting $\phi\theta=T$ from the third equation of (31) we obtain a quadratic equation for $T$

[TABLE]

Subtracting $1$ from the first and second equations of (31) we obtain a system of two equations for $\psi-1$ and $\phi$, where $T$ is considered as a parameter.

[TABLE]

Due to $\lambda_{1}\neq\lambda_{2}$ the system has a solution

[TABLE]

Substituting (34) into (33) we obtain a quadratic equation for $T$

[TABLE]

The equation (35) has the discriminant $\Delta=(e^{\lambda_{1}h}-e^{\lambda_{2}h})^{2}>0$, hence it has two distinct roots $T_{1},T_{2}$ given by

[TABLE]

[TABLE]

After finding $T=\phi\theta$, we can calculate $\psi$ and $\phi$ by the formulas (34), and $\theta=T/\phi$.

The following proposition implies that only the value $T=T_{1}$ is consistent with the assumption of non-degeneration of the system (29)

Proposition 1

For $T_{1}$ and $T_{2}$ determined by (35) and (36) we have

[TABLE]

Proof

The proposition is easily proved by the use of L’Hospital’s rule for the indeterminate form of the type $0/0$.

Thus, in the case of the Jordan form matrix $J_{3}$ the parameters of the difference scheme are determined as follows

[TABLE]

Theorem 3.4

The linear system (3) with the coefficient matrix

[TABLE]

has an exact difference scheme of the form (4), where the parameters $\psi,\phi,\theta$ are given by (38), where $T_{1}$ is determined by (36) if $\lambda_{1}\neq 0$ and by (32) if $\lambda_{1}=0$.

Remark 1

The system of conditions (31) is obtained from (25) by adding the third equation. Therefore, the parameters satisfying (31) also satisfies (25), i.e., the difference scheme with these parameters is exact for the linear system (3) having the matrix $A$ similar to $J_{2}$. Thus, in the case if $A$ has the set of eigenvalues \sigma(A)=\big{\{}\lambda_{1},\lambda_{1},\lambda_{2}\big{\}} then the exact difference scheme may be determined by Theorem 3.4 not depending on the similar Jordan form matrices.

3.3 The case $A$ has eigenvalues $\lambda_{1}=\lambda_{2}=\lambda_{3}=\lambda$

In this case $A$ is similar to one of the following three Jordan form matrices

[TABLE]

**(i). When $A$ is similar to $J_{4}$

**In this subcase the linear system with the coefficient matrix $J_{4}$ has the exact solution

[TABLE]

or equivalently

[TABLE]

Applying the difference scheme (4) to the system $\textbf{x}^{\prime}=J_{4}\textbf{x}$ we obtain

[TABLE]

Identifying (41) and (40) we obtain

[TABLE]

Remark 2

The equation (42) contains three unknowns. So, it has infinitely number of solutions, namely, the set of its solutions is a two-dimensional linear space. A simple solution of it is $\theta=0,\psi=1$ and $\phi=\dfrac{e^{\lambda h}-1}{\lambda}$ ($\lambda\neq 0$).

Theorem 3.5

The linear system (3) with the coefficient matrix

[TABLE]

has an exact difference scheme of the form(4), where the parameters $\psi,\phi,\theta$ satisfy (42).

**(ii). When $A$ is similar to $J_{5}$

**

In this subcase the linear system with the coefficient matrix $J_{5}$ has the exact solution

[TABLE]

or equivalently

[TABLE]

Applying the difference scheme (4) to the system $\textbf{x}^{\prime}=J_{5}\textbf{x}$ we obtain

[TABLE]

Identifying (44) and (43) we obtain the system of conditions for $\psi,\phi,\theta$

[TABLE]

Remark 3

The system (45) is of two equations with three unknowns, so, it has infinitely many solutions, namely, its set of solutions is an one-dimensional linear space. A simple solution of it is $\theta=0$, $\phi(h)=he^{\lambda h}$, $\psi=(1-\lambda h)e^{\lambda h}.$

Theorem 3.6

The linear system (3) with the coefficient matrix

[TABLE]

has an exact difference scheme of the form(4), where the parameters $\psi,\phi,\theta$ satisfy (45).

**(iii). When $A$ is similar to $J_{6}$

**In this subcase (4) with the coefficient matrix $J_{6}$ has the exact solution

[TABLE]

or equivalently

[TABLE]

Applying the difference scheme (4) to the system $\textbf{x}^{\prime}=J_{6}\textbf{x}$ we obtain

[TABLE]

Under the condition $\phi\lambda\theta\neq 1$ we have

[TABLE]

Identifying (48) and (47) we obtain the system of conditions for $\psi,\phi,\theta$

[TABLE]

By some elementary calculations we have found

[TABLE]

Clearly, $\phi\lambda\theta\neq 1$ if $h\neq-2/\lambda$. Notice that when $\lambda_{1}=\lambda_{2}=\lambda_{3}=0$ we obtain

[TABLE]

Theorem 3.7

The linear system (3) with the coefficient matrix

[TABLE]

has an exact difference scheme of the form (4), where the parameters $\psi,\phi,\theta$ are determined by (50).

Remark 4

The equation (42) for determining $\psi,\phi,\theta$ in the case of $J_{4}$ is the first equation in (49), and the system (45) in the case of $J_{5}$ is the first two equations in (49) in the case of $J_{6}$. Therefore, the parameter $\psi,\phi,\theta$ found from (49) are applicable for all cases of $J_{4}$, $J_{5}$ and $J_{6}$. It means that in the case if $A$ has the set of eigenvalues \sigma(A)=\big{\{}\lambda,\lambda,\lambda\big{\}} then the exact difference scheme can be determined by Theorem 3.7.

3.4 Summary of results

The obtained above results of exact difference schemes for the three-dimensional linear system (3) with constant coefficient matrix $A$ can be summarized in Table 2 below

Concerning the parameters $\psi$, $\phi$ and $\theta$ as functions of the stepsize of discretization it is easy to verify

Lemma 1

*For the functions $\psi,\phi$ and $\theta$ determined by theorems in Table 3 we have

$(i).\lim_{h\to 0^{+}}\psi(\lambda_{1},\lambda_{2},\lambda_{3})=1,$

$(ii).\lim_{h\to 0^{+}}\phi(\lambda_{1},\lambda_{2},\lambda_{3})/h=1,$

$(iii).\lim_{h\to 0^{+}}\theta(\lambda_{1},\lambda_{2},\lambda_{3})=\dfrac{1}{2}$.*

4 Explicit exact difference schemes (EEDS) for $\textbf{x}^{\prime}=A\textbf{x}$

In this section we construct explicit EDS for the system $\textbf{x}^{\prime}=A\textbf{x}$ in the form (9). From the previous section and Roeger4 we conclude that there exist an EDS for systems with the matrix having a same set of eigenvalues, i.e. EDS does not depend on the Jordan structure of the matrix $A$. Specifically, for all matrices with a same set of eigenvalues, if a difference scheme is exact for a system with matrix having higher minimal polynomial then it is also the exact scheme for a system with matrix having lower minimal polynomials. Therefore, when constructing EDS for systems with matrices having a same set of eigenvalues it suffices to consider the case of matrix having highest minimal polynomial. It is why in this section instead of 6 cases of the matrix $J$ as in Section 3 we consider only 3 cases of the matrix $J$ as follows.

4.1 The case when $A$ has 3 distinct eigenvalues

In this case matrix $A$ is similar to $J_{1}$. Applying the difference scheme (9) for the system $\textbf{x}^{\prime}=J_{1}\textbf{x}$ we obtain

[TABLE]

Identifying (53) and (12) we come to the system for finding the parameters $\phi,\psi,\theta$:

[TABLE]

Suppose that $\lambda_{1}+\lambda_{2}\neq 0$. Then by consecutive eliminations it is easy to obtain the solution of the system (54):

[TABLE]

4.2 The case $A$ has eigenvalues $\lambda_{1}=\lambda_{2}\neq\lambda_{3}$

In this case it suffices to consider the case when matrix $A$ similar to $J_{3}$. Analoguously as above we obtain a system of conditions for determining the parameters $\phi,\psi,\theta$:

[TABLE]

If $\lambda_{1}=0$ then from the system (56) we obtain

[TABLE]

Otherwise, if $\lambda_{1}\neq 0$ it is easy to get the solution of the system (56)

[TABLE]

4.3 The case $A$ has eigenvalues $\lambda_{1}=\lambda_{2}=\lambda_{3}=\lambda$

In this case it suffices to consider the case when the matrix $A$ is similar to $J_{6}$. Then the system of conditions for the parameters $\phi,\psi,\theta$ is

[TABLE]

If $\lambda=0$ from (59) we obtain $\psi=1,\phi=h,\theta=\dfrac{1}{2}$. Otherwise, if $\lambda\neq 0$ the solution of (59) is

[TABLE]

Remark 5

From Section 3 and Section 4 we see that the system of conditions for determining the parameters $\phi,\psi,\theta$ for implicit EDS is much more complicated than for explicit EDS. Specifically, the system of conditions for implicit EDS contains rational expressions while the system of conditions for explicit EDS contains polynomial expressions.

5 Perturbation analysis

Since the parameters of EDS contain exponential and trigonometric functions, in the process of computation rounding errors arise. Suppose that instead of the exact parameters $\psi,\phi,\theta$ we obtain only their approximate values $\hat{\psi},\hat{\phi},\hat{\theta}$. Notice that the explicit EDS and implicit EDS for the system $\textbf{x}^{\prime}=A\textbf{x}$ can be written in the form $\textbf{x}_{k+1}=Q(\psi,\phi,\theta)\textbf{x}_{k}:=Q\textbf{x}_{k}$. Due to the fact that the iterative parameters are computed approximately instead of $Q$ we only have $\hat{Q}:=(\hat{\psi},\hat{\phi},\hat{\theta}).$ Suppose

[TABLE]

Therefore, we obtain only the approximation $\hat{\textbf{x}_{k}}$ but not $\textbf{x}_{k}$ . It is easy to obtain the difference between $\hat{\textbf{x}_{k}}$ and $\textbf{x}_{k}$

[TABLE]

From here it follows

[TABLE]

We see that the error of computed solution mainly depends on the number of iterations and the rounding errors in computation of the parameters. When the number of iterations is large the error of EDS may be large. However, this error slightly depends on $h$, therefore we can overcome this phenomena as follows: instead of computing $\textbf{x}_{k+1}$ through $\textbf{x}_{k}$ by the formula $\textbf{x}_{k+1}=Q(\psi(h),\phi(h),\theta(h))\textbf{x}_{k}$ with grid size $h=T/N$ we can compute $\textbf{x}_{k}$ through $\textbf{x}_{0}$ with grid size $h^{*}=t_{k}$, i.e. use the formula $\textbf{x}_{k}=Q(\psi(h^{*}),\phi(h^{*}),\theta(h^{*})\textbf{x}_{0}$. Equivalently, instead of solving the problem on interval $[0,t_{k}]$ with $h=T/N$ via a large number of steps we can compute $x_{k}$ via $x_{0}$ with the step $h^{*}=t_{k}$. Thus, after only one iteration we obtain $\textbf{x}_{k}$. This is the advantage of EDS compared with high-order numerical methods because for ensuring the accuracy these methods must use small grid sizes. From the numerical examples in the next section it will be seen that even in the presence of rounding errors EDS are more efficient than high-order numerical methods.

6 Numerical simulations

In this Section we perform some numerical simulations for confirming the validity of theoretical results obtained in the previous sections. The numerical simulations for a stiff problem and problems with specific properties demonstrate the advantage of EDS over high-order numerical methods.

Example 1

Consider the system (3) with the coefficient matrix

[TABLE]

The set of eigenvalues of $A$ is $\sigma(A)=\{-1,\pm i\}$. For the initial conditions $x(0)=0,y(0)=-50,z(0)=50$, the system has the exact solution

[TABLE]

The implicit exact difference scheme for the system are determined from Theorem 3.2. The exact solution of the system and the the solution of the exact difference schemes are depicted in Figures 1.

Example 2

Consider the system (3) with the coefficient matrix

[TABLE]

The set of eigenvalues of $A$ is $\sigma(A)=\{0,-1\}$. For the initial conditions $x(0)=0,y(0)=-40,z(0)=50$, the system has the exact solution

[TABLE]

The explicit exact difference schemes for the system are determined by Subsection 4.2. The exact solution of the system and the the solution of the exact difference schemes are depicted in Figures 2.

From the two above numerical examples we see that the solution of the constructed EDS almost coincide with the exact solution of the system of differential equations at grid points due to the insignificant rounding errors in computing the parameters of EDS. In the ideal case when the rounding errors are absent the solution of EDS must coincide with the solution of the system of differential equations for any grid size $h$.

Example 3

Consider the system $\textbf{x}^{\prime}=A\textbf{x}$, $t\in[0,T]$ with the coefficient matrix

[TABLE]

where $\lambda>0$. The set of eigenvalues of $A$ is $\sigma(A)=\{\pm i,\lambda\}$. For the initial conditions $x(0)=1,y(0)=0,z(0)=1$ and $\lambda=1$, the system has the exact solution $x(t)=\cos(t)$, $y(t)=\sin(t)$, $z(0)=e^{\lambda t}$. The components $x(t)$, $y(t)$ are periodic and $x^{2}(t)+y^{2}(t)=1$ for all $t\in[0,T]$. It is easy to proved that (see (Ascher, , Problem 1, Section 3.9)):

The solution $(x_{k},y_{k})$ obtained by the explicit Euler method spirals out. 2. 2.

The solution $(x_{k},y_{k})$ obtained by the implicit Euler method spirals in. 3. 3.

The solution $(x_{k},y_{k})$ obtained by the trapezoidal method forms an approximate circle as desired.

In general, many numerical methods of higher order of accuracy such as Runge-Kutta or Taylor methods cannot preserve invariant properties of the differential problems although their global errors tend to zero as $h\rightarrow 0$. It means that the choice of small grid size only ensures the accuracy of the methods but not ensure the invariant properties of the problems. For the methods of higher order of accuracy including one step methods and multistep methods when solving the problems numerically it is necessary to take grid size small because all convergence theorems are stated for $h\rightarrow 0$. Therefore, when the final time $T>>1$ it is impossible to choose grid size very small because the number of steps of computation becomes extremely large and this may cause difficulty with computer memory and computation time. Moreover, when $h$ is very small the accuracy may decrease due to rounding errors.

Now we compare the accuracy of EDS with some typical higher order methods such as Runge-Kutta and Taylor methods Ascher ; Hairer1 ; Hairer2 .

In this example $\lambda>0$, therefore the problem is unstable and it is not necessary to use implicit A-stable or L-stable Runge-Kutta methods. Moreover, since the system is linear, implicit schemes are easily reduced to explicit ones, therefore, we shall use the classical four-stage Runge-Kutta method. Besides, we consider the Taylor method of $5$ order of accuracy

Suppose, we need to find the approximate value for the exact solution at the time $t=T$. The comparison of errors of the methods are given in Table 3, where $error=|x_{N}-x(t_{N})|+|y_{N}-y(t_{N})|+|z_{N}-z(t_{N})|$ is used as a measure of accuracy of methods, $x_{N},y_{N},z_{N}$ are computed solution by methods with the grid size $h=T/N$, $x(t_{N}),y(t_{N}),z(t_{N})$ are values of the exact solution at $t_{N}=Nh=T$. For methods of higher order of accuracy we take small grid size for guaranteeing accuracy and convergence, but for EDS it is not needed. Hence, we use the grid size $h=T$ for avoiding the decrease of accuracy after a large number of iterations due to rounding errors.

From Table 3, where IEDS and EEDS stand for Implicit and Explicit Exact Difference Schemes, respectively, we see that the methods of higher order of accuracy have small errors when $T$ and grid size $h$ are small. But when $T$ is large, despite small grid size $h$, the accuracy of these methods decreases due to the accumulation of rounding errors after a large number of iterations. This occurs because for computing the approximate value of the exact solution at the final time $T$ it is needed to compute consecutively the approximate values of the exact solution at every time before $T$. For example, for the fourth order Runge-Kutta method, theoretically, in order to obtain the approximate value of the solution at $T=1$ with the accuracy $10^{-16}$ we have to choose the grid size $h=10^{-4}$, and perform $10.000$ iterations. Nevertheless, in practice, the actual accuracy reached is only $10^{-14}$ due to the decrease of accuracy as the result of accumulation of rounding errors. Similar situation also occurs with the Taylor methods and other methods of higher order accuracy. Meanwhile, for EDS, if taking grid size $h=T=1$ then after exactly one step we obtain the approximate solution at the time $t_{N}=T$ with the accuracy $10^{-16}$. This completely agrees with the analysis in Section 5.

Besides, from Table 3 it is easily seen that for other values of $T\geq 10$ the accuracy of EDS depends slightly on grid sizes $h$ and it is best if $h=T$. In nature, it depends on the rounding errors of computation of the parameters of the schemes and the number of iterations. Meanwhile, from the table we also see that for large time interval the higher order methods are inefficient, even are impossible.

From the above example we can conclude that for the problems on large time intervals the exact difference schemes (implicit or explicit) are more efficient than higher order methods.

Example 4

(Stiff problem)

Consider the system $\textbf{x}^{\prime}=A\textbf{x}$, $t\in[0,T]$ with the coefficient matrix

[TABLE]

The set of eigenvalues of $A$ is $\sigma(A)=\{-1,-2,-100\}$. For the initial conditions $x(0)=1,y(0)=1,z(0)=1$, the system has the exact solution $x(t)=e^{-t}$, $y(t)=e^{-2t}$, $z(t)=e^{-2017t}$. Obviously, all the components of the solution monotonically tend to zero with the exponential rate. As is well known, for efficiently solving stiff problems it is necessary to use methods having stability properties such as A-stability or L-stability (see Ascher ; Hairer1 ; Hairer2 . Some implicit Runge-Kutta methods possess these stability propoerties, while explicit Runge-Kutta methods cannot have L-stability since they have bounded stability regions. In general for stiff problems explicit methods are inefficient.

In this example we compare EDS with the classical four-stage Runge-Kutta method, five-order Taylor method and five-order Radau IIA method (Hairer1, , Table II.7.7). The results of computation are reported in Table 4, where

[TABLE]

is a measure of accuracy of methods.

Analogously as Example 3, from Table 4 we see that EDS are much more efficient than RK4 and Taylor methods. Although Radau IIA gives errors better than RK4 and Taylor but it is implicit method, therefore, it requires more computational cost due to the multiplication of matrix by vector for determining stages of the method. In general, the computational cost of higher order methods is much more than one of EDS.

Example 5

(Nonstandard finite difference scheme of combined type for quasi-nonlinear system of differential equations)

Consider the quasi-nonlinear system of equations (see Problem 22.9 Agarwal1 )

[TABLE]

where g\in C\Big{[}[t_{0},\infty]\times\mathbb{R}^{n},\mathbb{R}^{n}\Big{]} and $||g(t,v)||\leq\lambda(t)||v||$, where $\lambda(t)$ is a nonnegative continuous function in $[x_{0},\infty]$. Additionally, suppose that the function $\lambda(t)$ satisfies the condition $\int^{\infty}\lambda(t)dt<\infty$, $\lambda(t)\to 0$ as $t\to\infty$ and the matrix $A$ satisfies $\lambda<0$ for any $\lambda\in\sigma(A)$. Then it is easy to prove that Agarwal1 every solution of (61) is bounded and the trivial solution is asymptotically stable.

Our objective is to construct difference scheme preserving the properties of (61) for any grid size $h>0$. It should be emphasized that standard finite difference schemes cannot preserve properties of differential equations for any $h>0$ mickens1 ; mickens2 ; mickens4 . In this example the properties of the problem are decided by the matrix $A$, therefore, in the simplest way we propose NSFD scheme for (61) in the form

[TABLE]

where $U(A,v_{k},h)$ is determined so that the scheme $\dfrac{v_{k+1}-v_{k}}{\phi}=U(A,v_{k},h)$ is exact for the system $v^{\prime}=Av$. Then, by Theorem 5.3.1 in Agarwal2 the properties of the problem are preserved for any $h>0$. Of course, for ensuring the accuracy of the scheme it is needed to choose $h<<1$. The scheme (63) has only first order of accuracy, nevertheless it prompts us of a way for constructing NSFD schemes of higher order of accuracy for (61) in the form

[TABLE]

where $U(A,v_{k},h),V(f,v_{k},h)$ is determined so that the scheme $\dfrac{v_{k+1}-v_{k}}{\phi}=U(A,v_{k},h)$ and the scheme $\dfrac{v_{k+1}-v_{k}}{\phi}=V(f,v_{k},h)$ consecutively are EDS for $v^{\prime}=Av$ and a scheme of higher order of accuracy for $v^{\prime}=f(t,v)$. In near future we will develop this idea for constructing NSFD scheme of higher order of accuracy preserving the properties of the general quasi-linear system of differential equations.

7 Conclusion

In this paper, based on the technique of Mickens and Roeger of exact difference schemes, we have constructed implicit and explicit exact difference schemes (EDS) for system of three linear differential equations with constant coefficients $\textbf{x}^{\prime}=A\textbf{x}$. We have done the perturbation analysis for showing the advantage of EDS over higher order methods when solving problems on large time intervals. Numerical experiments for several problems, especially, a stiff problem and a periodic problem on large time intervals confirm the advantages of constructed EDS over higher order accuracy methods. In the future we will extend the obtained results to general linear system of $n$ equations with constant coefficients and for constructing NSFD scheme of higher order of accuracy preserving the properties of the general quasi-linear system of differential equations $x^{\prime}=Ax+f(x)$.

Acknowledgments

We would like to thank the reviewers for their helpful comments and suggestions for improving the quality of the paper.

This work is supported by Vietnam National Foundation for Science and Technology Development (NAFOSTED) under the grant number 102.01-2014.20.