Numerical Simulations of Some Real-Life Problems Governed by ODEs

Abstract

In this chapter, some real-life model problems that can be formulated as ordinary differential equations (ODEs) are introduced and numerically studied. These models are the variable-order fractional Hodgkin–Huxley model of neuronal excitation (VOFHHM) and other models with the variable-order fractional (VOF) time delay, such as the 4-year life cycle of a population of lemmings model, the enzyme kinetics with an inhibitor molecule model, and the Chen system model. A class of numerical methods is used to study the above-mentioned models such as non-standard finite difference (NSFD) and Adams-Bashforth-Moulton (ABM) methods. Numerical test examples are presented.

1. Numerical simulations for VOFHHM of neuronal excitation

1.1. Introduction

Hodgkin and Huxley proposed the famous Hodgkin–Huxley equations which quantitatively describe the generation of action potential of squid giant axon [1].

$C_m\frac{d{V}_m}{dt}=I_{ext}-I_{ionic},$E1

where C_m is the membrane capacitance, V_m is the intracellular potential, t is the time, Iionic is the net ionic current flowing across the membrane, and I_ext is the externally applied current [2].

The ionic current described in Hodgkin and Huxley’s giant squid axon can be divided into three different cases: sodium current, potassium current, and leakage current, for example, voltage-gated persistent K⁺ current with four activation gates (resulting in the term n⁴, where n is the activation variable for K⁺), voltage-gated transient Na⁺ current with three activation gates and one inactivation gate (the term m³ ħ) [3].

The total ionic current is therefore

$I_{ionic}={\displaystyle \stackrel{I_{Na}}{{\displaystyle \overbrace{G_{Na}(V-E_{Na})}}}}+{\displaystyle \stackrel{I_K}{{\displaystyle \overbrace{G_K(V-E_K)}}}}+{\displaystyle \stackrel{I_L}{{\displaystyle \overbrace{G_L(V-E_L)}}}},$E2

where I_Na and I_K are the currents through Na⁺ and K⁺ channels, respectively. The current I_L is the leak current and denotes all the residue currents through a cell membrane other than Na⁺ and K⁺ currents. G denotes conductance value for each individual ionic component, and E is an equilibrium potential.

Hodgkin and Huxley defined the probability for a gate that can be found in its permissive state to depend on the membrane potential, therefore incorporating the voltage-dependent conductance [4, 5]. A probability p ranging from 0 to 1 can be defined for each gate, but since a large population of channels (and therefore gates) has to be considered, this leads to the following:

$1-p(t)\rightleftharpoons {}_{{\alpha }_p}{}^{{\beta }_p}p(t)$

where 1−p(t) is the fraction in non-permissive state, p(t) is the fraction in permissive state, and α_p and β_p are called rate functions.

Transitions between permissive and non-permissive states are assumed to obey the relation:

$\frac{dp}{dt}={\alpha }_p\left(V\right)\left(1-p\right)-{\beta }_p\left(V\right)p,p=m,n,\hslash ,$E3

where m and n are the activation variables for the sodium and the potassium, respectively, and represent the sodium de-inactivation.

In the literature, Alan Lloyd Hodgkin and Andrew Huxley authored a succession of five papers describing the nonlinear ODEs that model how action potentials can be initiated and propagated through an axon [1–5]. HHM can be considered one of the most successful mathematical models in quantity lively describing biological phenomena, and it is the method which can be used in deriving the model of a squid that is directly applicable to many kinds of neurons and other excitable cells [5–7].

Action potentials occur in excitable cells, including neurons, muscle cells, and endocrine cells. It is known that the brain is a complicated network of a tremendous number of neurons. A neuron has a very special shape, which is much different from usual sphere-shaped or disk-like cells [8]. A soma is the main body of neuron from which a long cable called an axon is extended. Neurons transmit and exchange electric signals called action potentials or spikes. Then, the electric signals or information is transmitted in the direction from a dendrite to an axon, see Figure 1 [8].

Figure 1.

Neurons cells.

Doi et al. [8] have briefly explained the framework of the Hodgkin–Huxley formalism to model the action potential generation of neuron sand of other excitable cells. In 2012, Siciliano [9] has studied different numerical methods to solve the HHM. Six different numerical methods are first introduced and compared using a simple and arbitrary ODE. In 2014, Sweilam and Nagy [6] have presented numerical method for solving fractional Hodgkin–Huxley model (FHHM), described as follows:

$C{D}^{\eta }V=I_{ext}-{{\displaystyle \overline{\mathcal{G}}}}_{Na}{m}^{3}h\left(V-E_{Na}\right)-{{\displaystyle \overline{\mathcal{G}}}}_K{n}^{4}h\left(V-E_K\right)-{\mathcal{G}}_L\left(V-E_L\right),D^{\eta }p={\alpha }_p\left(V\right)\left(1-p\right)-{\beta }_p\left(V\right)p,p=m,n,\hslash ,$

where 0 < η ≤ 1.

In this work, the above-described model in [6] will be extended to VOFHHM, which is more general than fractional model.

The main aim of studying the VOFHHM is to show the behavior of the action potential when the derivative expressed as the VOF in order to explain the extent of this impact on the gating variables.

1.2. VOFHHM

It is well known that HHM is based on the parallel thought of a simple circuit with batteries, resistors, and capacitors [9]. Current can be carried through the circuit as ions passing through the membrane (resistors) or by charging the capacitors of the membrane. In this circuit, current flow across the membrane has two major components: one, I_C (μA cm⁻²) associated with charging the membrane capacitance C(μF cm⁻²) and, the other, I_ionic associated with the movement of specific types of ions across the membrane.

For the VOF model, we suggest that:

$I_C=C{D}^{\eta \left(t\right)}V,$E4

where V is the membrane voltage ((mV), measured in volts), 0 < η (t) ≤ 1, and D^η(t) is defined by Grünwald–Letinkov approximation for the VOF derivative as follows:

$D^{\eta \left(t\right)}y(t)={\lim }_{h\to 0}{h}^{-\eta (t)}{\displaystyle {\sum }_{i=0}^{\frac{[t]}{h}}{(-1)}^i\left(\begin{array}{c}\eta (t)\\ i\end{array}\right)y\left(t-ih\right),}$E5

where [t] denotes the integer part of t and h is the step size.

Substituting (4) into (1), we obtain:

$C{D}^{\eta (t)}V=I_{ext}-I_{ionic}(t).$E6

Transitions between permissive and non-permissive states are assumed to obey the relation:

$D^{\eta (t)}p={\alpha }_p(V)(1-p)-{\beta }_p(V)p,p=m,n,\hslash ,$E7

The conductance for each current component can be written in the form [5]:

$G_{Na}={{\displaystyle \overline{\mathcal{G}}}}_{Na}{m}^{3}h,G_{Na}={{\displaystyle \overline{\mathcal{G}}}}_K{n}^4,G_L={\mathcal{G}}_L,$E8

where ${{\displaystyle \overline{g}}}_{Na}$ and ${{\displaystyle \overline{g}}}_K$ are the maximal conductance of sodium and potassium, respectively, and g_L is a constant. Substituting (2) and (8) into (6), we obtain:

$C{D}^{\eta (t)}V=I_{ext}-{{\displaystyle \overline{\mathcal{G}}}}_{Na}{m}^{3}h\left(V-E_{Na}\right)-{{\displaystyle \overline{\mathcal{G}}}}_K{n}^{4}h\left(V-E_K\right)-g_L\left(V-E_L\right).$E9

Eqs. (9) with (7) are a generalization of HHM [5], and it can be described as follows:

$C{D}^{\eta (t)}V=I_{ext}-{{\displaystyle \overline{g}}}_{Na}{m}^{3}h\left(V-E_{Na}\right)-{{\displaystyle \overline{g}}}_K{n}^{4}h\left(V-E_K\right)-g_L\left(V-E_L\right),$

$D^{\eta (t)}p={\alpha }_p\left(V\right)\left(1-p\right)-{\beta }_p\left(V\right)p,p=m,n,\hslash .$

The variable p approaches the steady-state value [5]:

$p^{*}\left(V\right)=\frac{{\alpha }_p}{{\alpha }_p+{\beta }_p},p=m,n,\hslash$E10

when the membrane potential V remains constant.

1.3. Discretizations of the model using NSFD method

The aim of this part is to introduce numerical discretization for Eqs. (7) and (9). The NSFD method [10–17] is used where the step size h in the FDM is replaced by a function ψ(h). Also, the Grünwald–Letinkov definition will be used with V(t_k)=V_k, p(t_k), p(t_k) = p_k, and then, we claim

${\displaystyle {\sum }_{i=0}^{k=1}{c}_i^{\eta \left(t\right)}{p}_{k+1-i}={\alpha }_p{V}_k\left(1-p_k\right)-{\beta }_p{V}_k{p}_k,p=m,n,\hslash ,}$E11

$C{\displaystyle {\displaystyle \sum }_{i=0}^{k+1}}{c}_i^{\eta \left(t\right)}{V}_{k+1-i}=I_m-{{\displaystyle \overline{g}}}_{Na}{m}_{k-1}^2{m}_k\hslash \left(V_k-E_{Na}\right)-{{\displaystyle \overline{g}}}_K{n}_{k-1}^3{n}_k\hslash \left(V_k-E_K\right)-g_L\left(V_k-E_L\right).$E12

Doing some algebraic manipulation to Eqs. (11) and (12) yields the following equations:

${\text{p}}_{\text{k}+1}=\frac{-{{\displaystyle \sum }}_{\text{i}=1}^{\text{k}+1}{\text{c}}_{\text{i}}^{\text{η}\left(\text{t}\right)}{\text{p}}_{\text{k}+1-\text{i}}+{\text{α}}_{\text{p}}{\text{V}}_{\text{k}}\left(1-{\text{p}}_{\text{k}}\right)-{\text{β}}_{\text{p}}{\text{V}}_{\text{k}}{\text{p}}_{\text{k}}}{{\text{c}}_0^{\text{η}(t)}},\text{p}=\text{m},\text{n},\hslash$E13

$V_{k+1}=\frac{I_{ext}-{{\displaystyle \sum }}_{i=1}^{k+1}{c}_i^{\eta \left(t\right)}{V}_{k+1-i}-{{\displaystyle \overline{g}}}_{Na}{m}_{k-1}^2{m}_k\hslash \left(V_k-E_{Na}\right)-{{\displaystyle \overline{g}}}_K{n}_{k-1}^3{n}_k\hslash \left(V_k-E_K\right)-g_L\left(V_k-E_L\right)}{C{c}_0^{\eta (t)}}$E14

where $c_0^{\eta (t)}=\psi {(h)}^{-\eta (t)},\psi \left(h\right)=e^h-1$.

1.4. Numerical experiments

VOFHHM can be summarized neatly into four separate variable-order fractional ODEs with some supporting functions. These equations are described in Eqs. (7) and (9), where the rate functions are listed below [5]:

${\alpha }_n\left(\text{V}\right)=0.01\left(\text{V}+1\right){\left[exp\left(\frac{V+10}{10}\right)\right]}^{-1},{\beta }_n\left(V\right)=0.125exp\left(\frac{V}{80}\right),$

${\alpha }_m\left(\text{V}\right)=0.1\left(\text{V}+25\right){\left[exp\left(\frac{V+25}{10}\right)\right]}^{-1},{\beta }_m\left(V\right)=4exp\left(\frac{V}{18}\right),$

${\alpha }_{\hslash }\left(\text{V}\right)=0.07exp\left(\frac{V}{20}\right),{\beta }_{\hslash }\left(V\right)={\left[exp\left(\frac{V+30}{10}\right)\right]}^{-1},$

and

${{\displaystyle \overline{g}}}_{\text{Na}}=120\text{mS}/{\text{cm}}^2,{{\displaystyle \overline{g}}}_{\text{K}}=36\text{mS}/{\text{cm}}^2,g_L=0.3\text{mS}/{\text{cm}}^2,E_{Na}=120\text{mV},E_K=-12\text{mV},E_L=10.6\text{mV},$ with the initial conditions V(0) = V₀, and p(0) = p₀ = p^* (V₀). The behavior of the neuron can be simulated for different initial values of V. Figure 2 shows an action potential simulated by the solver ode45, i.e., fourth-order Runge–Kutta method, and NSFD method at η(t) = 1, of the variable-order Hodgkin–Huxley equations for zero I_ext, and V₀ = −40 mV. Also, in this figure, all parts of the action potential are presented including the rapid upraise, downfall, and unexcitable phase, according to the number sequence as follows:

1. Number 1 means that the action potential begins to fire, the sodium channels are now open, and an influx of sodium occurs.

2. Number 2 means that the membrane becomes more positive at this point due to sodium entering the cell. Now, the potassium channels open and leave the cell.

3. Number 3 means that the sodium channels have now become refractory and no more sodium is allowed to enter.

4. Number 4 means that the sodium channel is still in its refractory stage and only the potassium ions are passing through the membrane. As the potassium ions continue to leave the cell, the membrane potential moves toward the resting value.

5. Number 5 means that the potassium channels close and the sodium channels begin to leave the refractory phase and reset to its resting phase.

Figure 2.

An action potential created by using ode45 and NSFD method at η (t) = 1.

Finally, the potassium returns to its resting value where it can wait another action potential (for more details see [8, 9]). Moreover, in Figure 3, the gating variables are correctly constrained between 0 and 1 as expected using ode45 and NSFD method, respectively.

Figure 3.

Plot of the different gating variables using ode45 (left) and NSFD method (right) at η (t) = 1.

Figure 4.

The approximated action potential (left) and the gating variables (right) using the NSFD method at (t) = 0.99 – (0.01/100)t.

Figures 4–8 show the approximated action potential and the gating variables using the NSFD method obtained for different values of (t), where η(t) = 0.99 − (0.01/100)t, η(t) = 0.87 − (0.01/100)t, η(t) =0.79 − (0.001)t, and η(t) = 0.67 − (0.001)t, respectively. These figures show an example of a numerically solved solution of the Hodgkin–Huxley equations. Figures 4–8 (left) are a waveform of the membrane potential. A pulsatile input applied at a time t = 5 ms ((ms) time measured in second) induces an action potential. Figures 4–7 (right) show the dynamics of all three gating variables m, n, and h. Sodium activation, m, changes much more rapidly than either h or n. Figure 9 (left) describes the relationship between the variables n and h, where n is the activation variable of potassium channel expressed in a dark color, while h is inactivation variable of sodium channel expressed in a light color. The action potential, n, be in stillness stage and begins to enter the cell until it reaches the highest activity, unlike the action potential, h, at the top of its activity. Step by step, the action potential, h, decreases its activities and continues to leave the cell and moves toward the resting value. At a certain moment in time, the intersection between them are occurred, where n increased while h decreased. Figure 9 (right) describes the relationship between the variables m and h, where m is the activation variable of potassium channel expressed in a dark color, while h is inactivation variable of sodium channel expressed in a light color. The dynamics variables m and h have the same dynamics similar to n and h, where m increased while h decreased. Figure 10 describes the dynamics of gating variables n and m in 3D. The variables n and m are the activation variables of K⁺ and Na⁺ ionic channels, respectively. They are expressed in a dark color and a light color, respectively. They are in stillness stage and begin to enter the cell until they reach the highest activity, where they are in the case of increasing.

Figure 5.

The approximated action potential (left) and the gating variables (right) using the NSFD method at (t) = 0.87 – (0.01/100)t.

Figure 6.

The approximated action potential (left), and the gating variables (right) using the NSFD method at (t) = 0.79 – (0.001)t.

Figure 7.

The approximated action potential (left) and the gating variables (right) using the NSFD method at (t) = 0.67 – (0.001)t.

Figure 8.

An action potential created from the NSFD method using different values of η (t).

Figure 9.

The waveforms of the gating variables n and h (left), and m and h (right) in 3D.

Figure 10.

The waveforms of the gating variables n and m in 3D.

2. Real-life models governed by nonlinear delay differential equations (DDEs)

2.1. Introduction

DDEs are differential equations in which the derivatives of some unknown functions at present time are dependent on the values of the functions at previous times. In real-world systems, delay is very often encountered in many practical systems, such as control systems [18], lasers, traffic models [19], metal cutting, epidemiology, neuroscience, population dynamics [20], and chemical kinetics [21]. Recent theoretical and computational advancements in DDEs reveal that DDEs are capable of generating rich and plausible dynamics with realistic parameter values. Naturally, occurrence of complex dynamics is often generated by well-formulated DDE models. This is simply due to the fact that a DDE operates on an infinite-dimensional space consisting of continuous functions that accommodate high-dimensional dynamics.

For example, the Lotka–Volterra predator prey model [22] with crowding effect does not produce sustainable oscillatory solutions that describe population cycles. However, the Nicholson’s blowflies model [23] can generate rich and complex dynamics. Delayed fractional differential equations (DFDEs) are correspondingly used to describe dynamical systems [23]. In recent years, DFDEs begin to arouse the attention of many researchers [7, 19, 20, 25–27]. Simulating these equations is an important technique in the research, and accordingly, finding effective numerical methods for the DFDEs is a necessary process. Several methods based on Caputo or Riemann–Liouville definitions [28] have been proposed and analyzed. For instance, based on the predictor–corrector scheme, Diethelm et al. introduced the ABM algorithm [29–31] and mean while some error analysis was presented to improve the numerical accuracy. In 2011, Bhalekar and Daftardar-Gejji [25] have been extended the ABM algorithm to solve DDEs of fractional order and presented numerical illustrations to demonstrate utility of the method. In 2012, Ma et al. [32] have presented the numerical solution of a VOF financial system which is calculated by using the ABM method.

This section presents class of numerical method for solving the variable-order fractional nonlinear delay differential equations (VOFDDEs). The main aim of this part is to study VOFDDEs numerically.

2.2. ABM method for the VOFDDEs

In this section, the ABM algorithm has been extended to study the VOFDDEs, where the derivative is defined in the Caputo VOF sense. Special attention is given to prove the error estimate of the proposed method. Numerical test examples are presented to demonstrate utility of the method. Chaotic behaviors are observed in variable-order one-dimensional delayed systems.

In the following, we apply the ABM predictor–corrector method to implement the numerical solution of VOFDDEs.

Let us consider the following VOF system:

$D_t^{\alpha \left(t\right)}y\left(t\right)=f\left(t,y\left(t\right),y\left(t-\tau \right)\right),t\in \left[o,T\right],0<\alpha (t)\le 1.$E15

$y\left(t\right)=g\left(t\right),t\in \left[-\tau ,0\right],$E16

where f is in general a nonlinear function.

Also, consider a uniform grid {t_n = nh: n = −k, −k + 1,…, −1,0,1,…, N} where k and N are integers

such that h = τ/k. Let

$y_h\left(t_j\right)=g\left(t_j\right),j=-k,-k+1,\dots ,-1,0,$E17

and note that

$y_h\left(t_j-\tau \right)=y_h\left(jh-kh\right)=y_h\left(t_{j-k}\right),j=0,1,2,\dots ,N.$E18

Applying the integral $J_{t_{n+1}}^{\alpha (t_{n+1})}$, which is defined by the Riemann–Liouville variable-order integral as

$J_x^{\alpha \left(x\right)}f\left(x\right)=\frac{1}{\text{Γ}(\alpha \left(x\right))}{\int }_a^x{(x-t)}^{\alpha \left(x\right)-1}f\left(t\right)dt,$E19

on both sides of (15) and using (16), we claim to:

$y\left(t_{n+1}\right)=g\left(0\right)+\frac{1}{\text{Γ}\left(\alpha \left(t_{n+1}\right)\right)}{\int }_0^{t_{n+1}}{\left(t_{n+1}-\zeta \right)}^{\alpha \left(t_{n+1}\right)-1}f\left(\zeta ,y\left(\zeta \right),y\left(\zeta -\tau \right)\right)d\zeta .$E20

Further, the integral in Eq. (20) is evaluated using product trapezoidal quadrature formula [29–31, 33]. Then, we have the following corrector formula:

$y_h\left(t_{n+1}\right)=g\left(0\right)+\frac{h^{\alpha \left(t_{n+1}\right)}}{\text{Γ}\left(\alpha \left(t_{n+1}\right)+2\right)}f\left(t_{n+1},y_h\left(t_{n+1}\right),y_h\left(t_{n+1}-\tau \right)\right)+\frac{h^{\alpha \left(t_{n+1}\right)}}{\text{Γ}\left(\alpha \left(t_{n+1}\right)+2\right)}{\displaystyle {\displaystyle \sum }_{}^{}}j=0n{a}_{j,n+1}f(t_j,y_h\left(t_j\right),y_h\left(t_j-\tau \right),$

or

$y_h\left(t_{n+1}\right)=g\left(0\right)+\frac{h^{\alpha \left(t_{n+1}\right)}}{\text{Γ}\left(\alpha \left(t_{n+1}\right)+2\right)}f\left(t_{n+1},y_h\left(t_{n+1}\right),y_h\left(t_{n+1-k}\right)\right)+\frac{h^{\alpha \left(t_{n+1}\right)}}{\text{Γ}\left(\alpha \left(t_{n+1}\right)+2\right)}{\displaystyle {\sum }_{j=0}^n{a}_{j,n+1}f(t_j,y_h\left(t_j\right),y_h\left(t_{j-k}\right)}$E21

where

$a_{j,n+1}=\left\{\begin{array}{c}{n}^{\alpha \left(t_{n+1}\right)+1}-\left(n-\alpha \left(t_{n+1}\right)\right){\left(n+1\right)}^{\alpha \left(t_{n+1}\right)},j=0,\\ {(n-j+2)}^{\alpha \left(t_{n+1}\right)+1}+{(n-j)}^{\alpha \left(t_{n+1}\right)+1}-2{\left(n-j+1\right)}^{\alpha \left(t_{n+1}\right)+1},1\le j\le n,\\ 1,j=n+1,\end{array}\right\}$E22

$y_h\left(t_{j-k}\right)\approx v_{n+1}=\left\{\begin{array}{c}{\delta }_{y_{n-m+1}}+{(1-\delta )}_{y_{n-m+1}},ifm>1,\\ {\delta }_{y_{n+1}^p}+{(1-\delta )}_{y_n},ifm=1,\end{array}\right\}$E23

0 ≤ δ < 1 and the unknown term y_h(t_n+1) appears on both sides of (21).

Due to nonlinearity off, Eq. (21) cannot be solved explicitly for y_h(t_n+1), so we replace the term y_h(t_n+1) on the right-hand side by an approximation $y_h^p\left(t_{n+1}\right),$ which is called predictor. The product rectangle rule is used in (20) to evaluate the predictor term

$y_h^p\left(t_{n+1}\right)=g\left(0\right)+\frac{1}{\text{Γ}\left(\alpha \left(t_{n+1}\right)\right)}{\displaystyle {\displaystyle \sum }_{j=0}^n}{b}_{j,n+1}f(t_j,y_h\left(t_j\right),y_h\left(t_j-\tau \right),$

or

$y_h^p\left(t_{n+1}\right)=g\left(0\right)+\frac{1}{\text{Γ}\left(\alpha \left(t_{n+1}\right)\right)}{\displaystyle {\displaystyle \sum }_{j=0}^n}{b}_{j,n+1}f(t_j,y_h\left(t_j\right),y_h\left(t_{j-k}\right)$E24

where $b_{j,n+1}=\frac{h^{\alpha (t_{n+1})}}{\alpha (t_{n+1})}({(n-j+1)}^{\alpha \left(t_{n+1}\right)}+{(n-j)}^{\alpha \left(t_{n+1}\right)}.$

2.3. Error analysis of the algorithm

In this subsection, we aim to introduce the following lemma, which will be used in the proof of main theorem.

Lemma 2.1. [31] (a) Let z ∈ C¹[0, T], then

$\left|{\displaystyle \underset{0}{\overset{t_{n+1}}{{\displaystyle \int }}}}{\left(t_{n+1}-t\right)}^{\alpha \left(t_{n+1}\right)-1}z\left(t\right)dt-{\displaystyle {\displaystyle \sum }_{j=0}^n}{b}_{j,n+1}z(t_j)\right|\le \frac{1}{\alpha \left(t_{n+1}\right)}{\Vert z^{\prime }\Vert }_{\infty }{t}_{n+1}^{\alpha \left(t_{n+1}\right)}h.$E25

(b) Let z ∈ C²[0, T], then

$\left|{\displaystyle \underset{0}{\overset{t_{n+1}}{{\displaystyle \int }}}}{\left(t_{n+1}-t\right)}^{\alpha \left(t_{n+1}\right)-1}z\left(t\right)dt-{\displaystyle {\displaystyle \sum }_{j=0}^{n+1}}{a}_{j,n+1}z(t_j)\right|\le C_{\alpha \left(t_{n+1}\right)}^{T_r}{z}^{\text{'}\text{'}}{}_{\infty }{t}_{n+1}^{\alpha \left(t_{n+1}\right)}h.$E26

Proof. To prove statement (a), by construction of the product rectangle formula, we find that the quadrature error has the representation [31]

${\displaystyle \underset{0}{\overset{t_{n+1}}{{\displaystyle \int }}}}{\left(t_{n+1}-t\right)}^{\alpha \left(t_{n+1}\right)-1}z\left(t\right)dt-{\displaystyle {\displaystyle \sum }_{j=0}^n}{b}_{j,n+1}z(t_j)={\displaystyle {\displaystyle \sum }_{j=0}^n}{\displaystyle \underset{jh}{\overset{\left(j+1\right)h}{{\displaystyle \int }}}}{\left(t_{n+1}-t\right)}^{\alpha \left(t_{n+1}\right)-1}\left(z\left(t\right)-z\left(t_j\right)\right)dt.$

Apply the mean value theorem of differential calculus to the second factor of the integrand on the right-hand side of the above equation

$\left|{\displaystyle \underset{0}{\overset{t_{n+1}}{{\displaystyle \int }}}}{\left(t_{n+1}-t\right)}^{\alpha \left(t_{n+1}\right)\text{-}1}z\left(t\right)dt-{\displaystyle {\displaystyle \sum }_{j=0}^n}{b}_{j,n+1}z\left(t_j\right)\right|$

$\le {\Vert z^{\prime }\Vert }_{\infty }{\displaystyle {\displaystyle \sum }_{j=0}^n}{\displaystyle \underset{jh}{\overset{\left(j+1\right)h}{{\displaystyle \int }}}}{\left(t_{n+1}-t\right)}^{\alpha \left(t_{n+1}\right)-1}\left(t-jh\right)dt,$

$={\Vert z^{\prime }\Vert }_{\infty }\frac{h^{1+\alpha \left(t_{n+1}\right)}}{\alpha \left(t_{n+1}\right)}{\displaystyle {\displaystyle \sum }_{j=0}^n}(\frac{1}{1+\alpha \left(t_{n+1}\right)}[{\left(n+1-j\right)}^{1+\alpha \left(t_{n+1}\right)}$

$-{\left(n-j\right)}^{1+\alpha \left(t_{n+1}\right)}]-{\left(n-j\right)}^{\alpha \left(t_{n+1}\right)},$

$={\Vert z^{\prime }\Vert }_{\infty }\frac{h^{1+\alpha \left(t_{n+1}\right)}}{\alpha \left(t_{n+1}\right)}(\frac{{\left(n+1\right)}^{1+\alpha \left(t_{n+1}\right)}}{1+\alpha \left(t_{n+1}\right)}-{\displaystyle {\displaystyle {\displaystyle \sum }_{j=0}^n}{j}^{\alpha \left(t_{n+1}\right)}),}$

$={\Vert z^{\prime }\Vert }_{\infty }\frac{h^{1+\alpha \left(t_{n+1}\right)}}{\alpha \left(t_{n+1}\right)}({\displaystyle \underset{o}{\overset{t_{n+1}}{{\displaystyle \int }}}}{j}^{\alpha \left(t_{n+1}\right)}dt-{\displaystyle {\displaystyle \sum }_{j=0}^n}{j}^{\alpha \left(t_{n+1}\right)}).$

Here, the term in parentheses is simply the remainder of the standard rectangle quadrature formula, applied to the function $t^{\alpha \left(t_{n+1}\right)}$ and taken over the interval [0, n + 1]. Since the integrand is monotonic, we may apply some standard results from quadrature theory to find that this term is bounded by the total variation of the integrand, thus

$\left|{\displaystyle \underset{0}{\overset{t_{n+1}}{{\displaystyle \int }}}}{\left(t_{n+1}-t\right)}^{\alpha \left(t_{n+1}\right)-1}z\left(t\right)dt-{\displaystyle {\displaystyle \sum }_{j=0}^n}{b}_{j,n+1}z\left(t_j\right)\right|\le \frac{1}{\alpha \left(t_{n+1}\right)}{\Vert z^{\prime }\Vert }_{\infty }{\left(n+1\right)}^{\alpha \left(t_{n+1}\right)h^{1+\alpha \left(t_{n+1}\right)}.}$

Similarly, we can prove (b).

Now, let us consider that f(⋅) in (15) satisfies the following Lipschitz conditions with respect to its variables

$\left|f\left(t,y_1,u\right)-f\left(t,y_2,u\right)\right|\le L_1\left|y_1-y_2\right|,$

$\left|f\left(t,y,u_1\right)-f\left(t,y,u_2\right)\right|\le L_2\left|u_1-u_2\right|,$E27

where L₁ and L₂ are positive constants.

Theorem 1. Suppose the solution y(t) ∈ C² [0, T], of the Eqs. (15) and (16), satisfies the following two conditions:

$\left|{\displaystyle \underset{0}{\overset{t_{n+1}}{{\displaystyle \int }}}}{\left(t_{n+1}-t\right)}^{\alpha \left(t_{n+1}\right)-1}{D}_t^{\alpha \left(t_{n+1}\right)}y\left(t\right)dt-{\displaystyle {\displaystyle \sum }_{j=0}^n}{b}_{j,n+1}{D}_t^{\alpha \left(t_{n+1}\right)}y\left(t_j\right)\right|\le C{t}_{n+1}^{{\gamma }_1}{h}^{{\delta }_1}$E28

$\left|{\displaystyle \underset{0}{\overset{t_{n+1}}{{\displaystyle \int }}}}{\left(t_{n+1}-t\right)}^{\alpha \left(t_{n+1}\right)-1}{D}_t^{\alpha \left(t_{n+1}\right)}y\left(t\right)dt-{\displaystyle {\displaystyle \sum }_{j=0}^n}{a}_{j,n+1}{D}_t^{\alpha \left(t_{n+1}\right)}y\left(t_j\right)\right|\le C{t}_{n+1}^{{\gamma }_2}{h}^{{\delta }_2},$

with some γ₁, γ₂ ≥ 0, and δ₁, δ₂ < 0, then for some suitable T > 0, we have

${\displaystyle \underset{0\le j\le N}{\max }}\left|y(t_j)-y_h(t_j)\right|=\text{O}(h^q),$

where q = min (δ₁+ α (t), δ₂), N = [T/h], and C is a positive constant.

Proof. We will use the mathematical induction to prove the result. Suppose that the conclusion is true for j = 0, 1,…, n, then we have to prove that the inequality also holds for j = n + 1. To do this, we first consider the error of the predictor $y_{n+1}^p$. From (27), we have

$\left|f\left(t_j,y(t_j),y(t_j-\tau )\right)-f\left(t_j,y_j,{\upsilon }_j\right)\right|\le \left|f\left(t_j,y(t_j),y(t_j-\tau )\right)-f\left(t_j,y_j,y(t_j-\tau )\right)\right|+\left|f\left(t_j,y_j,y(t_j-\tau )\right)-f\left(t_j,y_j,{\upsilon }_j\right)\right|\le L_1{h}^q+L_2{h}^q=\left(L_1+L_2\right)h^q.$

So

$\left|y(t_{n+1})-y_{n+1}^p\right|=\frac{1}{\text{Γ}(\alpha \left(t_{n+1}\right))}\left|{\displaystyle \underset{0}{\overset{t_{n+1}}{{\displaystyle \int }}}}{\left(t_{n+1}-t\right)}^{\alpha \left(t_{n+1}\right)-1}f\left(t,y(t),y(t-\tau )\right)dt-{\displaystyle {\displaystyle \sum }_{j=0}^n}{b}_{j,n+1}f\left(t_j,y_j,{\upsilon }_j\right)\right|$

$\begin{array}{l}\le \frac{1}{\text{Γ}(\alpha \left(t_{n+1}\right))}\left|{\displaystyle \underset{0}{\overset{t_{n+1}}{{\displaystyle \int }}}}{\left(t_{n+1}-t\right)}^{\alpha \left(t_{n+1}\right)-1}{D}_t^{\alpha \left(t_{n+1}\right)}y\left(t\right)dt-{\displaystyle {\displaystyle \sum }_{j=0}^n}{b}_{j,n+1}{D}_t^{\alpha \left(t_{n+1}\right)}y\left(t_j\right)\right|+\\ \frac{1}{\text{Γ}(\alpha \left(t_{n+1}\right))}{\displaystyle {\displaystyle \sum }_{j=0}^n}{b}_{j,n+1}\left|f\left(t_j,y(t_j),y(t_j-\tau )\right)-f\left(t_j,y_j,y(t_j-\tau )\right)\right|+\\ \left|f\left(t_j,y_j,y(t_j-\tau )\right)-f\left(t_j,y_j,{\upsilon }_j\right)\right|\\ \le \frac{C{t}_{n+1}^{{\gamma }_1}}{\text{Γ}(\alpha \left(t_{n+1}\right))}+\frac{1}{\text{Γ}(\alpha \left(t_{n+1}\right))}{\displaystyle {\displaystyle \sum }_{j=0}^n}{b}_{j,n+1}\left(L_1+L_2\right)C{h}^q,\end{array}$

and from

${\displaystyle {\displaystyle \sum }_{j=0}^n}{b}_{j,n+1}={\displaystyle {\displaystyle \sum }_{j=0}^n}{\displaystyle \underset{\text{jh}}{\overset{\left(\text{j}+1\right)\text{h}}{{\displaystyle \int }}}}{\left(t_{n+1}-t\right)}^{\alpha \left(t_{n+1}\right)-1}dt=\frac{1}{\alpha \left(t_{n+1}\right)}{t}_{n+1}^{\alpha \left(t_{n+1}\right)}\le \frac{1}{\alpha \left(t_{n+1}\right)}{T}^{\alpha \left(t_{n+1}\right)},$

we have

$\left|y(t_{n+1}-y_{n+1}^p)\right|\le \frac{C{T}^{{\gamma }_1}}{\text{Γ}\left(\alpha \left(t_{n+1}\right)\right)}{h}^{{\delta }_1}+\frac{C{T}^{\alpha \left(t_{n+1}\right)}}{\text{Γ}\left(\alpha \left(t_{n+1}\right)+1\right)}{h}^q.$

Since

${\displaystyle {\displaystyle \sum }_{j=0}^n}{a}_{j,n+1}\le {\displaystyle {\displaystyle \sum }_{j=0}^n}\frac{h^{\alpha \left(t_{n+1}\right)}}{\alpha \left(t_{n+1}\right)(\alpha \left(t_{n+1}\right)+1)}[{\left(k-j+2\right)}^{\alpha \left(t_{n+1}\right)+1}-{\left(k-j+1\right)}^{\alpha \left(t_{n+1}\right)+1}-{(k-j+1)}^{\alpha \left(t_{n+1}\right)+1}+{(k-j)}^{\alpha \left(t_{n+1}\right)+1}],$

$\begin{array}{l}={\displaystyle \underset{0}{\overset{t_{n+1}}{{\displaystyle \int }}}}{\left(t_{n+2}-t\right)}^{\alpha \left(t_{n+1}\right)}dt-{\displaystyle \underset{0}{\overset{t_{n+1}}{{\displaystyle \int }}}}{\left(t_{n+1}-t\right)}^{\alpha \left(t_{n+1}\right)}dt,\\ =\frac{1}{\alpha \left(t_{n+1}\right)}{\displaystyle \underset{0}{\overset{t_{n+1}}{{\displaystyle \int }}}}\left[{\left(t_{n+1}-t\right)}^{\alpha \left(t_{n+1}\right)}\right](t)\prime dt,\\ =\frac{1}{\alpha \left(t_{n+1}\right)}{t}_{n+1}^{\alpha \left(t_{n+1}\right)}\le T^{\alpha \left(t_{n+1}\right)},\end{array}$

we have

$\left|y(t_{n+1})-y_{n+1}^p\right|=\frac{1}{\text{Γ}(\alpha \left(t_{n+1}\right))}\left|{\displaystyle \underset{0}{\overset{t_{n+1}}{{\displaystyle \int }}}}{\left(t_{n+1}-t\right)}^{\alpha \left(t_{n+1}\right)-1}f\left(t,y(t),y(t-\tau )\right)dt-{\displaystyle {\displaystyle \sum }_{j=0}^n}{a}_{j,n+1}f\left(t_j,y_j,{\upsilon }_j\right)-a_{n+1,n+1}f\left(t_{n+1},y_{n+1}^p,{\upsilon }_{n+1}\right)\right|,$

$\begin{array}{l}\le \frac{1}{\text{Γ}\left(\alpha \left(t_{n+1}\right)\right)}[\left|\begin{array}{l}{\displaystyle \underset{0}{\overset{t_{n+1}}{{\displaystyle \int }}}}{\left(t_{n+1}-t\right)}^{\alpha \left(t_{n+1}\right)-1}f\left(t,y\left(t\right),y\left(t-\tau \right)\right)dt\\ -{\displaystyle {\displaystyle \sum }_{j=0}^n}{a}_{j,n+1}f\left(t_j,y_j,y\left(t_j-\tau \right)\right)\end{array}\right|+\\ {\displaystyle {\displaystyle \sum }_{j=0}^n}{a}_{j,n+1}\left|f\left(t_j,y_j,y\left(t_j-\tau \right)\right)-f\left(t_j,y_j,{\upsilon }_j\right)\right|+\\ a_{n+1,n+1}\left|f\left(t_{n+1},y(t_{n+1}),y(t_{n+1}-\tau )\right)-f\left(t_{n+1},y_{n+1}^p,{\upsilon }_{n+1}\right)\right|]\\ \le \frac{1}{\text{Γ}\left(\alpha \left(t_{n+1}\right)\right)}[C{t}_{n+1}^{{\gamma }_2}{h}^{{\delta }_2}+C{h}^q{\displaystyle {\displaystyle \sum }_{j=0}^n}{a}_{j,n+1}+\\ a_{n+1,n+1}L\left(h^q+\frac{C{T}^{{\gamma }_1}}{\text{Γ}\left(\alpha \left(t_{n+1}\right)\right)}{h}^{{\delta }_1}+\frac{C{T}^{\alpha \left(t_{n+1}\right)}}{\text{Γ}\left(\alpha \left(t_{n+1}\right)+1\right)}{h}^q\right)]\end{array}$

$\le \frac{1}{\text{Γ}\left(\alpha \left(t_{n+1}\right)\right)}[C{T}^{{\gamma }_2}+\frac{L{h}^h}{\alpha \left(t_{n+1}\right)\left(\alpha \left(t_{n+1}\right)+1\right)}+\frac{CL{T}^{{\gamma }_1}}{\text{Γ}\left(\alpha \left(t_{n+1}\right)+2\right)}+\frac{CL{T}^{\alpha \left(t_{n+1}\right)}{h}^{\alpha \left(t_{n+1}\right)}}{\alpha \left(t_{n+1}\right)\text{Γ}\left(\alpha \left(t_{n+1}\right)+2\right)}]h^q=C{h}^q.$

2.4. Numerical examples

We present the purpose of this subsection to show that the proposed scheme designed provides good approximations for VOFDDEs. Throughout this subsection, we discuss four examples and their numerical solutions.

Example 2.1. Consider a VOFDDE:

$D_t^{\alpha \left(t\right)}y\left(t\right)=\frac{2y(t-2)}{1+y{(t-2)}^{9.65}}-y\left(t\right),$

$y\left(t\right)=0.5,t\le 0$E29

Figure 11.

The numerical behavior of system (29) and chaotic attractors at α = 1.

Figure 12.

The numerical behavior of system (29) and chaotic attractors at α(t) = 0.97.

Figure 13.

The numerical behavior of system (32) and chaotic attractors at α(t) = 0.99 – (0.01/100)t.

In Figures 11 and 1(a, b) show the solutions y(t) and y(t − 2) of the system (29), for α = 1 and h = 0.1, whereas Figure 1(c) shows phase portrait of the system, i.e., plot of y(t) versus y(t − 2) for the same value of α. In this figure, it observed that the system (29) shows aperiodic chaotic behavior. Moreover, Figure 12 shows the plot of the numerical solutions y(t) and the plot of y(t) versus y(t − 2), respectively, at α(t) = 0.97. In the following, we choose different cases for α(t). Figures 13 and 14 show the solutions y(t) and y(t − 2) of the given system for a(t) = 0.99 − (0.01/100)t and α(t) = 0.95 − (0.01/100)t, respectively. The chaotic portrait for these values of α(t) is shown. In Figure 15, we have decreased the value of α(t), where α(t) = 0.85 − (0.01/100)t, and observed that the system becomes periodic.

Figure 14.

The numerical behavior of system (29) and chaotic attractors at α(t) = 0.95 – (0.01/100)t.

Figure 15.

The numerical behavior of system (29) and chaotic attractors at α(t) = 0.85 – (0.01/100)t.

2.5. The four-year life cycle of a population of lemmings model

Lemmings are small rodents, usually found in or near the Arctic, in tundra biomes. It used to be that every three to 4 years, there were massive numbers of lemmings in the mountains and then the next year it was gone. Therefore, it has been an interesting thing to try to find out why. Are there factors that affect the lemming population such as climate, temperature, precipitation, and the like, that is, what we refer to as the lemming cycle. The modern scientific study of lemmings started with work carried out by the Norwegian Professor of Zoology Robert Collett, who at the end of the nineteenth century gathered a great deal of information about lemmings. To try to understand why lemmings fluctuate both regularly and extensively is indeed an important problem in ecology [34]. In [27], Tavernini has solved a model of the 4-year life cycle of a population of lemmings in an integer order. In [25], Bhalekar and Daftardar-Gejji have studied the model in a fractional order. In the following example [27], we study the extension of the lemming model:

Example 2.2. Consider the variable-order version of the four-year life cycle of a population of lemmings

$D_t^{\alpha \left(t\right)}y\left(t\right)=3.5y\left(t\right)\left(1-\frac{y\left(t-0.74\right)}{19}\right),y\left(0\right)-19.00001,y\left(t\right)=19,t<0.$E30

Figure 16.

The numerical behavior of system (33) and plots of the (t), y(t – 0.74) relation for α = 1.

Figure 17.

The numerical behavior of system (33), and plots of the y(t), y(t – 0.74) relation for α(t) = 0.97.

Figure 18.

The numerical behavior of system (33) and the stretching phenomena for α(t) = 0.99 – (0.01/100)t.

Figure 16 shows the solutions y(t) and y(t − 0.74) of system (33) for α = 1 and h = 0.1 and shows phase portrait of y(t) versus y(t − 0.74) for the same value of α. Figure 17 shows the plot of the numerical solutions y(t) and the plot of y(t) versus y(t − 0.74), respectively, at (t) = 0.97. The numerical results for VOFDDEs at different values of α(t) are given in Figures 18–22. These figures show the stretching phenomena between the numbers of lemmings in y(t) versus y(t − 2), for the values α(t) = 0.99 − (0.01/100)_t, α(t) = 0.98 − (0.02/100)_t; α(t) = 0.95 − (0.01/100)_t, α(t) = 0.87 − (0.02/100)_t, and α(t) = 0.85 − (0.01/100)_t, respectively. It is observed that the phase portrait gets stretched as the value of α(t) decreases, and this stretching is towards positive side of the axis.

Figure 19.

The stretching.

Figure 20.

The stretching phenomena for α(t) = 0.98 − (0.02/100)t. Phenomena for α(t) = 0.95 − (0.01/100)t.

Figure 21.

The stretching phenomena for α(t) = 0.87 − (0.02/100)t.

Figure 22.

The stretching phenomena for α(t) = 0.85 − (0.01/100)t.

2.6. The enzyme kinetics with an inhibitor molecule model

An enzyme inhibitor is a molecule that binds to an enzyme and decreases its activity or completely inhibits the enzyme catalytic activity. It is well known that all these inhibitors follow the same rule to interplay in enzyme reaction. Furthermore, there are many factors that affect enzyme’s activity, such as temperature and pH. Many drug molecules are enzyme inhibitors, so their discovery and improvement are an active area of research in biochemistry and pharmacology to protect enzyme from any change. Therefore, studying the enzyme kinetics and structure–function relationship is vital to understand the kinetics of enzyme inhibition that in turn is fundamental to the modern design of pharmaceuticals in industries [35]. The frequency conversion mechanism in enzymatic feedback systems has been investigated with computer simulations in 1984 by Okamoto and Hayashi [36]. In [25], Bhalekar and Daftardar-Gejji have been study the enzyme kinetics with an inhibitor molecule model in a fractional order. In the following, we will study the model in [25] in general form where the derivative is given in a VOF.

Example 2.3. Consider the variable-order version of four-dimensional enzyme kinetics with an inhibitor molecule

$D_t^{\alpha \left(t\right)}{y}_1\left(t\right)=10.5-\frac{y_1(t)}{1+0.0005y_4^3(t-4)},$

$D_t^{\alpha \left(t\right)}{y}_2\left(t\right)=\frac{y_1(t)}{1+0.0005y_4^3(t-4)}-y_2\left(t\right),$

$D_t^{\alpha \left(t\right)}{y}_3\left(t\right)=y_2\left(t\right)-y_3\left(t\right),$

$D_t^{\alpha \left(t\right)}{y}_4\left(t\right)=y_3\left(t\right)-0.5y_4\left(t\right),$

$y(t)={[60,10,10,20]}^T,t\le 0.$E31

Figure 23.

The numerical behavior of system (34) at α = 1 (left) and α = 95 (right).

Figure 23 shows the solutions y(t), (1≤ i ≤ 4), with step size h= 1, for = 1, and α = 95, respectively. Whereas Figures 24–26 show the numerical results in case of variable order at (t) = 0.99 − (0.01/100)t, α(t) = 0.95 − (0.01/100)t, α(t) = 0.91 − (0.01/100)t, and α(t) = 0.83 − (0.01/100)t, respectively. For 0.90 < α(t) ≤ 1, the height of oscillations of the solutions increases as t increases as shown in Figures 24 and 25, while in Figure 26, the system settles down for sufficiently large t, for α(t) < 90.

Figure 24.

The numerical behavior of system (34) at α(t) = 0.99 – (0.01/100)t (left) and α(t) = 0.95 – (0.01/100)t (right).

Figure 25.

The numerical behavior of system (34) at α(t) = 0.91 – (0.01/100)t.

Figure 26.

The numerical behavior of system (34) at α(t) = 0.83 – (0.01/100)t.

2.7. The Chen system model

The atmosphere is a layer of gases surrounding the planet Earth that has, at each altitude above each point of the Earth’s surface, a density, a pressure, a temperature, etc., and all these vary over time. It is of course unthinkable to know all this infinite number of data and to understand something about it. Approximations must be made. Edward Norton Lorenz (1917–2008) is a pioneer of chaos theory. He introduced the strange attractor notion and coined the term butterfly effect. In 1963 [37], he developed a simplified mathematical model for atmospheric convection when he studied the atmospheric convection. His model of the atmosphere was reduced to only three numbers x, y, and z and the evolution of the atmosphere to a tiny equation, where each point x, y, and z in space symbolizes a state of the atmosphere and the evolution follows a vector field. Later, several dynamical systems exhibiting chaos have been presented in various branches of science [24]. For example, in 1999, Chen and Ueta found another simple three-dimensional autonomous system, which is not topologically equivalent to Lorenz’s system and which has a chaotic attractor [38]. In [39], Chen proved that “The Chen system is a special case of the Lorenz system.” In 2004, Li and Peng [40] have presented Chen system with a fractional order. In 2012, Bhalekar et al. [41] proposed the fractional-order Chen system with time delay. In the following example, we will extend Chen system to a VOF, which can be more general.

Example 2.4. Consider the generalization of the delay fractional-order version of the Chen system [3] which involves the variable order:

$D_t^{\alpha \left(t\right)}x\left(t\right)=a(y\left(t\right)-x\left(t-\tau )\right),$

$D_t^{\alpha \left(t\right)}y\left(t\right)=\left(c-a\right)x\left(t-\tau \right)-x\left(t\right)z\left(t\right)+cy(t),$

$D_t^{\alpha \left(t\right)}z\left(t\right)=x\left(t\right)y\left(t\right)-bz\left(t-\tau \right),$

$x\left(t\right)=0.2,y\left(t\right)=0,z\left(t\right)=0.5fort\in \left[-\tau ,0\right],$E32

on the interval [0,30] and with step size h = 0.001.

At = 1, Figure 27 shows the numerical solution for τ = 0.005 and shows chaotic yz phase portrait for this case. Figure 28 shows the numerical solution and chaotic yz phase portrait at α = 0.97 and τ = 0.005. Moreover, for the variable order, we choose different cases for α(t):

At α(t) = 0.97 – (0.01/100)t, Figure 29 shows the numerical solution y(t) and chaotic yz phase portrait of this example for τ = 0.005. When we increase the value of τ to 0.015, the results of the numerical solution and chaotic yz phase portrait become stable as shown in Figure 30. At α(t) = 0.94 – (0.01/100)t, Figure 31 shows the numerical solution for τ = 0009, and a limit cycle phase portrait was observed. When we increase the value of τ to 0.011, the chaotic yz phase portrait becomes stable as shown in Figure 32. Figure 33 shows the numerical solution for τ = 0.005 and chaotic yz phase portrait at α(t) = 89 − (0.01/100)t. It is observed that even in one-dimensional delayed systems of variable order, chaotic behavior can be shown, and subjected to some critical order, the system changes its nature and becomes periodic. In some cases, it is observed that the phase portrait gets stretched as the order of the derivative is reduced. According to the numerical test examples, it can be concluded that the proposed numerical technique is a powerful technique to calculate approximate solution of VOFDDEs.

Figure 27.

The numerical behavior and chaotic attractors for = 1, and τ = 0.005.

Figure 28.

The numerical behavior and chaotic attractors for α = 0.97 and τ = 0.005.

Figure 29.

The numerical behavior and chaotic attractors for α(t) = 0.97 – (0.01/100)t and τ = 0.005.

Figure 30.

The numerical behavior and chaotic attractors for α(t) = 0.97 – (0.01/100)t and τ = 0.015.

Figure 31.

The numerical behavior and chaotic attractors for α(t) = 0.94 – (0.01/100)t and τ = 0.009.

Figure 32.

The numerical behavior and chaotic attractors for α(t) = 0.94 – (0.01/100)t and τ = 0.011.

Figure 33.

The numerical behavior and chaotic attractors for α(t) = 0.89 – (0.01/100)t and τ = 0.005.