1. Introduction
There has been a great deal of interest in working with nonlinear ordinary differential equations (ODEs) and dynamical systems over the past decades. Diverse physical as well as nonphysical phenomena including economical, epidemiological, and biological models belong to these category. Among others, nonlinear systems with chaotic behavior are of particularly of interest. Chaos theory of nonlinear dynamics is a branch of mathematical science that deals with highly complex systems, where small changes now in their inputs cause large changes in their outputs later [
1]. The notion of chaos was first introduced to the modern world by Edward Lorenz in 1972 with conceptualization of “Butterfly Effect”. By understanding this theory we are able to predict the behavioral representation of a complex system more accurately. Chaotic systems are unstable since they tend not to resist any outside disturbances but instead react in significant manner. As a meteorologist, Lorenz developed a mathematical model for the air movement in the atmosphere. It caused various differences in the outcome of the given model. In this manner he discovered the principle of sensitive dependence on initial conditions, which is now considered as a key element in any chaotic system [
2]. Indeed, systems with chaotic behavior are common in nature and diverse real-life phenomena can also be classified in the framework of a chaotic system. Beside of their appearance in meteorology, they can be found also in solar system, chemical reactions, climate dynamic, population dynamic, economics, and so on.
In the economic setting, Ma and Chen [
3,
4] were one of the pioneering who studied and analyzed the chaotic economic systems. The study of memory effects and fractional-order derivative on a chaotic financial system was given by Chen [
5]. The concept of fractional calculus provides more insight in predicting the behavior of complex dynamical systems especially those with chaotic characteristics, see [
6]. Various fractional order systems with chaos have been considered in literature such as Chua [
7], Chen [
8], Rössler [
9], Lorentz [
10] to name a few. Finding the exact or analytical solutions of the most fractional-order systems cannot be possible. Alternately, the researchers have been tried to develop effective semi-analytical as well as numerical procedures to deal with such systems of fractional order. Among diverse existing proposed methods, we refer to the recent published works including the spectral collocation approaches [
11,
12,
13], the fractional Adams-Bashforth approach [
14], the residual power series method and Laplace technique [
15], the variational iteration scheme [
16], and the wavelet technique [
17] to name a few.
The main objective of the current research is to investigate a real-life model arising in the financial market. The model equation consists of three differential equations with quadratic nonlinearities given by [
5]
where
and
for
are constants. Here,
denotes the saving amount rate,
represents the costs per investment, and
shows the elasticity of demand of commercial markets. In the model (
1) by
we denote the fractional derivative operator of order
in the sense of Liouville–Caputo (LC). Let
, and
be given real numbers. The supplemented initial conditions are as follow
It should be emphasized that the dynamic analysis and chaotic intrinsic nature of the financial system (
1) have already been investigated by researchers, see cf. [
5].
Based on our knowledge, only a few works have been proposed to solve the system (
1) numerically. The Adams-Bashforth-Moulton scheme used in [
18] to deal with the variable-order fractional counterpart of (
1). The integer-order of the above model solved by a robust spectral Chebyshev method [
19]. The fractional-order of model (
1) in the sense of both Caputo and Atangana-Baleanu-Caputo (ABC) derivatives considered in [
20] and solved by the modified homotopy analysis transform method. They also proved the existence and uniqueness of the solutions based on Picard-Lindelof’s approach. The variable-order of the underlying model in the sense of ABC derivative invsetiagted in [
21]. A computational stochastic approach based on the artificial neural networks introduced recenlty in [
22]. Furthermore, the authors in [
23] proposed a new robust nonlinear controller that stabilizes the chaotic finance system (
1). Some other financial models related to (
1) with some proposed approximate methods were given in [
24,
25,
26,
27].
In this work, we solve the nonlinear system (
1) with two spectral collocation approaches based on (novel) modified Bessel functions. To preserve the high accuracy of the proposed methods on a long computational domain, we further employ the domain-splitting methodology. In the first and direct collocation procedure, we convert the financial model into an algebraic system of nonlinear equations. On the other hand, in the next and more efficient algorithm, we utilize first the methodology of quasilinearization to (
1) and then employ the former direct technique to the resulting sequence of the linearized system of equations. In the latter technique we arrive at a family of linear matrix equations to be solved easily. Using the domain-splitting methodology enables us to keep the accuracy high on a long time domain. For this purpose, it is sufficient to divide the domain into several subdomains and then employ the traditional collocation schemes on each subintervals more accurately.
The balance of the current work is given next. A review of fundamental facts on fractional integral and derivatives is provided in
Section 2.
Section 3 contains of introducing the fundamental definitions and main properties of the modified Bessel functions. The generalized form of these polynomials is introduced next. The error estimate of the given basis functions is also established. In
Section 4, we first develop the direct GMBPs collocation based on domain-splitting strategy. Afterward, the hybrid QLM-GMBPs procedure relied on the quasilinearization process is implemented in detail. We also testify the accuracy of two proposed collocation technique via defining the technique of residual error functions in this section. Experiments and computational results are carried out in
Section 5. A concise conclusion of the current manuscript is discussed in
Section 6.
2. Basic Preliminaries on Fractional Calculus
Here and in this section, our aim is to provide some facts associated to fractional calculus. The definition of Liouville–Caputo (LC) is first given and some of its properties will be mentioned. At the end, we provide a theorem for the generalized Taylor’s formula with LC operator.
Definition 1 ([
6]).
Let suppose that is a continuously differentiable function r-times. The LC fractional operator of order of is defined as followswhere Let
and
be two constants. The LC operator has linear property in the sense that
where two functions
for
are
r-times continuously differentiable. In addition to the above linear property, we utilize the next properties of LC operator will be used frequently throughout the paper as follow
where
,
and
denote the floor and ceiling functions of
respectively. The next Theorem states the generalized form of Taylor series expansion when we have the LC fractional derivative operator. To see a proof, we refer the readers to [
28]. Towards this end, let us have the convention that
Theorem 1. Let and the functions are continuous on for . Then, the series representation of is written as followsfor some . The immediate consequence of the former Theorem is given next.
Corollary 1 (Fractional Taylor inequality).
Let assume that the hypotheses of Theorem 1 hold and we have . Then, we get the upper bound for remainder of Taylor’s formula as 4. Domain-Splitting QLM-GMBPs/GMBPs Collocation Techniques
As previously mentioned in the
Section 1, the spectral collocation procedure may be failed to converge if the computational domain of interest is large. To overcome this inaccuracy, the idea is to split the given (large) domain into a sequence of sub-domains and then apply the spectral approach in each subdomain. To this end, we consider a (uniform) partitioning of the interval
into
subintervals
where
and
. The length of each
will denote by
and is equal to
. When
, we recover the traditional spectral collocation strategy, see cf. [
31,
32,
33,
34,
35,
36]. In each subdomain
we now consider the approximate solution to be in the form (
15) as
where
is the vector of unknown coefficients and
is defined in (
15). After obtaining the approximate solutions for
, the desired approximation on the whole domain
can be expressed as
where the characteristic function
is given by
Furthermore, at each subinterval
, we utilize the followinng collocation nodes
Finally, we emphasize that the initial conditions on the first subinterval
are taken as the original initial conditions for our model problem (
1) whilst on the remaining subintervals
,
we will use the approximate solutions obtained on the previous subintervals
. The fundamental idea of the spectral collocation for the model (
1) will be illustrated next. Towards this end, we only consider the proposed approach on an arbitrary domain
for
.
4.1. The GMBPs Collocation Approach
Let us first represent the financial system (
1) into a matrix format on the subinterval
. It can be expressed as
Here, we have utilized the following notations
Furthermore, the nonlinear terms
for
are
Now, the task is that to state the unknown solutions
and
as a finite combination of cutted series form (
17) using (
) GMBPs bases. It follows that
To find the coefficients
for
, we apply the spectral matrix collocation procedure by using the GMBPs. A decomposition of the vector
is provided next:
Lemma 1. We can write the vector as the followingwhererepresents the vector of standard bases and the second matrix is a non-singular matrix whose components are given in (12). It has a triangular structure and given by Proof. By expanding the relation (
12) in terms of powers of
t followed by using the definitions of
and
we conclude the result easily. □
By inserting the expression
in (
21) into (
20) we render
The next job is to compute the
-derivative of the aforementioned approximate solutions. It is evident that we need to calculate only the
-derivative of
. Towards this end, we exploit the property (
4) as well property (5) to calculate
, which is denoted by
. Note that we can calculate it algorithmically with linear complexity
, see Algorithm 3.1 in [
38]. Therefore, we get
In the vector format, we can approximate the true solutions of system (
19) in the form
By introducing the following matrix notations
we are able to rewrite
and
in the forms of matrices. It follows that
where we have used the notation
.
We now insert the collocation nodes (
18) of the subinterval
into the matrix represenation form (
19). The next set of equations are obtained
If we introduce the next notations
we are able to represent the system of equations in a concise form as
Our goal is to write the involved vectors in (
27) in a matrix format.
Lemma 2. In matrix representation form, two vectors and in (25) at the collocation nodes (18) can be written asHere, two matrices and are defined aswhere three matrices as well as the vector are introduced in (25). The proof of the foregoing Lemma can be done easily by just considering relations (
25) followed by inserting the nodes (
18) into them. Similarly, two nonlinear terms
for
still need to be written in the matrix expressions. We have
Lemma 3. We can express in (27) for in the formsHere, we have the representations , and so thatAlso, we have , where . Proof. To prove the first assertion, let us note the following relation
In accordance to relations (
22), it is no difficult to show that the first diagonal matrix is written as
and the second column matrix can be expressed as
. Similar arguments hold for
, where the second column vector is terms of
instead of
. □
We now turn our attention to (
27) and insert all obtained matrix forms into it. To be more precise, we replace
, and
by the corresponding representation form in relations (
28) and (
29). We reach at the fundamental matrix equation, which has the representation
where
It should be noted that the resultant matrix Equation (
30) is a nonlinear system with
unknowns
for
and
to be determined by our proposed algorithm.
Next, we incorporate the intial conditions (
2) into (
30). Towards this end, we provide the matrix form of (
2). By tending
in the first relation in (
25) we reach at the next form
Let us emphasize that the three constants
and
are at hand by (
2). The rows one to third of the augmented matrix
will now be substituted by the above three rows
. Thus, we reach at the modified and final version of fundamental matrix equation as
We solve (
31) to acquire the approximate solution of the financial model (
1) on the subintervals
for
. As stated before, for the first subinterval
, we use the original initial conditions (
2). By iterating on
m, the approximate solutions on
are obtained by using the previously approximate solutions on
as the initial conditions evaluated at the starting point of
.
4.2. Computing Errors via Residual Error Function Technique
The true exact form of the solutions of the financial system (
1) are not yet known. To check the quality of approximations, we compute the residual error functions (REFs) rather than the difference between the exact and obtained numerical solutions. Towards this end, we first calculate the three approximations
,
, and
via the proposed GMBPs matrix collocation algorithm as described above. Afterwards, we insert these approximations into the financial system (
1) followed by defining the REFs as the difference between the left and right hand sides of the involved equations. In the other words, we set the REFs for (
1) on the subintervals
as
4.3. The QLM-GMBPs Collocation Approach
The direct matrix collocation procedure based on the GMBPs is illustrated in the last part for the nonlinear chaotic system (
1). We expect that a combination of the proposed approach with the domain-splitting leads to a better accuracy of the obtained approximate solutions. However, using a large value for the number of basis function (
P) may cause of some inefficiencies during the execution of the algorithm, which comes from the intrinsic nonlinearity of the model under consideration. The idea of linearization has usually been worked to accelerate the running time of proposed algorithms. Applying the technique of quasi-lineariztion method (QLM) has been successfully proved by many published works, see cf. [
39,
40,
41]. The QLM aimed to not only linearize the model but also keep the accuracy of the proposed algorithm in the same level as the directly employed algorithm.
Let us we rewrite the original model (
1) in the form
where
By starting the initial rough guess
for the solution
of system (
33), we propose the following iterative QLM as follows
where the matrix
represents the Jacobian matrix of size
. Applying the linearization, we get
Here, we have utilized the notations
The above system has also the initial conditions
in accordance to (
2).
Now, the direct GMBPs matrix approach can be applied to the set of linearized equations (
34) without any difficulties stem from the nonlinear terms in (
19). Supposedly, we have already the approximate solution
available for the solution
of (
34) in the iteration number
. Then, we take the next approximation on
as
where
for
are the vectors of unknowns in the iteration
. Moreover, the vector
and matrix
are defined in (
21). We continue by inserting the collocation nods (
18) into (
34) followed by introducing the following matrix and vector notations
By employing the relations (
28), we reach at the next fundamental matrix equation
where
for
on
. Here, the matrices
,
, and
are given in (
28). In analogue to the direct approach, we enter the initial
into the matrix equation (
36). By solving the alternate matrix equation we arrive at the unknown coefficients
efficiently and our job is finished.
Finally, let us mention that the REFs related to the QLM-GMBPs approach are obtained as (
32). We denote them by
for
defined on
and on the iteration
. For a fixed value of
e, we further calculate the numerical order of convergence related to the achieved REFs in the infinity norm as
Hence, we set
5. Experimental and Graphical Results
The financial system (
1) will now be solved by using the GMBPs/QLM-GMBPs collocation techniques as described in the last parts. In order to validate our results numerically, different model parameters will be considered to show the performance of the proposed spectral collocation approaches. Moreover, various values of fractional order
are used in the simulations to show the applicability of the presented algorithms. All the experimental simulations are performed by utilizing Matlab 2021a executed on a personnel laptop. Also, note that employing
gives us sufficiently accurate results in the QLM-GMBPs collocation procedure. The starting vector function
is taken as the zero function.
Example 1. To begin, we consider the parameters related to the chaotic financial system (1) taking from [22] as followThe initial conditions are given by We set
and take
in the computations. We first focus on the performance of the novel QLM-GMBPs collocation technique based on domain-splitting strategy compared to traditional version of this algorithm. For this purpose, we consider
and
and
in these two algorithms respectively. The result of approximations are depicted in
Figure 1 and
Figure 2 for the three solutions
, and
. For validation, we also plot the results obtained via ode45 routine in Matlab. It is evident that the outcomes of domain-splitting are accurate than those obtained by using
especially for the interest rate
. The gap between the results of both approaches will increase if one takes a longer domain of computation.
We next show that the results of direct and indirect collocation schemes i.e., the GMBPs and QLM-GMBPs matrix algorithms are coincided. Let us take
and set the parameters as above. We only report the results of the direct GMBPs collocation matrix approach for
as follow
It should be noted that the outcomes of QLM-GMBPs technique are very close to the above approximations and we omit them to save space. However, the achieved REFs related to both approaches are visualized in
Figure 3 graphically. On the left plot, we show the REFs of direct GMBPs collocation method using
while on the right picture we depict the REFs utilizing
in the QLM-GMBPs matrix technique. It is clearly seen that in the direct approach we have problem with the convergence of the method if we increase
P. The reason may be the fact that the used nonlinear solver fails to give us a reasonable approximation. On the other hand, we can utilize larger values of
P in the QLM-GMBPs technique to reach at more accurate results. Moreover, the latter iterative method has more efficacy than the direct method. For instance, for
, the required time to solve the fundamental matrix Equation (
31) is about
s while the elapsed time for the corresponding matrix equation via QLM-GMBPs procedure is approximately
s. Consequently, in the following experiments, we usually run the latter numerical approach due to its efficacy.
We further show the utility of QLM-GMBPs approach by using
and a relatively large interval
. The results of approximations on the first and last subibtervals
and
are obtained as follow
for
and
for
. In
Figure 4 and
Figure 5 we depict all approximate solutions on the domains
for
. The results of the ode45 are also plotted to show the validity of our approximations. By looking at these figures we find that the results of both numerical strategies are overlapped. We get a higher order accuracy if we increase the number of bases or use a lager value of subdivisions
M. We next take
. The 3D pictorial representation of uvw-plot as well as the 2D depiction of uv-plot are displayed in
Figure 6, which show the chaotic dynamical behavior of the financial system (
1). Furthermore, the similar results related to 2D representations of uw/vw-plot are shown in the next
Figure 7. For validation, the numerical approximations using the ode45 are also displayed in this figure for comparison purposes. It is clear that the outcomes of both approaches are coincided.
We now pay attention to the non-integer cases and consider
. Using the domain-splitting approach along with the QLM-GMBPs strategy we plot the approximate solutions in
Figure 8. We utilize the values of parameter
to be both 1 and equal to
. As shown in this figure, the results obtained to
and
are very close together. However, the REFs achieved by using
are smaller in magnitude compared to those obtained by
as presented in
Figure 9. Note that one gets more accurate results by increasing
P. Towards this end, we employ different values of
P in the next experiments and report the maximum values of REFs and their associated order of convergence, i.e.,
for
. These simulation results are tabulated in
Table 1. We use
,
, and
. The parameter
is taken as
. Clearly, the exponential order of convergence is visible for the results presented in
Table 1. Of course, a slightly higher order of accuracy is seen for the case
compared to
.
Finally, we employ diverse values of fractional orders
and visualize the numerical simulations for the three approximated solutions of system (
1). We utilize
,
, and the parameter
is equal to one. These experimental results are depicted in
Figure 10 and
Figure 11. We also plotted the results related to
in these figures.