The Bernstein Operational Matrices for Solving the Fractional Quadratic Riccati Differential Equations with the Riemann-Liouville Derivative

Abstract

We obtain the approximate analytical solution for the fractional quadratic Riccati differential equation with the Riemann-Liouville derivative by using the Bernstein polynomials (BPs) operational matrices. In this method, we use the operational matrix for fractional integration in the Riemann-Liouville sense. Then by using this matrix and operational matrix of product, we reduce the problem to a system of algebraic equations that can be solved easily. The efficiency and accuracy of the proposed method are illustrated by several examples.

1. Introduction

The Riccati differential equation is named after the Italian Nobleman Count Jacopo Francesco Riccati (1676–1754). The book of Reid [1] contains the fundamental theories of the Riccati equation, with applications to random processes, optimal control, and diffusion problems. Moreover, it is well known that the one-dimensional static Schrödinger equation is closely related to a Riccati differential equation [2]. Solitary wave solution of a nonlinear partial differential equation can be represented as a polynomial in two elementary functions satisfying a projective Riccati equation [3].

In this paper, we are dealing with the fractional quadratic Riccati differential equation as follows: _𝑎𝐷^𝛼_𝑡𝑦(𝑡)=𝑓(𝑡)+𝑏(𝑡)𝑦(𝑡)+𝑎(𝑡)𝑦²(𝑡),𝑎<𝑡≤𝑏, (1) subject to the initial conditions 𝑦^(𝑘)(𝑎)=𝑦^𝑘_𝑎,𝑘=0,1,…,⌈𝛼⌉−1, (2) where 𝑦^𝑘_𝑎??(𝑘 = 0,1,…,⌈𝛼⌉−1) are constants and 𝑎(𝑡), 𝑏(𝑡), and 𝑓(𝑡) are known functions. We can see [4] to guarantee the existence and uniqueness of the solution of (1) with initial conditions (2).

The general response expression (1) contains a parameter 𝛼, the order of the fractional derivative that can be varied to obtain various responses. In the case that 𝛼 is integer, then (1) is reduced to the classical Riccati differential equation.

This problem to develop the analytical and numerical method to solve the Riccati differential equation with standard derivative, the Caputo fractional derivative and the Riemann-Liouville fractional derivative, has attracted much attention and has been studied by many authors [5–18].

The aim of this work is using the Bernstein polynomials for solving the problem (1) and (2). We notice that the problem presented [18] was in the Caputo sense but in our work, the problem is with the Riemann-Liouville derivative; therefore we considered a more general space of functions. Also, in [18], the authors used the polynomials in the form of 𝐵^𝛼_𝑖,𝑚(𝑥) = (^𝑚_𝑖)𝑥^𝛼𝑖(1−𝑥^𝛼)^𝑚−𝑖??(𝑖 = 0,1,…,𝑚) that is different from the standard Bernstein polynomials. So, the operational matrices in this work are different from those in [18].

The organization of this paper is as follows. In Section 2, the Bernstein polynomials are introduced. Some basic definitions and properties of the fractional calculus and also the BPs operational matrix for the Riemann-Liouville fractional integration are presented in Section 3. In Section 4, by BPs operational matrices, we solve the fractional quadratic Riccati differential equation. In Section 5, we discuss the convergence of the proposed method. In Section 6, several examples are considered to evaluate the power and effectiveness of the presented method. Some conclusions are summarized in the last section.

2. The Bernstein Polynomials and Their Properties

On the interval [0,1] we define the Bernstein polynomials (BPs) of mth degree as follows [19]: 𝐵_𝑖,𝑚(𝑥)=𝑚𝑖𝑥^𝑖(1−𝑥)^𝑚−𝑖,𝑖=0,1,…,𝑚. (3) Set {𝐵_0,𝑚(𝑥),𝐵_1,𝑚(𝑥),…,𝐵_𝑚,𝑚(𝑥)} in the Hilbert space 𝐿²[0,1] is a complete basis. We can write Φ_𝑚(𝑥)=𝐴𝑇_𝑚(𝑥), where 𝐴 is a matrix upper triangular, 𝑇_𝑚(𝑥)=[1,𝑥,…,𝑥^𝑚]^𝑇, and Φ_𝑚(𝑥)=[𝐵_0,𝑚(𝑥),𝐵_1,𝑚(𝑥),…,𝐵_𝑚,𝑚(𝑥)]^𝑇 [20].

As a result, any polynomial of degree 𝑚 can be expanded in terms of linear combination of 𝐵_𝑖,𝑚(𝑥)??(𝑖=0,1,…,𝑚) as given below: 𝑃(𝑥)=^𝑚_𝑖=0𝑐_𝑖𝐵_𝑖,𝑚(𝑥)=𝑐^𝑇Φ_𝑚(𝑥), (4) where 𝑐=¹₀Φ(𝑥)Φ(𝑥)^𝑇𝑑𝑥⁻¹¹₀𝑃(𝑥)Φ(𝑥)𝑑𝑥. (5) The approximation of functions within the Bernstein polynomials and convergence analysis can be found in [20, 21].

3. BPs Operational Matrix for the Riemann-Liouville Fractional Integration

In this section, firstly, we give some basic definitions and properties of the fractional calculus which are used further in this paper.

Definition 1 (see [4, 22–24]). Let 𝛼≥0; the operator _𝑎𝐼^𝛼_𝑡, defined on 𝐿₁[𝑎,𝑏] by _𝑎𝐼^𝛼_𝑡𝑓(𝑡)=1Γ(𝛼)^𝑡_𝑎(𝑡−𝑥)^𝛼−1𝑓(𝑥)𝑑𝑥,𝑎≤𝑡≤𝑏,_𝑎𝐼⁰_𝑡𝑓(𝑡)=𝑓(𝑡), (6) is called the Riemann-Liouville fractional integral operator of order 𝛼.

Definition 2 (see [4, 22–24]). Let 𝛼≥0 and 𝑛=⌈𝛼⌉; the operator _𝑎𝐷^𝛼_𝑡, defined by _𝑎𝐷^𝛼_𝑡𝑓(𝑡)=1Γ()𝑑𝑛−𝛼𝑑𝑡^𝑛𝑡_𝑎(𝑡−𝑥)^{𝑛−𝛼−1}𝑓(𝑥)𝑑𝑥,𝑎≤𝑡≤𝑏,_𝑎𝐷⁰_𝑡𝑓(𝑡)=𝑓(𝑡), (7) is called the Riemann-Liouville fractional derivative operator of order 𝛼.

Definition 3 (see [4, 22–24]). Let 𝛼≥0, 𝑛=⌈𝛼⌉, and 𝑑^𝑛𝑓(𝑥)/𝑑𝑥^𝑛∈𝐿₁[𝑎,𝑏]. The operator ^𝑐_𝑎𝐷^𝛼_𝑡, defined by ^𝑐_𝑎𝐷^𝛼_𝑡𝑓(𝑡)=1Γ(𝑛−𝛼)^𝑡_𝑎(𝑡−𝑥)^{𝑛−𝛼−1}𝑑^𝑛𝑓(𝑥)𝑑𝑥^𝑛𝑑𝑥,𝑎≤𝑡≤𝑏,_𝑐𝑎𝐷⁰_𝑡𝑓(𝑡)=𝑓(𝑡), (8) is called the Caputo fractional derivative operator of order 𝛼.

Lemma 4. If 𝛼≥0, 𝑛=⌈𝛼⌉, and 𝑎≤𝑡≤𝑏, then (1)^𝑐_𝑎𝐷^𝛼_𝑡𝑎𝐼^𝛼_𝑡𝑓(𝑡)=𝑓(𝑡),(2)_𝑎𝐼^𝛼_𝑡^𝑐_𝑎𝐷^𝑎_𝑡𝑓(𝑡)=𝑓(𝑡)−^𝑛−1_𝑘=0𝑓^(𝑘)(𝑎)()𝑘!𝑥−𝑎^𝑘,(3)^𝑐_𝑎𝐷^𝛼_𝑡𝑓(𝑡)=_𝑎𝐷^𝛼_𝑡𝑓(𝑡)−^𝑛−1_𝑘=0𝑓^(𝑘)(𝑎)Γ()()𝑘−𝛼+1𝑥−𝑎^𝑘−𝛼. (9)(10)(11)

Proof. See [22–24].

Theorem 5. One can get BPs operational matrix 𝐹_𝛼 from order (𝑚+1)×(𝑚+1) for the Riemann-Liouville fractional integral as _𝑎𝐼^𝛼_𝑡Φ_𝑚(𝑡)=1Γ(𝛼)^𝑡_𝑎(𝑡−𝑥)^𝛼−1Φ_𝑚(𝑥)𝑑𝑥≈𝐹_𝛼Φ_𝑚(𝑡). (12)

Proof. See [21].

4. BPs for Solving the Fractional Quadratic Riccati Differential Equation

Firstly, we use the initial conditions to reduce a given initial-value problem to a problem with zero initial conditions.

So, we define 𝑦(𝑡)=𝑦(𝑡)+𝑧(𝑡), (13) where 𝑦(𝑡) is some known function that satisfied the initial conditions (2) and 𝑧(𝑡) is a new unknown function.

Substituting (13) in (1) and (2), we have an initial-value problem as follows: _𝑎𝐷^𝛼_𝑡𝑧(𝑡)=𝑓(𝑡)+𝑏̃(𝑡)𝑧(𝑡)+𝑎(𝑡)𝑧²(𝑡),𝑎<𝑡≤𝑏, (14) subject to the initial conditions 𝑧^(𝑘)(𝑎)=0,𝑘=0,1,…,⌈𝛼⌉−1. (15) On the other hand, by (11) in Lemma 4 we can write ^𝑐_𝑎𝐷^𝛼_𝑡𝑧(𝑡)=_𝑎𝐷^𝛼_𝑡𝑧(𝑡). (16) Also, by using Lemma 3.3 in [20] the inputs 𝑓(𝑡), 𝑎(𝑡), 𝑏̃(𝑡), and _𝑎𝐷^𝛼_𝑡𝑧(𝑡) can be approximated as follows: 𝑓(𝑡)≈𝐹^𝑇Φ_𝑚(𝑡),_𝑎𝐷^𝛼_𝑡𝑧(𝑡)≈𝐶^𝑇Φ_𝑚(𝑡),𝑎(𝑡)≈𝐴^𝑇Φ_𝑚(𝑡),𝑏̃(𝑡)≈𝐵^𝑇Φ_𝑚(𝑡), (17)(18)(19)(20) where 𝐹, 𝐴, and 𝐵 are known (𝑚+1) column vectors and 𝐶 is an unknown (𝑚+1) column vector.

From (10), (15), (16), (18), and (12), we have 𝑧(𝑡)₍₎=_𝑎10𝐼^𝛼_𝑡^𝑐_𝑎𝐷^𝛼_𝑡𝑧(𝑡)+^𝑛−1_𝑘=0𝑧^(𝑘)(𝑎)()𝑘!𝑥−𝑎^𝑘=_()𝑎15𝐼^𝛼_𝑡^𝑐_𝑎𝐷^𝛼_𝑡𝑧(𝑡)₍₎=_𝑎16𝐼^𝛼_𝑡𝑎𝐷^𝛼_𝑡𝑧(𝑡)≈_()𝑎18𝐼^𝛼_𝑡𝐶^𝑇Φ_𝑚(𝑡)=𝐶^𝑇_𝑎𝐼^𝛼_𝑡Φ_𝑚(𝑡)≈₍₎₁₂𝐶^𝑇𝐹_𝛼Φ_𝑚(𝑡)=𝐶^𝑇_𝛼Φ_𝑚(𝑡), (21) where 𝐶^𝑇_𝛼=𝐶^𝑇𝐹_𝛼.

Now, by substituting (17)–(21) into (14), we obtain 𝐶^𝑇Φ_𝑚(𝑡)=𝐹^𝑇Φ_𝑚(𝑡)+𝐵^𝑇Φ_𝑚(𝑡)Φ_𝑚(𝑡)^𝑇𝐹^𝑇_𝛼𝐶+𝐴^𝑇Φ_𝑚(𝑡)×𝐶^𝑇𝐹_𝛼Φ_𝑚(𝑡)Φ_𝑚(𝑡)^𝑇𝐹^𝑇_𝛼𝐶=𝐹^𝑇Φ_𝑚(𝑡)+𝐵^𝑇Φ_𝑚(𝑡)Φ_𝑚(𝑡)^𝑇𝐶_𝛼+𝐴^𝑇Φ_𝑚(𝑡)×𝐶^𝑇_𝛼Φ_𝑚(𝑡)Φ_𝑚(𝑡)^𝑇𝐶_𝛼. (22) Then, from Lemma??3.5 in [20] we have 𝐶^𝑇_𝛼Φ_𝑚(𝑥)Φ_𝑚(𝑥)^𝑇≈Φ_𝑚(𝑥)^𝑇𝐶_𝛼,𝐴^𝑇Φ_𝑚(𝑥)Φ_𝑚(𝑥)^𝑇≈Φ_𝑚(𝑥)^𝑇𝐴,𝐵^𝑇Φ_𝑚(𝑥)Φ_𝑚(𝑥)^𝑇≈Φ_𝑚(𝑥)^𝑇𝐵. (23)(24)(25) Therefore we can reduce (22) by (23)–(25) as 𝐶^𝑇Φ_𝑚(𝑡)=𝐹^𝑇Φ_𝑚(𝑡)+Φ_𝑚(𝑡)^𝑇𝐵𝐶_𝛼+Φ_𝑚(𝑡)^𝑇𝐴𝐶_𝛼𝐶_𝛼. (26) Finally, we obtain the following nonlinear system of algebraic equation: 𝐶=𝐹+𝐵𝐶_𝛼+𝐴𝐶_𝛼𝐶_𝛼, (27) such that by solving this system we can obtain the vector 𝐶. Then, we can get 𝑦(𝑡)≈𝑦(𝑡)+𝐶^𝑇𝐹_𝛼Φ_𝑚(𝑡). (28)

5. Convergence Analysis

In this section, we investigate the convergence analysis for the method presented in Section 4.

The problem (14) changes to the following problems _𝑎𝐷^𝛼_𝑡𝑧(𝑡)=𝑓(𝑡)+𝑏̃(𝑡)_𝑎𝐼^𝛼_𝑡𝑎𝐷^𝛼_𝑡𝑧(𝑡)+𝑎(𝑡)_𝑎𝐼^𝛼_𝑡𝑎𝐷^𝛼_𝑡𝑧(𝑡)²,𝑎<𝑡≤𝑏, (29) since 𝑧(𝑡)=_()𝑎10𝐼^𝛼_𝑡^𝑐_𝑎𝐷^𝛼_𝑡𝑧(𝑡)+^𝑛−1_𝑘=0𝑧^(𝑘)(𝑎)()𝑘!𝑥−𝑎^𝑘₍₎=_𝑎15𝐼^𝛼_𝑡^𝑐_𝑎𝐷^𝛼_𝑡𝑧(𝑡)=_()𝑎16𝐼^𝛼_𝑡𝑎𝐷^𝛼_𝑡𝑧(𝑡). (30) By taking 𝑢(𝑡)=_𝑎𝐷^𝛼_𝑡𝑧(𝑡) we obtain the following fractional integral equation: 𝑢(𝑡)=𝑓(𝑡)+𝑏̃(𝑡)_𝑎𝐼^𝛼_𝑡𝑢(𝑡)+𝑎(𝑡)_𝑎𝐼^𝛼_𝑡𝑢(𝑡)²,𝑎<𝑡≤𝑏. (31) If we use the approximation 𝑢(𝑡)≈𝐶^𝑇Φ_𝑚(𝑡), then the problem (31) from space 𝐶¹[0,1] reduces to the following problem in space 𝑆_𝑚=Span{𝐵_0,𝑚(𝑡),𝐵_1,𝑚(𝑡),…,𝐵_𝑚,𝑚(𝑡)}: 𝐶^𝑇Φ_𝑚(𝑡)=𝑓(𝑡)+𝑏̃(𝑡)_𝑎𝐼^𝛼_𝑡𝐶^𝑇Φ_𝑚(𝑡+𝑎(𝑡)_𝑎)𝐼^𝛼_𝑡𝐶^𝑇Φ_𝑚(𝑡)²,𝑎<𝑡≤𝑏. (32) Now, similar to Theorem??6.1 in [21], we propose the next theorem.

Theorem 6. Suppose that 𝑢^∗(𝑡)∈𝐶¹[0,1] is the exact solution of (31) and 𝜇_𝑚=𝐽[𝑢_𝑚]=Min_{𝑢∈𝑆𝑚}𝐽[𝑢], where 𝐽[𝑢]=𝑢(𝑡)−𝑓(𝑡)−𝑏̃(𝑡)_𝑎𝐼^𝛼_𝑡𝑢(𝑡)−𝑎(𝑡)_𝑎𝐼^𝛼_𝑡𝑢(𝑡)². (33) Then one has 𝜇_𝑚→0 as𝑚→∞ (i.e.,??𝑢_𝑚(𝑡)→𝑢^∗(𝑡) as 𝑚→∞).

Proof. By substituting (33) instead of (41) in [21], we can use the proof of Theorem??6.1 in [21].

Theorem 7. Suppose that 𝑢^∗(𝑡)∈𝐶¹[0,1] is the exact solution of (31) and 𝑢_𝑚(𝑡)∈𝑆_𝑚 is the obtained solution of (26). Then one has 𝑢_𝑚(𝑡)→𝑢^∗(𝑡) as𝑚→∞.

Proof. Substituting (17), (19), (20), and (12) in (32) we have 𝐶^𝑇Φ_𝑚(𝑡)=𝐹^𝑇Φ_𝑀(𝑡)+𝑒^𝑀𝑓+𝐵^𝑇Φ_𝑀(𝑡)+𝑒^𝑀̃𝑏𝐶^𝑇𝐹_𝛼Φ_𝑀(𝑡)+𝐸^𝑀_𝐼+𝐴^𝑇Φ_𝑀(𝑡)+𝑒𝑎^𝑀𝐶^𝑇𝐹_𝛼Φ_𝑀(𝑡)+𝐸^𝑀_𝐼². (34) From Lemmas 2.3 and??4.1 in [21], we have 𝑒^𝑀𝑓,𝑒^𝑀̃𝑏,𝑒𝑎^𝑀,𝐸^𝑀_𝐼→0 as 𝑀→∞. So we can observe that as 𝑀 increases, (34) gets close to (32). Now, by taking 𝑀=𝑚 we propose the following problem that gets close to (32) as 𝑚 increases: 𝐶^𝑇Φ_𝑚(𝑡)=𝐹^𝑇Φ_𝑚(𝑡)+𝐵^𝑇Φ_𝑚(𝑡)Φ_𝑚(𝑡)^𝑇𝐶_𝛼+𝐴^𝑇Φ_𝑚(𝑡)𝐶^𝑇_𝛼Φ_𝑚(𝑡)Φ_𝑚(𝑡)^𝑇𝐶_𝛼. (35) Then by (23) and (25), (35) reduces to the following equation: 𝐶^𝑇Φ_𝑚(𝑡)=𝐹^𝑇Φ_𝑚(𝑡)+Φ_𝑀(𝑡)^𝑇𝐵+𝐸^𝑀_𝐵𝐶_𝛼+𝐴^𝑇Φ_𝑚(𝑡)Φ_𝑀(𝑡)^𝑇𝐶_𝛼+𝐸^𝑀_𝐶𝛼𝐶_𝛼. (36) Equation (36) gets to (35) as 𝑀→∞, because from Lemma??3.1 in [21] 𝐸^𝑀_𝐵,𝐸^𝑀_𝐶𝛼→0 as 𝑀→∞. Then by deleting 𝐸^𝑀_𝐵, 𝐸^𝑀_𝐶𝛼, taking 𝑀=𝑚, and using (24) in (36), we have 𝐶^𝑇Φ_𝑚(𝑡)=𝐹^𝑇Φ_𝑚(𝑡)+Φ_𝑚(𝑡)^𝑇𝐵𝐶_𝛼+Φ_𝑀(𝑡)^𝑇𝐴+𝐸^𝑀_𝐴𝐶_𝛼𝐶_𝛼, (37) where, from Lemma??3.1 in [21], 𝐸^𝑀_𝐴→0 as 𝑀→∞. Now, by taking 𝑀=𝑚 and deleting 𝐸^𝑀_𝐴 in (37), we get (26). Obviously, if 𝑢_𝑚(𝑡) is solution of (32), then we have 𝑢_𝑚−??𝑢_𝑚→0 as 𝑚→∞.
On the other hand, from Theorem 6 we obtained 𝑢_𝑚→𝑢^∗ as 𝑚→∞. Therefore we can write 𝑢_𝑚→𝑢^∗ as 𝑚→∞ and the proof is complete.

6. Illustrative Numerical Examples

In this section, we apply our method with 𝑚=10 (BPs of degree 𝑚=10) to solve the following examples. We define 𝑦_𝑚(𝑡) and 𝑦(𝑡) for the approximate solution and the exact solution, respectively.

Example 1. Consider the nonlinear Riccati differential equation [14]: ₀𝐷^𝛼_𝑡𝑦(𝑡)=−𝑦(𝑡)²+1,0<𝑡≤1, (38) subject to the initial condition as 𝑦(0)=0. The exact solution of the equation for 𝛼=1 is given as 𝑦(𝑡)=𝑒^2𝑡−1𝑒^2𝑡.+1 (39) Numerical results compared to [14] are given in Table 1 and also Figure 1 shows the absolute error for our method for 𝛼=1 and Figure 2 shows behavior 𝑦₁₀(𝑡) for different values of 𝛼.

Table 1 

Numerical results for 𝛼=1 and 𝑚=10 in Example 1 with comparison to exact solution and [14].

Figure 1 

Behavior of the absolute error function in Example 1 for 𝛼=1 and 𝑚=10.

Figure 2 

Behavior of 𝑦10(𝑡) in Example 1 for different 𝛼 and exact solution for 𝛼=1.

Example 2. Consider the following quadratic Riccati differential equation of fractional order [14] ₀𝐷^𝛼_𝑡𝑦(𝑡)=2𝑦(𝑡)−𝑦(𝑡)²+1,0<𝑡≤1, (40) subject to the initial condition as 𝑦(0)=0. The exact solution of the equation for 𝛼=1 is given as 𝑦(𝑡)=1+√2tanh√12𝑡+2√log2−1√.2+1 (41) Numerical results compared to [14] are given in Table 2 and also, Figure 3 shows the absolute error for our method for 𝛼=1 and Figure 4 shows behavior 𝑦₁₀(𝑡) for different values of 𝛼.

Table 2 

Numerical results for 𝛼=1 and 𝑚=10 in Example 2 with comparison to exact solution and [14].

Figure 3 

Behavior of the absolute error function in Example 2 for 𝛼=1 and 𝑚=10.

Figure 4 

Behavior of 𝑦10(𝑡) in Example 2 for different 𝛼 and exact solution for 𝛼=1.

Example 3. Consider the nonlinear fractional differential equation: ₀𝐷^𝛼_𝑡𝑦(𝑡)=𝑦²(𝑡)+1,0<𝑡≤1,𝑛−1<𝛼≤𝑛, (42) with the initial conditions 𝑦^(𝑘)(0)=0,𝑘=0,1,…,𝑛−1. (43) This problem has been studied by using ADM [25], FDTM [26], and BPFs [27]. Our results with 𝛼=1.5, 𝛼=2.5 are compared to [25–27] in Tables 3 and 4. Therefore, we see that our method is very effective and obtained solutions that are in good agreement with the results in [25–27]. Also, Figure 5 shows behavior 𝑦_𝑚(𝑡) for different values of 𝛼.

Table 3 

Numerical results for 𝛼=1.5 in Example 3 with comparison to [25–27].

Table 4 

Numerical results for 𝛼=2.5 in Example 3 with comparison to [25, 26].

Figure 5 

Plot of 𝑦10(𝑡) for different 𝛼 in Example 3.

7. Conclusion

In this paper, we proposed a numerical method for solving the fractional quadratic Riccati differential equations by the operational matrices of the Bernstein polynomials. We applied operational matrix for fractional integration in the Riemann-Liouville sense. Then by using this matrix and operational matrix of product, we reduced the fractional quadratic Riccati differential equation to a system of algebraic equations that can be solved easily. Finally, examples have been simulated to demonstrate the high performance of the proposed method. We saw that the results were in good agreement with the analytical solutions and the solutions in the open literatures. Also, we observed that the solutions approach to classical solutions as the order of the fractional derivatives approaches 1, for fixed 𝑚.