A new third order convergent numerical solver for continuous dynamical systems

Sania Qureshi; Abdullahi Yusuf

doi:10.1016/j.jksus.2019.11.035

View/Download PDF

Buy Reprints

PDF

Translate this page into:

Original article

32 (

); 1409-1416

doi:

10.1016/j.jksus.2019.11.035

A new third order convergent numerical solver for continuous dynamical systems

Sania Qureshi^{^a,⁎}, Abdullahi Yusuf^{^b}

a Department of Basic Sciences and Related Studies, Mehran University of Engineering and Technology, Jamshoro 76062, Pakistan

b Federal University Dutse, Science Faculty, Department of Mathematics, 7156 Jigawa, Nigeria

⁎Corresponding author. sania.qureshi@faculty.muet.edu.pk (Sania Qureshi),

Received: 2018-11-12, Accepted: 2019-11-24, Published: 2019-03-01

Disclaimer:
This article was originally published by Elsevier and was migrated to Scientific Scholar after the change of Publisher.

Peer review under responsibility of King Saud University.

Abstract

A new numerical solver with third order convergence has been proposed.
The solver is consistent, stable, and possesses smaller error bounds.
Error analysis of the solver has been carried in detail.
The performance of the solver is better than many existing solvers having similar characteristics.

Abstract

Continuous dynamical systems are attempted to solve via proposed third order convergent numerical solver. The solver is shown to possess third order convergence on the basis of containing $O (h^{4})$ term in the leading coefficient of the local truncation error which has further been employed for the local and global error bounds whereas its increment function is found to be Lipchitz continuous. Consistency of the solver has extensively been discussed using Taylor’s series expansion. The necessary and sufficient conditions for the solver to be stable have been proved accompanying the linear stability function obtained via Dahlquist test problem. In order to illustrate the performance of the solver, few scalar and vector-valued dynamical systems of both linear and nonlinear nature have numerically been solved in comparison to two well-known numerical solvers having same order of convergence as that of the proposed solver. Whereupon, the proposed solver is found to contain the smallest errors. The analysis of comparison is based upon three kinds of errors; namely, maximum absolute global relative error, absolute relative error determined at the final nodal point along the integration interval and the $ℓ^{2}$ -error norm. MATLAB software having version 9.3.0.713579 (R2017b), using Intel(R) Core(TM) i3-4500U processor running on 1.70 GHz has been used.

Keywords

Differential equations

Local error

Global error

Lipchitz continuity

Stability

Show Related Articles from PubMed

1 Introduction

Mathematical models containing ordinary differential equations (ODEs) play a vital role in various fields of study. Areas including engineering, applied mathematics, bio-medical, statistical mechanics, quantum mechanics, economics, image processing, control theory, and many others are not considered to be complete without inclusion of subjects on differential equations and their applications in respective fields (Amann, 2011; Butcher, 2016; Hirsch et al., 2012; Jordan and Peter, 1999; Mermoud and Lausanne, 2015; Shampine, 2018; Stickler and Ewald, 2014; Teschl, 2000 ). Unfortunately, large quantity of such physical and natural models have no closed form solutions due to one or other type of non-linearity associated with the models.

Due to the universal acceptance of the importance of such mathematical models, many researchers have devised new solvers to get the desired solutions of the models under their study whereas others have substantially improved the performance of many existing solvers with regard to the convergence order, computational complexity, stability, consistency and error bounds.

In order to solve the scalar and continuous dynamical systems, the authors in (Ma and Simos, 2017; Rabiei and Fudziah, 2012; Rabiei et al., 2012a; Rabiei et al., 2012b ) have improved the efficiency for the well-known family of explicit linear Runge-Kutta numerical solvers by decreasing the number of slope evaluations in its each application while maintaining the convergence order whereas others in (Qureshi and Ramos, 2018; Ramos et al., 2010; Ramos and Jesús, 2008 ) have devised new nonlinear numerical solvers to deal with the autonomous and non-autonomous continuous dynamical systems having singular or singularly perturbed solutions at some point of the time interval $[0, T]$ . Rational block and trigonometrically fitted solvers have recently been proposed for the efficient solutions of initial value problems having oscillatory solutions, in particular (Ehigie et al., 2017; Jesús et al., 2017; Jesús et al., 2008; Kalogiratou, 2016; Tsitouras and Simos, 2018 ). Apart from this, the most recent work is discussed in (Bakodah, 2018; Liu, 2018; Moutsinga et al., 2018; Neštický et al., 2017; Nuruddeen and Nass, 2017; Qureshi and Emmanuel, 2018; Qureshi et al., 2018; Raza and Khan, 2019 ).

There are many reasons for one numerical solver to become superior over the other. One of the major reasons is the convergence order (a numerical solver is of order r if the error is proved to have $O (h^{r})$ as $h \to 0$ for a constant time step length $h > 0$ ) which shows the rate at which the numerical error in an application of a constant time step length solver along $[0, T]$ decreases as a consequence of decreasing the time step length h. It is therefore a numerical solver having first order convergence would experience halving of the numerical error when the time step length h is halved. Likewise, second order convergent numerical solver would yield the error being multiplied by a factor of $1 / 2^{2}$ when h is halved.

Nevertheless, the numerical errors observed in practice could be different (worse or better) from those predicted and it is because the convergence orders are asymptotic results and apply only when $h \to 0$ . There is always an asymptotic error constant that depends on the solver and the equation which decides the size of the true error. This is the reason for two numerical solvers of same convergence order to have entirely different behavior in errors when they are used for solving the same initial value problem.

The present study aims to propose a new numerical solver having third order convergence. The proposed solver is of explicit linear single-step structure but the inclusion of partial derivative with respect to the second variable $y (t)$ of the slope u makes it well comparable with existing third order convergent solvers of similar characteristics. Thus, it can be included in the family of explicit single-step numerical solvers which are considered to be cost effective and easy to be implemented within a computer code.

2 Derivation of the proposed numerical solver

Suppose that a well-defined initial value problem is given as

(1)

\frac{dy (t)}{dt} = u (t, y (t)), y (t_{0}) = y_{0}, y, u (t, y (t)) \in R, t \in [0, T] \subset R, T > 0 .

Further, it is assumed that the problem (1) possesses a unique continuously differentiable solution $y (t)$ where $y (t_{i}) \approx y_{i}$ for $y_{i}$ is taken to be the approximation to the analytical (exact) solution $y (t)$ at $t = t_{i}$ over the associated interval of integration $[0, T]$ with $h = T / n, (i = 0, 1, \dots n)$ as a fixed time step length used at each step of the integration.

A numerical solver having explicit nature for single-step structure is generally formulated as:

(2)

y_{i + 1} = y_{i} + h Φ_{u} (t_{i}, y_{i}; h),

where

Φ_{u} (t_{i}, y_{i}; h)

can be expressed in terms of Taylor series expansion of an arbitrary function

u (t, y)

. Moreover, the Taylor series expansion of

y (t_{i} + h)

is represented as:

(3)

y (t_{i} + h) = y (t_{i}) + hu + \frac{1}{2!} h^{2} (u_{t} + {uu}_{y}) + \frac{1}{3!} h^{3} (u_{tt} + 2 {uu}_{ty} + u^{2} u_{yy} + {uu}_{y}^{2} + u_{t} u_{y}) + \frac{1}{4!} h^{4} (u_{ttt} + 3 {uu}_{tty} + 3 u^{2} u_{tyy} + 5 {uu}_{y} u_{ty} + 3 u_{t} u_{ty} + u^{3} u_{yyy} + 4 u^{2} u_{y} u_{yy} + 3 {uu}_{t} u_{yy} + {uu}_{y}^{3} + u_{t} u_{y}^{2} + u_{tt} u_{y}) + O (h^{5}) .

Taking into consideration the general structure for the single-step numerical solver, we formulated the required solver as follows:

(4)

y_{i + 1} = y_{i} + h (c_{1} m_{1} + c_{2} m_{2} + c_{3} m_{3}) .

The slopes in (4) are structured as follows:

(5)

\begin{matrix} m_{1} = u (t_{i}, y_{i}), \\ m_{2} = u (t_{i} + p_{2} h, y_{i} + {hm}_{1} (q_{21} + {hr}_{21} u_{y})), \\ m_{3} = u (t_{i} + p_{3} h, y_{i} + h (q_{31} m_{1} + q_{32} m_{2}) + h^{2} r_{31} m_{1} u_{y}) . \end{matrix}

Expanding

m_{2}

and

m_{3}

in Taylor’s series, we obtain

\begin{matrix} m_{2} = u + ({uu}_{y} q_{21} + u_{t} p_{2}) h + \\ (\frac{1}{2} u_{yy} u^{2} q_{21}^{2} + u_{ty} {up}_{2} q_{21} + {uu}_{y}^{2} r_{21} + \frac{1}{2} u_{tt} p_{2}^{2}) h^{2} + \\ (\frac{1}{6} u_{yyy} u^{3} q_{21}^{3} + \frac{1}{2} u_{tyy} u^{2} p_{2} q_{21}^{2} + u_{yy} u^{2} u_{y} q_{21} r_{21} + \\ \frac{1}{2} u_{tty} {up}_{2}^{2} q_{21} + u_{ty} {uu}_{y} p_{2} r_{21} + \frac{1}{6} u_{ttt} p_{2}^{3}) h^{3} + O (h^{4}) . \end{matrix}

\begin{matrix} m_{3} = u + (u_{y} {uq}_{32} + u_{y} {uq}_{31} + u_{t} p_{3}) h + (\frac{1}{2} u_{yy} u^{2} q_{32}^{2} + u_{yy} u^{2} q_{31} q_{32} + \frac{1}{2} u_{yy} u^{2} q_{31}^{2} + \\ u_{ty} {up}_{3} q_{32} + u_{ty} {up}_{3} q_{31} + u_{y}^{2} {ur}_{31} + \frac{1}{2} u_{tt} p_{3}^{2} + u_{y}^{2} {uq}_{21} q_{32} + u_{t} u_{y} p_{2} q_{32}) h^{2} + \\ (\frac{1}{2} u_{yyy} u^{3} q_{31} q_{32}^{2} + \frac{1}{2} u_{yyy} u^{3} q_{31}^{2} q_{32} + \frac{1}{2} u_{tyy} u^{2} p_{3} q_{32}^{2} + \frac{1}{2} u_{tyy} u^{2} p_{3} q_{31}^{2} + \\ \frac{1}{2} u_{tty} {up}_{3}^{2} q_{32} + \frac{1}{2} u_{tty} {up}_{3}^{2} q_{31} + \frac{1}{2} u_{y} u^{2} u_{yy} q_{21}^{2} q_{32} + u_{y}^{3} {uq}_{32} r_{21} + \\ u_{y} {uu}_{ty} p_{2} q_{21} q_{32} + \frac{1}{2} u_{y} u_{tt} p_{2}^{2} q_{32} + \frac{1}{6} u_{yyy} u^{3} q_{32}^{3} + \frac{1}{6} u_{yyy} u^{3} q_{31}^{3} + \\ u_{tyy} u^{2} p_{3} q_{31} q_{32} + u_{yy} u^{2} u_{y} q_{31} r_{31} + u_{ty} {uu}_{y} p_{3} r_{31} + u_{y} u^{2} u_{yy} q_{21} q_{31} q_{32} + \\ u_{t} {uu}_{yy} p_{2} q_{31} q_{32} + \frac{1}{6} u_{ttt} p_{3}^{3} + u_{y} u^{2} u_{yy} q_{21} q_{32}^{2} + u_{t} {uu}_{yy} p_{2} q_{32}^{2} + \\ u_{y} {uu}_{ty} p_{3} q_{21} q_{32} + u_{t} u_{ty} p_{2} p_{3} q_{32} + u_{yy} u^{2} u_{y} q_{32} r_{31}) h^{3} + O (h^{4}) . \end{matrix}

Substituting the result of

m_{1}, m_{2}

and

m_{3}

into

(4)

and equating the coefficients of powers of h up to

h^{3}

with that of (3), one obtains:

\begin{matrix} c_{1} + c_{2} + c_{3} = 1 & p_{2} c_{3} q_{32} = \frac{1}{6} & p_{2} c_{2} + p_{3} c_{3} = \frac{1}{2} \\ \frac{1}{2} (p_{2}^{2} c_{2} + p_{3}^{2} c_{3}) = \frac{1}{6} & c_{2} q_{21} + c_{3} q_{31} + c_{3} q_{32} = \frac{1}{2} & p_{2} c_{2} q_{21} + p_{3} c_{3} q_{31} + p_{3} c_{3} q_{32} = \frac{1}{3} \\ c_{2} r_{21} + c_{3} r_{31} + c_{3} q_{21} q_{32} = \frac{1}{6} & \frac{1}{2} (c_{2} q_{21}^{2} + c_{3} q_{31}^{2} + c_{3} q_{32}^{2}) + c_{3} q_{31} q_{32} = \frac{1}{6} . \end{matrix}

This non-linear system contains eight equations and ten variables which has been solved using MATLAB to obtain one of its optimal solutions and finally we obtained the following numerical solver with third order accuracy:

(6)

\begin{matrix} m_{1} = u (t_{i}, y_{i}), \\ m_{2} = u (t_{i} + \frac{2}{3} h, y_{i} + {hm}_{1} (\frac{2}{3} + h \frac{1}{2} u_{y})), \\ m_{3} = u (t_{i} + \frac{2}{3} h, y_{i} + h (- \frac{5}{6} m_{1} + \frac{3}{2} m_{2}) - h^{2} \frac{7}{4} m_{1} u_{y}), \\ y_{i + 1} = y_{i} + \frac{h}{12} (3 m_{1} + 7 m_{2} + 2 m_{3}) . \end{matrix}

3 Analysis for the proposed solver

For the required analysis, some definitions and the related theorems have been stated.

3.1

3.1 Consistency

Definition 1

A numerical solver used for finding the solution of an initial value problem in ODEs is said to be consistent if

(7)

\lim_{h \to 0} Φ_{u} (t, y; h) = u (t, y (t)),

where

Φ_{u}

is the increment function for the numerical solver.

The increment function of the proposed solver is as follows:

(8)

Φ_{u} (t, y; h) = \frac{1}{4} (3 m_{1} + 7 m_{2} + 2 m_{3}),

which clearly satisfies the condition (7) of the above definition for the consistency. Thus, the proposed solver is found to be consistent with the given initial value problem (1).

3.2

3.2 Stability

Theorem 1

Suppose that $λ_{0}, λ_{1}, \dots, λ_{n}$ be real numbers satisfying

(9)

| λ_{j + 1} | ⩽ (1 + α) | λ_{j} | + β, α > 0, β ⩾ 0, j = 0, 1, 2, \dots n - 1 .

Thus,

(10)

| λ_{n} | ⩽ \exp (n α) | λ_{0} | + \frac{\exp (n α) - 1}{α} β .

Proof

Using the assumption presented in the above Theorem 1, we can write the following:

(11)

\begin{matrix} | λ_{1} | ⩽ (1 + α) | λ_{0} | + β, \\ | λ_{2} | ⩽ {(1 + α)}^{2} | λ_{0} | + (1 + α) β + β, \\ | λ_{3} | ⩽ {(1 + α)}^{3} | λ_{0} | + {(1 + α)}^{2} β + (1 + α) β + β, \\ ⋮ \\ | λ_{n} | ⩽ {(1 + α)}^{n} | λ_{0} | + \frac{{(1 + α)}^{n} - 1}{α} β, \\ ⩽ \exp (n α) | λ_{0} | + \frac{\exp (n α) - 1}{α} β, for α > - 1; 0 < 1 + α ⩽ \exp (α) . \end{matrix}

Theorem 2

Suppose that $y_{i + 1}$ and ${\hat{y}}_{i + 1}$ be two numerical solutions to the ODE $y^{'} (t) = u (t, y (t))$ produced by a numerical solver with prescribed initial conditions $y (t_{0}) = y_{0}$ and $\hat{y} (t_{0}) = \hat{y_{0}}$ respectively such that $| y_{0} - \hat{y_{0}} | < ∊, ∊ > 0$ . The following condition

(12)

| y_{i + 1} - {\hat{y}}_{i + 1} | ⩽ M | y_{0} - {\hat{y}}_{0} |, M > 0,

is considered to be the necessary and sufficient condition for the numerical solver to be stable.

Proof

Using the above Theorem 2, the proposed solver would give us the following:

(13)

\begin{matrix} y_{i + 1} - y_{i} = h Φ_{u} (t_{i}, y_{i}; h), \\ {\hat{y}}_{i + 1} - {\hat{y}}_{i} = h Φ_{u} (t_{i}, {\hat{y}}_{i}; h) . \end{matrix}

Using the triangle inequality, we have

(14)

\begin{matrix} | y_{i + 1} - {\hat{y}}_{i + 1} | ⩽ | y_{i} - {\hat{y}}_{i} | + h | Φ_{u} (t_{i}, y_{i}; h) - Φ_{u} (t_{i}, {\hat{y}}_{i}; h) |, \\ ⩽ | y_{i} - {\hat{y}}_{i} | + h | y_{i} - {\hat{y}}_{i} |, \\ = (1 + \hat{L} h) | y_{i} - {\hat{y}}_{i} |, \\ = {(1 + \hat{L} h)}^{i} | y_{0} - {\hat{y}}_{0} | . \end{matrix}

Employing the Theorem 1 while taking

α = \hat{L} h

and

β = 0

, we obtain the following

(15)

| y_{i + 1} - {\hat{y}}_{i + 1} | ⩽ M | y_{0} - {\hat{y}}_{0} |, M = \exp (\hat{L} hi) .

This confirms the proof for the proposed solver to be stable.

Further, the Dahlquist (Dahlquist, 1956) test problem $y^{'} (t) = μ y (t), y (0) = 1, Re (μ) < 0$ gives the following stability function for the proposed solver.

(16)

Ψ (z) = 1 + \frac{11}{12} z + \frac{1}{2} z^{2} + \frac{1}{6} z^{3} + \frac{1}{8} z^{4}, z = μ h .

3.3

3.3 Error analysis

Theorem 3

Suppose that $(t_{i}, y_{i})$ and $(t_{i}, {\hat{y}}_{i})$ be any two points in the region S defined by $S = {(t, y) \in R^{2} | t_{0} ⩽ t ⩽ t_{n}, - \infty < y < \infty}$ , and that u is a Lipschitz continuous function on S such that $| u (t_{i}, y_{i}) - u (t_{i}, {\hat{y}}_{i}) | ⩽ L | y_{i} - {\hat{y}}_{i} |$ , then the increment function $Φ$ is Lipschitz continuous, and $| Φ (t_{i}, y_{i}; h) - Φ (t_{i}, {\hat{y}}_{i}; h) | ⩽ \hat{L} | y_{i} - {\hat{y}}_{i} |$ where L and $\hat{L}$ are, respectively, the Lipschitz constants for u and $Φ$ .

Proof

Employing the assumptions above, we get the following

(17)

| m_{1} - \hat{m_{1}} | ⩽ L | y_{i} - {\hat{y}}_{i} | .

(18)

| m_{2} - \hat{m_{2}} | ⩽ (1 + \frac{2}{3} hL + \frac{1}{2} h^{2} L^{2}) | y_{i} - {\hat{y}}_{i} |, ∵ | u_{y} | ⩽ L .

(19)

| m_{3} - \hat{m_{3}} | ⩽ (1 + \frac{2}{3} hL - \frac{3}{4} h^{2} L^{2} + \frac{3}{4} h^{3} L^{3}) | y_{i} - {\hat{y}}_{i} | .

Thus, the increment function

Φ

for the proposed solver becomes

(20)

| Φ_{u} (t, y; h) - Φ_{u} (t, \hat{y}; h) | ⩽ \hat{L} | y_{i} - {\hat{y}}_{i} |,

where

\hat{L} = \frac{3}{4} + \frac{1}{4} L + \frac{1}{2} hL + \frac{1}{6} h^{2} L^{2} + \frac{1}{8} h^{3} L^{3}

This completes the proof for the increment function $Φ_{u}$ for the proposed solver to be Lipchitz continuous.

3.3.1

3.3.1 Error bounds

Definition 2

Suppose that $y (t_{i})$ be the analytical (exact) solution of the underlying initial value problem at $t = t_{0}$ and $y_{i}$ be its corresponding numerical (approximate) solution obtained via a numerical solver. The global error $Ω_{i}$ at $t = t_{i}$ is defined as:

(21)

Ω_{i} = y (t_{i}) - y_{i} .

Definition 3

The local error $Ξ_{i + 1}$ at $t = t_{i + 1}$ for a numerical solver like the one being studied in the present study is defined as:

(22)

Ξ_{i + 1} = y (t_{i + 1}) - y_{i + 1},

where it should be assumed that the history of the solver is free from any error, that is,

y (t_{i}) = y_{i}

Definition 4

The local error $Ξ_{i + 1}$ of a numerical solver with order of accuracy r is also defined as

(23)

Ξ_{i + 1} = Ψ (t_{i}, y_{i}) h^{r + 1} + O (h^{r + 2}),

where the first term in

Ξ_{i + 1}

is known as the leading local truncation error.

In order to obtain the local truncation error of the proposed integrator, a usual functional associated to the integrator has been considered, that is given below:

(24)

Ξ (v (t), h) = v (t + h) - y_{i + 1},

where

v (t)

is an arbitrary function defined along the integration interval

[0, T]

and differentiable as many times as required. Having expanded it into Taylor series about t and collecting the terms in h, the local truncation error under local assumption of the following form has been obtained that ensures at least third order accuracy of the proposed solver:

(25)

Ξ_{i + 1} = \frac{h^{4}}{216} (3 {uu}_{tty} + 3 u_{t} u_{ty} + u^{2} u_{tyy} + {uu}_{t} uyy - {uu}_{y} u_{ty} - u_{tt} u_{y} - 18 {uu}_{y}^{3} + 9 u_{t} u_{y}^{2} + u^{3} u_{yyy} + u_{ttt}) + O (h^{5}) .

Lotkin in (Lotkin, 1951) introduced the concept of determining the local bounds for the function u and its partial derivatives for

t \in [0, T]

and

y \in (- \infty, \infty)

as follows:

(26)

| u (t, y) | < P, | \frac{\partial^{l + m} u}{\partial t^{l} \partial u^{m}} | < \frac{Q^{l + m}}{P^{m - 1}}, (l + m) ⩽ r,

where P and Q are some positive real numbers and r is the order of accuracy of the numerical solver. Using this idea of the Lotkin’s bounds, a bound for the leading local error of the proposed solver is determined to be as follows:

(27)

| Ψ (t_{i}, y_{i}) | < \frac{47}{216} h^{4} {PQ}^{3} .

A bound for the leading local error of a single step explicit linear numerical solver is also considered to be a bound for its local truncation error. Using the fact, the local error bound for the proposed solver is finalized as follows:

(28)

| Ξ_{i + 1} | < {Ah}^{4}, A = \frac{47}{216} {PQ}^{3} .

3.4

3.4 Convergence

Theorem 4

The convergence for a numerical solver is guaranteed if for every initial value problem $y^{'} (t) = u (t, y (t)), y (t_{0}) = y_{0}$ with u being a Lipchitz continuous function, one obtains

(29)

\lim_{h \to 0} \sup_{0 ⩽ i ⩽ N} | | y (t_{i}) - y_{i} | | \to 0, \forall t > 0 .

Proof

Let us consider the proposed solver to be

(30)

y_{i + 1} - y_{i} - h Φ_{u} (t_{i}, y_{i}; h) = 0 .

Using the definition of the local error in conjunction with the above theorem (4), one obtains

(31)

| Ω_{i + 1} | ⩽ | Ω_{i} | + h | Φ_{u} (t_{i}, y (t_{i}); h) - Φ_{u} (t_{i}, y_{i}; h) | + | Ξ_{i} |,

where

Ω_{i}

is the global error at

t = t_{i}

Employing the Theorems 1 and 3, we obtain

(32)

| Ω_{i + 1} | ⩽ (1 + h \hat{L}) | Ω_{i} | + | Ξ_{i} | .

Suppose that

R = \sup_{0 ⩽ i ⩽ N} | Ξ_{i} |

, then from Theorem 1 one can obtain

(33)

| Ω_{i} | ⩽ \exp (i \hat{L} h) | Ω_{0} | + \frac{R}{\hat{L} h} (\exp (i \hat{L} h) - 1) .

This finally yielded the global error bound for the proposed solver on the basis of having

Ω_{0} = 0

and

h = t_{i} / i

(34)

| Ω_{i} | ⩽ \frac{R}{\hat{L} h} (\exp (\hat{L} t_{i}) - 1),

where R is the local error bound as found in (28).

Thus, the convergence of the proposed solver is verified for

(35)

\lim_{h \to 0} | Ω_{i} | = 0 \Rightarrow \lim_{h \to 0} | y (t_{i}) - y_{i} | = 0 .

4 Numerical experiments

Example 1. As shown in the Table 1 where the proposed solver shows the three-order decreasing behavior in the maximum absolute global relative errors $(E_{\max} = \max_{0 ⩽ i ⩽ T} | y (t_{i + 1}) - y_{i + 1} |)$ and in the errors computed at final nodal point $(E (t = T) = | y (T) - y_{N} |)$ over $[0, T]$ for every one-order decrease in the time step length h. The $ℓ^{2}$ -error norm $(ℓ^{2} = \sqrt{\sum_{k = 0}^{n} | y (t_{i + 1}) - y_{i + 1} |^{2}})$ is shown in the last row. In addition, the computed errors are smallest in every case upon comparison with Heun and RK3HM (Wazwaz, 1990).

(36)

\begin{matrix} y^{'} (t) = ty {(t)}^{3} - y (t), y (0) = 1, t \in [0, 2] \\ y (t) = \frac{2}{\sqrt{2 + 4 t + 2 \exp (2 t)}} \end{matrix}

Table 1 Maximum absolute global relative errors on

[0, 2]

(first row), absolute relative errors at

t = 2

(second row), and

ℓ^{2}

-error norm (third row) for the Example-1.

h	$10^{- 1}$	$10^{- 2}$	$10^{- 3}$	$10^{- 4}$
Heun	1.3048e−04 1.3048e−04 4.2260e−04	1.2425e−07 1.2425e−07 1.2441e−06	1.2352e−10 1.2352e−10 3.9015e−09	1.5984e−13 1.5970e−13 1.5036e−11

RK3HM	6.6502e−04 6.6502e−04 2.2030e−03	5.9464e−06 5.9464e−06 6.0417e−05	5.9024e−08 5.9024e−08 1.8898e−06	5.8979e−10 5.8979e−10 5.9693e−08

Proposed	2.3861e−05 8.2608e−06 8.1340e−05	2.6075e−08 1.3196e−08 2.8703e−07	2.6284e−11 1.3664e−11 9.1636e−10	5.8001e−14 4.9848e−14 5.6582e−12

Example 2. The observation for the linear problem with exponentially increasing solution based upon the Table 2 is more or less same as was in the former example.

(37)

\begin{matrix} y^{'} (t) = t^{2} y (t), y (0) = 1, t \in [0, 1] \\ y (t) = \exp (t^{3} / 3) \end{matrix}

Table 2 Maximum absolute global relative errors on

[0, 1]

(first row), absolute relative errors at

t = 1

(second row), and

ℓ^{2}

-error norm (third row) for the Example-2.

h	$10^{- 1}$	$10^{- 2}$	$10^{- 3}$	$10^{- 4}$
Heun	6.4697e−05 6.4697e−05 8.3000e−05	6.8998e−08 6.8998e−08 2.2883e−07	6.9400e−11 6.9400e−11 7.1171e−10	1.3111e−13 1.3110e−13 4.4807e−12

RK3HM	2.2335e−03 2.2335e−03 4.1400e−03	2.4716e−05 2.4716e−05 1.3587e−04	2.4971e−07 2.4971e−07 4.3116e−06	2.4998e−09 2.4998e−09 1.3639e−07

Proposed	2.0183e−05 2.0183e−05 2.8573e−05	1.8702e−08 1.8702e−08 7.7040e−08	1.8535e−11 1.8535e−11 2.3974e−10	7.3701e−14 7.3664e−14 2.6931e−12

Example 3. The third initial value problem is structured in the form of a well-known nonlinear first order Riccati ODE. The proposed solver as shown in the Table 3 provides smallest magnitude of errors when compared with Heun and RK3HM methods. With one-order decrease in the time step length h, there are three-order decrease in the maximum absolute global relative errors and in the absolute relative errors determined at the final nodal point over $[0, 0.5]$ .

(38)

\begin{matrix} y^{'} (t) = \frac{2 \cos^{2} (t) - \sin^{2} (t) + y^{2} (t)}{2 \cos (t)}, y (0) = - 1, t \in [0, 0.5] \\ y (t) = \sin (t) - \frac{1}{0.5 \sin (t) + \cos (t)} . \end{matrix}

Table 3 Maximum absolute global relative errors on

[0, 0.5]

(first row), absolute relative errors at

t = 0.5

(second row), and

ℓ^{2}

-error norm (third row) for the Example-3.

h	$10^{- 1}$	$10^{- 2}$	$10^{- 3}$	$10^{- 4}$
Heun	5.4644e−05 5.4644e−05 7.9834e−05	5.1896e−08 5.1896e−08 2.0943e−07	5.1603e−11 5.1603e−11 6.4899e−10	5.9233e−14 5.9172e−14 2.1293e−12

RK3HM	3.1635e−04 3.1635e−04 3.5014e−04	3.5178e−06 3.5178e−06 9.1739e−06	3.5477e−08 3.5477e−08 2.8310e−07	3.5505e−10 3.5505e−10 8.9296e−09

Proposed	6.4731e−06 3.2754e−06 1.0836e−05	8.3861e−09 1.9656e−09 4.2872e−08	8.3674e−12 2.1622e−12 1.3480e−10	1.2274e−14 1.2155e−14 4.2406e−13

Example 4. Here, the errors are computed for each unknown function $x (t)$ and $y (t)$ of the system as shown in the Tables 4–6 where decreasing behavior is clearly seen for all the solvers with the proposed solver having the smallest errors almost for every time step length h.

(39)

\begin{matrix} x^{'} (t) = x (t) - 10 y (t), x (0) = 0, \\ y^{'} (t) = 15 x (t) + y (t), y (0) = 1 \\ x (t) = - \sqrt{\frac{2}{3}} \exp (t) \sin (5 \sqrt{6} t) \\ y (t) = \exp (t) \cos (5 \sqrt{6} t) . \end{matrix}

Table 4 Maximum absolute global relative errors on

[0, 10]

for each unknown in the Example-4.

h	$10^{- 1}$	$10^{- 2}$	$10^{- 3}$	$10^{- 4}$
Heun	9.3411e+01 2.0974e+01	8.1516e−01 3.9783e+00	1.8376e−03 4.0405e−02	1.9256e−04 6.1855e−05

RK3HM	2.5717e+07 4.1029e+07	1.5990e+01 9.2617e+01	1.3092e+00 2.4304e+01	1.3639e+00 4.2652e−01

Proposed	2.6958e+01 1.7027e+01	3.6696e−01 1.7810e+00	5.1390e−04 1.1450e−02	5.3860e−05 1.7773e−05

Table 5 Absolute relative errors at

t = 10

for each unknown in the Example-4.

h	$10^{- 1}$	$10^{- 2}$	$10^{- 3}$	$10^{- 4}$
Heun	3.5681e+00 1.0769e+00	8.9169e−02 8.3767e−03	7.4569e−05 8.8139e−06	7.5540e−08 8.8355e−09

RK3HM	2.5571e+07 1.8402e+07	2.7526e+00 1.8522e−01	4.8144e−02 1.6401e−03	4.9232e−04 2.1372e−05

Proposed	1.3021e+00 1.0083e+00	4.5361e−02 9.2884e−03	2.8136e−05 9.5089e−06	2.8907e−08 9.5067e−09

Table 6

ℓ^{2}

-error norms on

[0, 10]

for each unknown in the Example-4 where NaN stands for Not a Number.

h	$10^{- 1}$	$10^{- 2}$	$10^{- 3}$	$10^{- 4}$
Heun	NaN 2.9205e+01	NaN 4.1695e+00	NaN 4.3136e−02	NaN 9.9840e−05

RK3HM	NaN 4.9949e+07	NaN 9.6513e+01	NaN 2.6034e+01	NaN 6.9477e−01

Proposed	NaN 2.1094e+01	NaN 1.8718e+00	NaN 1.2236e−02	NaN 2.7782e−05

Example 5. The angular displacement $x (t)$ for an undamped pendulum with a driving force $f (t) = \cos (4 t)$ can be modeled in terms of a second order ODE as shown below. (See Fig 1)

(40)

x^{″} (t) + \sin (x) = \cos (4 t), x (0) = 1, x^{'} (0) = 0 .

(41)

\begin{matrix} x_{1}^{'} (t) = x_{2} (t), x_{1} (0) = 1, \\ x_{2}^{'} (t) = - \sin (x_{1}) + \cos (4 t), x_{2} (0) = 0, \\ t \in [0, 20] . \end{matrix}

Due to the nonavailability of the exact solution for the above nonlinear dynamical system, we have made some graphical comparisons of the results obtained through the numerical solvers with that of the numerical solution obtained through ode45 solver while taking the integration steps

N = 100

and

N = 200

. However, for

N = 200

it is difficult to observe the discrepancy among the curves for which a zooming plot at the bottom right corner for the Figs. 4 and 5 help to see the performance of the solvers over a small domain. Upon closer examination of the Figs. 2–5 , it can be seen that the proposed and Heun solvers agree well with each other.

Stability region (shaded) for the proposed solver with | Ψ ( z ) | ⩽ 1 .

Comparison of solvers with N = 100 for the angular displacement x 1 ( t ) in the Example – 5.

Comparison of solvers with N = 100 for the angular velocity x 2 ( t ) in the Example – 5.

Comparison of solvers with N = 200 for the angular displacement x 1 ( t ) in the Example – 5.

Comparison of solvers with N = 200 for the angular velocity x 2 ( t ) in the Example – 5.

5 Conclusion

A numerical solver (6) having third order convergence has been proposed in the present study. Its error analysis led us to compute both local and global error bounds wherein the leading term $O (h^{4})$ in the local truncation error ensures the third order convergence of the solver. The increment function of the solver is shown to be Lipchitz continuous and the linear stability analysis has been carried out to get the stability region of the solver satisfying the condition $| Ψ (z) | ⩽ 1$ .

Upon comparison with two well-known solvers (Heun and RK3HM) having similar characteristics as that of the proposed solver, it has been observed that the latter yielded smallest errors in both scalar and continuous dynamical systems of first and higher order as shown in the Tables 1–6 . Thus, the proposed solver can be considered a member for the family of explicit linear single-step numerical solvers used to solve various continuous dynamical systems in the scientific fields.

Conflict of Interest

The authors declare no conflict of interest. Both authors have read and approved the final version of the present research work.

Acknowledgement

Sania Qureshi is grateful to Mehran University of Engineering and Technology, Jamshoro, Pakistan for the kind support and facilities provided to carry out this research work.

References

Amann, H., 2011. Ordinary differential equations: an introduction to nonlinear analysis. Walter degruyter.
Bakodah H.O., . An efficient modification of the decomposition method with a convergence parameter for solving Korteweg de Vries equations. J. King Saud 2018 University-Science
[Google Scholar]
Butcher J.C., . Numerical Methods for Ordinary Differential Equations. John Wiley & Sons; 2016.
Liu Chein-Shan, . Solving singular convection-diffusion equation by exponentially-fitted trial functions and adjoint Trefftz test functions. J. King Saud Univ. - Sci.. 2018;30:100-105.
[Google Scholar]
Dahlquist G., . Convergence and stability in the numerical integration of ordinary differential equations. Math. Scandinavica. 1956;4:33-53.
[Google Scholar]
Ehigie J.O., Jator S.N., Okunuga S.A., . A multi-point numerical integrator with trigonometric coefficients for initial value problems with periodic solutions. Numer. Analys. Appl.. 2017;10:272-286.
[Google Scholar]
Hirsch M.W., Stephen S., Robert L.D., . Differential Equations, Dynamical Systems, and an Introduction to Chaos. Academic Press; 2012.
Jesús V.A., Ramos H., Carmelo C., . A first approach in solving initial-value problems in ODEs by elliptic fitting methods. J. Comput. Appl. Math.. 2017;318:599-603.
[Google Scholar]
Jesús V.A., Vaquero J.M., Ramos H., . Exponential fitting BDF-Runge-Kutta algorithms. Comput. Phys. Commun.. 2008;178:15-34.
[Google Scholar]
Jordan D.W., Peter S., . Nonlinear Ordinary Differential Equations: An Introduction to Dynamical Systems. Oxford University Press; 1999.
Kalogiratou Z., . A new approach on the construction of trigonometrically fitted two step hybrid methods. J. Comput. Appl. Math.. 2016;303:146-155.
[Google Scholar]
Lotkin M., . On the accuracy of Runge-Kutta’s method. Math. Tables Other Aids Comput.. 1951;5:128-133.
[Google Scholar]
Ma J., Simos T.E., . An efficient and computational effective method for second order problems. JOMC. 2017;55:1649-1668.
[Google Scholar]
Mermoud, O., Lausanne B.Z.M., 2015. Ordinary Differential Equations and Dynamical Systems.
Moutsinga C.R.B., . Homotopy perturbation transform method for pricing under pure diffusion models with affine coefficients. J. King Saud Univ.-Sci.. 2018;30:1-13.
[Google Scholar]
Neštický M., . Analytical and numerical approach to defining the stability region of a dynamical system. J. Appl. Math., Stat. Inf.. 2017;13:55-62.
[Google Scholar]
Nuruddeen R.I., Nass A.M., . Exact solutions of wave-type equations by the Aboodh decomposition method. Stochastic Model. Appl.. 2017;21:23-30.
[Google Scholar]
Qureshi S., Ramos H., . L-stable explicit nonlinear method with constant and variable step-size formulation for solving initial value problems. Int. J. Nonlinear Sci. Numer. Simul.. 2018;19:741-751.
[Google Scholar]
Qureshi S., Emmanuel F.S., . Convergence of a numerical technique via interpolating function to approximate physical dynamical systems. J. Adv. Phys.. 2018;7:446-450.
[Google Scholar]
Qureshi S., Memon Z.A., Shaikh A.A., . Local accuracy and error bounds of the improved Runge-Kutta numerical methods. J. Appl. Math. Comput. Mech.. 2018;17:73-84.
[Google Scholar]
Rabiei F., Fudziah I., . Fifth-order improved Runge-Kutta methods with reduced number of function evaluations. Aust. J. Basic Appl. Sci.. 2012;6:97-105.
[Google Scholar]
Rabiei F., . Accelerated Runge-Kutta Nystrom method for solving autonomous second-order ordinary differential equations. WASJ. 2012;17:1549-1555.
[Google Scholar]
Rabiei F., . Construction of improved Runge-Kutta Nystrom method for solving second-order ordinary differential equations. WASJ. 2012;20:1685-1695.
[Google Scholar]
Ramos H., . Numerical solution of nonlinear singularly perturbed problems on nonuniform meshes by using a non-standard algorithm. JOMC. 2010;48:38-54.
[Google Scholar]
Ramos H., Jesús V.A., . A new algorithm appropriate for solving singular and singularly perturbed autonomous initial-value problems. Int. J. Comput. Math.. 2008;85:603-611.
[Google Scholar]
Raza A., Khan A., . Haar wavelet series solution for solving neutral delay differential equations. J. King Saud Univ.-Sci.. 2019;31:1070-1076.
[Google Scholar]
Shampine L.F., . Numerical Solution of Ordinary Differential Equations. Routledge; 2018.
Stickler B.A., Ewald S., . Ordinary Differential Equations: Initial Value Problems. Basic Concepts in Computational Physics. Cham: Springer; 2014. p. :61-79.
Teschl G., . Ordinary differential equations and dynamical systems. University of Vienna; 2000. Lecture Notes
Tsitouras C., Simos T.E., . Trigonometric-fitted explicit numerov-type method with vanishing phase-lag and its first and second derivatives. MedJM. 2018;15:168.
[Google Scholar]
Wazwaz A.M., . A modified third order Runge-Kutta method. Appl. Math. Lett.. 1990;3:123-125.
[Google Scholar]

Introduction

Derivation of the proposed numerical solver

Analysis for the proposed solver

Consistency
Stability
Error analysis
Convergence

Numerical experiments

Conclusion

Conflict of Interest

Fulltext Views
27

PDF downloads
23

View/Download PDF
Download Citations

BibTeX
RIS

Show Sections

[1] Amann, H., 2011. Ordinary differential equations: an introduction to nonlinear analysis. Walter degruyter.

[2] Bakodah H.O., . An efficient modification of the decomposition method with a convergence parameter for solving Korteweg de Vries equations. J. King Saud 2018 University-Science
[Google Scholar]

[3] Butcher J.C., . Numerical Methods for Ordinary Differential Equations. John Wiley & Sons; 2016.

[4] Liu Chein-Shan, . Solving singular convection-diffusion equation by exponentially-fitted trial functions and adjoint Trefftz test functions. J. King Saud Univ. - Sci.. 2018;30:100-105.
[Google Scholar]

[5] Dahlquist G., . Convergence and stability in the numerical integration of ordinary differential equations. Math. Scandinavica. 1956;4:33-53.
[Google Scholar]

[6] Ehigie J.O., Jator S.N., Okunuga S.A., . A multi-point numerical integrator with trigonometric coefficients for initial value problems with periodic solutions. Numer. Analys. Appl.. 2017;10:272-286.
[Google Scholar]

[7] Hirsch M.W., Stephen S., Robert L.D., . Differential Equations, Dynamical Systems, and an Introduction to Chaos. Academic Press; 2012.

[8] Jesús V.A., Ramos H., Carmelo C., . A first approach in solving initial-value problems in ODEs by elliptic fitting methods. J. Comput. Appl. Math.. 2017;318:599-603.
[Google Scholar]

[9] Jesús V.A., Vaquero J.M., Ramos H., . Exponential fitting BDF-Runge-Kutta algorithms. Comput. Phys. Commun.. 2008;178:15-34.
[Google Scholar]

[10] Jordan D.W., Peter S., . Nonlinear Ordinary Differential Equations: An Introduction to Dynamical Systems. Oxford University Press; 1999.

[11] Kalogiratou Z., . A new approach on the construction of trigonometrically fitted two step hybrid methods. J. Comput. Appl. Math.. 2016;303:146-155.
[Google Scholar]

[12] Lotkin M., . On the accuracy of Runge-Kutta’s method. Math. Tables Other Aids Comput.. 1951;5:128-133.
[Google Scholar]

[13] Ma J., Simos T.E., . An efficient and computational effective method for second order problems. JOMC. 2017;55:1649-1668.
[Google Scholar]

[14] Mermoud, O., Lausanne B.Z.M., 2015. Ordinary Differential Equations and Dynamical Systems.

[15] Moutsinga C.R.B., . Homotopy perturbation transform method for pricing under pure diffusion models with affine coefficients. J. King Saud Univ.-Sci.. 2018;30:1-13.
[Google Scholar]

[16] Neštický M., . Analytical and numerical approach to defining the stability region of a dynamical system. J. Appl. Math., Stat. Inf.. 2017;13:55-62.
[Google Scholar]

[17] Nuruddeen R.I., Nass A.M., . Exact solutions of wave-type equations by the Aboodh decomposition method. Stochastic Model. Appl.. 2017;21:23-30.
[Google Scholar]

[18] Qureshi S., Ramos H., . L-stable explicit nonlinear method with constant and variable step-size formulation for solving initial value problems. Int. J. Nonlinear Sci. Numer. Simul.. 2018;19:741-751.
[Google Scholar]

[19] Qureshi S., Emmanuel F.S., . Convergence of a numerical technique via interpolating function to approximate physical dynamical systems. J. Adv. Phys.. 2018;7:446-450.
[Google Scholar]

[20] Qureshi S., Memon Z.A., Shaikh A.A., . Local accuracy and error bounds of the improved Runge-Kutta numerical methods. J. Appl. Math. Comput. Mech.. 2018;17:73-84.
[Google Scholar]

[21] Rabiei F., Fudziah I., . Fifth-order improved Runge-Kutta methods with reduced number of function evaluations. Aust. J. Basic Appl. Sci.. 2012;6:97-105.
[Google Scholar]

[22] Rabiei F., . Accelerated Runge-Kutta Nystrom method for solving autonomous second-order ordinary differential equations. WASJ. 2012;17:1549-1555.
[Google Scholar]

[23] Rabiei F., . Construction of improved Runge-Kutta Nystrom method for solving second-order ordinary differential equations. WASJ. 2012;20:1685-1695.
[Google Scholar]

[24] Ramos H., . Numerical solution of nonlinear singularly perturbed problems on nonuniform meshes by using a non-standard algorithm. JOMC. 2010;48:38-54.
[Google Scholar]

[25] Ramos H., Jesús V.A., . A new algorithm appropriate for solving singular and singularly perturbed autonomous initial-value problems. Int. J. Comput. Math.. 2008;85:603-611.
[Google Scholar]

[26] Raza A., Khan A., . Haar wavelet series solution for solving neutral delay differential equations. J. King Saud Univ.-Sci.. 2019;31:1070-1076.
[Google Scholar]

[27] Shampine L.F., . Numerical Solution of Ordinary Differential Equations. Routledge; 2018.

[28] Stickler B.A., Ewald S., . Ordinary Differential Equations: Initial Value Problems. Basic Concepts in Computational Physics. Cham: Springer; 2014. p. :61-79.

[29] Teschl G., . Ordinary differential equations and dynamical systems. University of Vienna; 2000. Lecture Notes

[30] Tsitouras C., Simos T.E., . Trigonometric-fitted explicit numerov-type method with vanishing phase-lag and its first and second derivatives. MedJM. 2018;15:168.
[Google Scholar]

[31] Wazwaz A.M., . A modified third order Runge-Kutta method. Appl. Math. Lett.. 1990;3:123-125.
[Google Scholar]

A new third order convergent numerical solver for continuous dynamical systems

Abstract

Abstract

Keywords

Differential equations

Local error

Global error

Lipchitz continuity

Stability

1 Introduction

2 Derivation of the proposed numerical solver

3 Analysis for the proposed solver

3.1 Consistency

3.2 Stability

3.3 Error analysis

3.3.1 Error bounds

3.4 Convergence

4 Numerical experiments

5 Conclusion

Conflict of Interest

Acknowledgement

References

Suggested read for related articles: