A new operational matrix method to solve nonlinear fractional di ﬀ erential equations

: This study aims to propose novel Zernike wavelets and a new method based on the operational matrices for solving nonlinear fractional di ﬀ erential equations. First, non-orthogonal Zernike wavelets are introduced using the Zernike polynomials. Then, a new technique based on combining these wavelets with the block pulse functions is presented to derive the operational matrix of fractional integration and to solve nonlinear fractional di ﬀ erential equations. Moreover, an error analysis is conducted by providing required theorems. Besides, the proposed method is employed to solve a nonlinear fractional competition model of breast cancer. Finally, a parametric study is performed to consider the e ﬀ ect of fractional order on the population of healthy, cancer stem, tumour, and immune cells, as well as the excess estrogen.


Introduction
Bio-mathematical modelling of diseases has recently received significant attention from researchers to improve the process of treatments.Specifically, various types of cancer have been modelled using differential equations.Solving these equations enables researchers to predict the effect of drugs on the population of the immune system and tumour cells, as well as the interactions between them.
While integer-order time derivatives have been widely used in the aforementioned differential equations in primary studies, recently, fractional-order ones have been employed instead.To consider the memory of the biological systems and the history of their reactions, it has been frequently suggested that time fractional differential equations be employed [1,2].In this article, a competition model of breast cancer, previously proposed and numerically analysed in [3], is investigated in which integer-order derivatives have been replaced with timefractional ones.Fractional derivatives have also been used in several studies on the COVID-19 outbreak [4,5].
Moreover, several researches have been conducted on nonlinear fractional differential equations, e.g., solution of nonlinear fractional models using two kinds of fractional dual-function methods [6].Furthermore, the nonlinear fractional partial differential Schrödinger equation was studied using fractional mapping and fractional bi-function methods [7].Besides, conformable fractional derivative was used for solving the conformable fractional discrete complex cubic Ginzburg-Landau equation [8].In addition, the wick-type stochastic nonlinear fractional Schrödinger equation was investigated employing the fractional F-expansion method in the study by Wang and Wang [9].Additional bio-mathematical simulations have been presented in studies by Noeiaghdam et al. [10][11][12].
On the other hand, recently, breast cancer as a challenging area has drawn researchers' attention.This issue has led to conducting numerous researches focusing on providing accurate mathematical models describing various aspects of the disease treatment.In the study by Bitsouni and Tsilidis [13], the intricate interactions among the tBregs, CD8 + T, CD4 + T, and B cells and the effect of concentration of rituximab as a monoclonal antibody were studied through a computational approach.In the study by Marino et al. [14], a numerical technique was developed to obtain a clear decision support system for breast carcinoma surgeries.Besides, the accuracy of predicting tumour size was improved by a new scheme of integrating a machine-learning algorithm with a mathematical model describing a chemotherapy treatment of breast cancer [15].An experimental-mathematical basis, combining numerical magnetic resonance imaging data with a biophysical model, was employed in the study by Jarrett et al. [16] to simulate how locally advanced breast cancer responds to neoadjuvant therapy.A detailed mathematical modelling review can be found in the study by Simmons et al. [17] on the environmental factors in breast cancer invasion.In addition to the aforementioned research works on the mathematical modelling of breast cancer, numerous studies have been conducted to provide bio-mathematical simulations in other areas.The researchers interested in reading in this field might be referred to the previous studies [18][19][20][21].
In this research, radial Zernike polynomials are employed to create novel wavelets.The Zernike polynomials have been successfully used in various areas, including optical imaging.For instance, they have been employed as base functions of image moments as shape descriptors to classify benign and malignant breast masses [22].In addition, they have been utilized in analysing the surface of vibrating disks [23].Thus, using Zernike polynomials can enable researchers to model highly fluctuating phenomena.This article examines a nonlinear time fractional model of breast cancer proposed in Solís-Pérez et al. [3].First, the Zernike polynomials are employed to construct Zernike wavelets.Then, these wavelets are combined with block pulse functions to create a numerical method for solving nonlinear time fractional differential equations.Besides, error estimation is performed.Moreover, by employing this technique, an operational matrix of fractional integration is derived.Then, the mentioned breast cancer model is solved, and the effects of the fractional order on the response of the immune system, as well as the population of various cells, are investigated.
The article is organized as follows.In Section 2, some preliminaries and definitions regarding fractional calculus are presented.In Section 3, the properties of wavelets and block pulse functions are briefly reviewed.Moreover, Zernike wavelets are defined and proposed in Section 4. A novel scheme for solving fractional differential equations and the numerical solution of the breast cancer model are presented in Section 5. Section 6 presents the results, and a comprehensive discussion on the convergence, accuracy of the solution, computational costs, and the effects of the fractional orders on the field variables are provided.Finally, in Section 7, concluding clarifications are summarized.

Preliminaries
Definition 2.1.[24] The Riemann-Liouville fractional integral operator of order α is defined as follows: Definition 2.2.[25] The Caputo fractional derivative operator of order α  is defined as follows: where ≥ t 0 and ∈ n .Definitions 2.1 and 2.2 demonstrate the following properties: 3 A brief review of wavelets and block pulse functions

Wavelets
Wavelets, as special functions, rapidly oscillate and then vanish.These properties have enabled researchers to extensively employ them as accurate tools for function approximation in recent studies.Wavelets are constructed from a single function x Ψ( ) called the mother wavelet [26]: Note that for = a 2 , ( ) form an ortho- normal basis.

Block pulse functions (BPFs)
The M -set of BPFs on 0, 1 [ ) is defined as follows: where M is an integer and = i M 1, 2, 3,…, .The vector of block pulse functions, B t M ( ), is defined as follows: The following properties of b t i ( ) are considered: , , where ⨂ is the tensor product operator defined as follows: [ ] : For the functions 11) directly leads to the following properties (2): Kilicman and Al Zhour [27] proposed the block pulse operational matrix of the fractional integration, F α , as follows: where and: where 4 Zernike wavelets

The Zernike polynomials
The Zernike polynomials form orthogonal basis functions on the continuous unit disk.These polynomials are two dimensional and constructed from a product of two separate angular and radial parts defined in polar coordinate system as follows [28].
where Z r θ , n m ( ) is the Zernike polynomial of degree n and azimuthal order m.R r n m ( ) and mθ Θ( ) are, respectively, radial and angular parts.In addition, δ represents Kronecker delta function.Moreover, n and m are integers in a way that − n m ( )is an even element of ∪ 0 { }.The following proper- ties are considered for radial Zernike polynomials: where P is the Jacobi polynomial.Eq. ( 18) implies that for a certain order m, the set } can be separately employed to approximate a function in L 0, 1 2 [ ], whereas the union of these sets does not form an orthogonal set.

Zernike wavelets
In this article, the following non-orthogonal wavelets are defined as follows.First, the domain 0, 1 ) equidistant subdomains.Then, the wavelet ψ t i j , ( ) on the (i + 1)th subdomain is defined as follows: (19)   where Furthermore, the vector of all wavelets is defined as follows: .
It is obvious from Eq. ( 21) that t Λ( ) includes = M N 2 k wavelets whose polynomials' degree varies from 0 to − N 1.

Function approximation by Zernike wavelets
Since the proposed set of the Zernike wavelets, defined in Eq. ( 21), is non-orthogonal, a collocation scheme is employed to approximate a function f t ( ) as follows: where C̅ f is the vector of constant coefficients, , to be determined later.To calculate C̅ f , both sides of Eq. ( 22) are collocated at points Solving Eq. ( 23) as a set of linear algebraic equations results in the approximation of the function f t ( ) using Eq. ( 22).

Hybrid Zernike-block pulse wavelets (HZBW) method
In this section, the proposed Zernike wavelets are combined with block pulse functions to develop a new technique for solving fractional differential equations.HZBW method is based on approximating a function using block pulse functions and Zernike wavelets, described in Section 5.1.

Function approximation by HZBW
First, the vector of Zernike wavelets is approximated by the vector of the block pulse function: where ϕ is an × M M unknown block pulse coefficient matrix to be determined by using collocation points Applying Eq. ( 7) on the right-hand side of Eq. ( 25) leads to directly obtaining the block pulse coefficient matrix ϕ given as follows: where Λ i is the ith element of the vector Λ, which can be expressed, according to Eq. ( 21), as follows: Besides, according to Eq. ( 24), the vector of the block pulse function is approximated by the vector of Zernike wavelets as follows: In a similar way, an arbitrary function f t ( ) can be approximated by block pulse functions via a collocation scheme: Hence, on the basis of Eqs. ( 28) and ( 29), the function f t ( ) can be approximated as follows: Eq. ( 30) states that the function f t ( ) can be expressed and approximated in terms of Zernike wavelets Λ t ( ) using block pulse coefficient matrix ϕ.This proposed method is considered the HZBW approximation scheme.
This technique is used to derive the Zernike-block pulse operational matrix of fractional integration in Section 5.2.

HZBW method to solve set of nonlinear fractional differential equations (SNFDEs)
Consider the following SNFDEs: where d is the number of equations and i is a nonlinear operator representing the ith fractional differential equation.Besides, U t d ( ) demonstrates the vector of unknown functions: Moreover, i is a nonlinear function to which the fractional derivative of the order α i is applied.It is defined as follows: To solve SNFDEs Eq. ( 33), unknown functions, u t i ( ), their fractional derivatives, and nonlinear functions, i , are approximated by Eq. ( 29): where

and C D T
i are unknown constants to be determined later.It should be mentioned that without loss of generality in Eq. ( 35), the function . This decision enables us to simplify the following relations.Now, the mentioned constants C T i and C D T i can be calculated as follows.Applying I α to both sides of Eq. (36) and using Eqs.( 4), (14), and (35) lead to: Eq. (37) reveals how C D T i is dependent on and can be substituted with Besides, C i can be obtained in terms of C ui as follows.According to Eq. ( 29), .
On the other hand, U τ d j ( ) can be obtained as follows.
According to Eqs. ( 29) and (34), u τ i j ( ) equals to the jth element of the vector C ui : where C u j i ( ) is the jth element of the vector C ui .Thus, Eq. ( 40) is changed to: Hence, Eq. ( 38) is simplified to: .
1 i , Eq. ( 42) can be simply presented as follows: where X is the vector of all unknowns.In addition, sub- stituting Eq. (44) into Eq. (37)gives C D i in terms of X : Substituting Eqs. ( 34)-(36) as well as Eqs.( 44) and (45) into Eq.( 33) results in Now, all arguments of i in Eq. ( 46) have been rewritten in terms of the unknowns vector X and block pulse functions B t M ( ).Then, the left-hand side of Eq. ( 46) is approximated using Eq. ( 29): Similar to the procedure performed through Eqs.(38)-( 44), C T i can be proved to be in terms of Eq. ( 48) demonstrates a set of Md nonlinear algebraic equations to be solved using a proper numerical method to obtain unknowns C ui .After calculating C ui , a HZBW solution of Eq. ( 33) can be expressed as follows: This method can be employed to solve numerous sets of nonlinear fractional differential equations arising in biomechanics, fluid dynamics, heat transfer, and more.Especially since fractional derivatives have been widely used in bio-mathematical modelling, this technique will definitely be utilized in future research in this area.
Operational matrix to solve nonlinear fractional differential equations  5

Error analysis of HZBW method
In this section, error bounds of approximating functions and fractional operational matrix are provided.Here, the absolute and L 2 norm of a function f x ( ) throughout its domain, a b , [ ], are denoted, respectively, by The following theorem, rewritten for the functions of one variable using [29], constructs a basis for error analysis by calculating the error of interpolation of a function by + n 1 ( ) points.
Theorem 1. Suppose that ] such that [29]: ] and . Then, the absolute maximum of H t n ( ) occurs at = t a and Proof.Eq. ( 51) can be proved by employing mathematical induction on n as follows.
Base case: concedes that the proposition holds for = n 0. Induction step: Suppose that the proposition is true for = > n k 0 providing that for any real numbers a and b, 2 .Thus: , Eq. (52) leads to: In addition, it is obvious that , By considering Eqs. ( 52)-( 54), it can be written for Therefore, the proposition holds for = + n k 1. Conclusion: As the base case and the induction step have been proved true, using mathematical induction, the proposition Eq. ( 51) holds for all nonnegative integer n.
, then there exists an approximation of f t ( ) by Zernike wavelets Λ t ( ) according to Eq. ( 22) with the following error upper bound: ).
Proof.Since Λ t ( ) includes polynomials of degree at most − N 1, Eqs. ( 19)- (22) suggest that f t ( ) can be approximated . Thus, according to Theorem 1: where is the error of f t ( ) Hence, on the basis of Lemma 1, the upper bound of the error on this interval can be given as follows: Considering the maximum error over the whole interval 0, 1 Max .
( ) with the error upper bound given in Eq. (55).
( ) is a function with conditions defined in Theorem 1, then there exists an approximation of f t ( ) by block pulse functions B t M ( ) according to Eq. ( 29) with the following error upper bound: Proof.Since B t M ( ) includes polynomials of degree zero, the proof is simply obtained by considering a procedure similar to that in Theorem 1 for = N 1.
Theorem 4.There exists an approximation of Λ t ( ) by block pulse functions B t M ( ) according to Eq. ( 24) with the fol- lowing error upper bound: Proof.Λ t ( ) includes polynomials of degree at most − N 1, each of which satisfies the required conditions of f t ( ) in Theorem (3), leading to the proof of the Eq.(59).

HZBW solution to the breast cancer competition model
The proposed method is now employed to study a fractional competition model of breast cancer, previously investigated in the study by Solís-Pérez et al. [3].This model was first introduced by Abernathy et al. [31] as a system of five ordinary differential equations governing cancer stem cells (C), tumour cells (T ), healthy cells (H ), immune cells (I ), and excess estrogen (E).Then, to include memory of the biological system, integer-order fractional derivatives were replaced with Caputo fractional derivatives in [3].These SNFDEs were given as follows: Under the following initial conditions: The values of all parameters in Eq. ( 71) are presented in Table 1.
To solve Eq. ( 71), the temporal domain must be changed from T 0, [ ] to 0, 1 [ ].Thus, the time domain is normalized by the total time, T .In other words, the non-dimensional time = t* t T is defined, and the time t must be replaced with Tt*.It can be proved that: Eq. ( 73) implies that to convert Eq. ( 71) into a nondimensional form, it necessitates multiplying the fractional derivative terms by T 1 α i and rewriting all unknown functions in terms of t*.After applying these changes, for simplicity and without loss of generality, the superscript "*" is omitted.Moreover, to provide a unified notation and more simple presentation, unknown functions are renamed as follows: Hence, Eq. ( 71) is changed to: Applying I α in both sides of Eq. ( 76) and employing Eq. ( 14) yield: Comparing Eq. ( 74) with Eq. (77) leads to: Therefore, Eqs. ( 75) and ( 76) are changed to: Now, after inserting u t i ( ) and D u t α i i { ( )} from Eqs. (79) and (80) into Eq.(74) as well as employing the properties of Eq. ( 13), Eq. ( 74) is simplified to:  Operational matrix to solve nonlinear fractional differential equations  9 (81) Note that the initial values, u 0 i ( ), are substituted from Eq. (72) in Eq. ( 81).Besides, in Eq. ( 81), the vector B t M ( ) can be cancelled out.Finally, the following set of M 5 nonlinear algebraic equations is obtained to be solved for the unknown constants Cu i : Now, an HZBW solution can be obtained using Eq. ( 49):

Results and discussion
To demonstrate the capability of the proposed method for investigating a nonlinear fractional model of breast cancer, the set of equations in Eq. ( 82) is solved using Wolfram Mathematica software for = T 50 days.The initial values and fractional orders, respectively, arising in Eqs. ( 72) and (82), are given as follows [3]: First, in this section, an experimental approach to considering the convergence of the solution is employed.As the exact solution is not available, the accuracy of the method is verified by comparing the results for the special case: = α 1 i ( = i 1, 2, 3,…, 5) with the results obtained by the command "NDSolve" of the software Wolfram Mathematica using the Runge-Kutta (RK4) scheme.Then, the CPU running time, as the computational cost of the presented method, is reported.Finally, the effect of the fractional orders on the trend of breast cancer therapy is studied.

Convergence of the solution
In this section, Eq. ( 82) is solved for various values of k and N , and the results are tabulated to demonstrate the con- vergence trend.
Besides, to analyse the rate of convergence (ROC), it is calculated according to the following formula: where f stands for each of the functions u t i ( )    and N .To study the ROC with respect to k, the problem is solved for = N 7 and = k 3,4,5,6.Then, the results are tabu- lated in columns 2-6 of the aforementioned tables.Similarly, columns 7-11 demonstrate the effect of N on the convergence of the solutions.It can be realized that by increasing the order of approximation, the solutions converge with higher precision.Besides, it is evident that ≈ ROC 1 with respect to k, while it is estimated as ≈ ROC 0.8 concerning N .Furthermore, it can be inferred that by increasing t, the solu- tion demonstrates more rapid convergence due to the asymptotic behaviour of the functions C t H t T t I t , , , ( ) ( ) ( ) ( ), and E t ( ).These tables confirm that the solution is convergent with high accuracy and an acceptable convergence rate.Hence, for the next studies in this article, = k 5 and = N 6 have been considered as the proper values of k and N , respectively.

Verification of the solution
Since this problem cannot be solved exactly, there is no absolute and definite criterion to ensure the level of accuracy.However, comparing the results of the HZBW method with a numerical method, e.g., RK4, may relatively verify the obtained results.In this section, the proposed solution procedure is verified by comparing the results of the problem with integer-order derivatives ( ) with those obtained by the RK4 method using the command "NDSolve" of the software "Wolfram Mathematica." Figures 1-5 depict the results of HZBW and RK4 methods.It is evident that these solutions are in great agreement, proving the trustworthiness of the HZBW method.To provide a more precise view of these results, they are listed in Table 7.It can be deduced that the results of these two methods are closely approximated as t increases.This table proves the high accuracy of the HZBW method as a promising scheme for solving similar nonlinear fractional differential equations in future works.Operational matrix to solve nonlinear fractional differential equations  13

Computational cost
The computational cost of numerical methods for solving algebraic equations can be generally estimated analytically by calculating the number of algebraic operations.For those methods, e.g., the HZBW scheme, in which the focus is only on converting differential equations into algebraic equations, computational cost can be analysed by measuring CPU running time.
The CPU running time depends upon the programming language, hardware configuration, and the algorithms used in the computer code.Thus, comparing CPU running time between two different methods requires these factors to be the same for both methods.In this section, CPU running time for solving the problem with fractional derivatives for various values of k and N are presented in Table 8 using "Wol- fram Mathematica" to be compared with similar methods in future research works.However, this computed CPU running time cannot be compared with those presented in former studies by employing other programming languages.Table 8 indicates that this method imposes an acceptable computational cost and can be successfully employed to solve nonlinear fractional differential equations.

The effects of fractional orders
Since fractional derivatives are used to model the memory of a biological system, proper values for the fractional orders should be determined.In this section, the responses of the system to various values of the fractional orders are investigated for the case: = = = = = α α α α α α.In this study, C t ( ), H t ( ), T t ( ), I t ( ), and E t ( ) are depicted versus time in Figures 6-10 for five different values of α, respec- tively.It can be inferred from these figures that increasing the fractional orders leads to more severe reactions in the response of the system.This implies larger values of fractional orders increase C t ( ), T t ( ), and E t , ( ) whereas increasing α results in a significant decrease in H t ( ).On the other hand, changing the fractional orders totally alters I t ( ).It is evident that for = α 0.2, I t ( ) is an increasing func- tion, whereas for = α 0.4, 0.6, 0.8, 1, it forms an absolute maximum, rising for greater values of fractional orders.

Conclusion
In this article, non-orthogonal Zernike wavelets were introduced based on Zernike polynomials.Since these polynomials have demonstrated their superiority over similar functions and capability in accurate approximations, they were utilized to construct novel wavelets in this research.Besides, a new hybrid method based on combining these wavelets with block pulse functions was proposed to solve nonlinear fractional differential equations.The limitation of the presented Zernike wavelets was that they were not orthogonal.Although this issue was resolved by combining them with block pulse functions, creating orthogonal Zernike wavelets might be the objective of future works.Moreover, this technique was employed to investigate a fractional model of breast cancer treatment.Error analysis was also conducted.Furthermore, convergence analysis, verification, and computational cost of the results were provided.Finally, the effects of fractional orders on the population of healthy, cancer stem, tumour, and immune cells, as well as the excess estrogen, were studied.It was demonstrated that this method is a promising technique for obtaining an accurate solution of nonlinear fractional differential equations.Moreover, it was illustrated that considering appropriate values of fractional derivatives may remarkably affect the trend of breast cancer treatment.Operational matrix to solve nonlinear fractional differential equations  15

| | ( 5 )
where a and b are dilation and translation parameters, respectively.Continuously varying a and b creates a family of wavelets.Setting = − n and k generate a family of discrete wavelets serving as an orthogonal basis of L 2 ( ):

Table 2 :
The

Table 3 :
The values of H t ( ) and ROC at various time

Table 4 :
The values of T t ( ) and ROC at various time =

Table 5 :
The values of I t ( ) and ROC at various time

Table 6 :
The values of E t ( ) and ROC at various time

Table 7 :
Comparison of the results obtained by HZBW and RK4