# A Numerical Approach Based on Taylor Polynomials for Solving a Class of Nonlinear Differential Equations

Authors: C. Guler; S. O. Kaya
Journal: Mathematical Problems in Engineering
 (2018)

---

## Abstract

In this study, a matrix method based on Taylor polynomials and collocation points is presented for the approximate solution of a class of nonlinear differential equations, which have many applications in mathematics, physics and engineering. By means of matrix forms of the Taylor polynomials and their derivatives, the technique we have used reduces the solution of the nonlinear equation with mixed conditions to the solution of a matrix equation which corresponds to a system of nonlinear algebraic equations with the unknown Taylor coefficients. On the other hand, to illustrate the validity and applicability of the method, some numerical examples together with residual error analysis are performed and the obtained results are compared with the existing results in literature.


---

## Body

## 1. Introduction

In this study, we consider the high order differential equation with nonlinear terms in the form(1)∑k=0mPkxykx+∑p=02∑q=0pQpqxypxyqx=gxsubject to the mixed conditions;(2)∑k=0m-1akjyka+bkjykb+ckjykc=λja≤c≤b,j=0,1,…,m-1Here Pkx,Qpqx and gx are the functions defined on the interval a≤x≤b;akj,bkj,ckj,c and λj are suitable real constants; y(x) is an unknown function to be determined.In recent years, it is well known that there exists an increasing interest in the application of the models to nonlinear problems of (1)-(2) type in biology, physics and engineering. Also, the numerical solution methods of them have been developed very rapidly and intensively by many authors [1–10]. In [11], they present a new analytical technique for obtaining the analytical approximate solutions for system of Fredholm-Integral equations based on the use of residual power series method (RPSM). A new general form fractional power series is introduced in the sense of the Caputo fractional derivative in [12]. They implement a relatively new analytic iterative technique to get approximate solutions of differential algebraic equations system based on generalized Taylor series formula and the analysis proposes an analytical-numerical approach for providing solutions of a class of nonlinear fractional Klein-Gordan equation subjected to appropriate initial conditions in Caputo sense by using the Fractional Reduced Differential Transform Method (FRDTM) in [13, 14], respectively. In this study, by means of the matrix methods based on collocation points given by M. Sezer and coworkers [15–17], we developed a new numerical method to find the approximate solutions of (1) in the truncated Taylor series form(3)yx≅yNx=∑n=0Nanxn,a≤x≤bwhere an,n=0,1,…,N≥m, are unknown coefficients to be determined; yNx is a Taylor polynomial of degree N on the interval a≤x≤b.


## 2. Operational Matrix Relations

We now consider (1) and find the matrix forms of each term in the equation. We first convert the Taylor polynomial equation defined by (3) to matrix form(4)yx≅XxAwhere(5)Xx=1x⋯xN,A=a0a1⋯aNT.Besides, it is clearly seen that the recurrence relation between the matrix Xx and its derivatives Xk(x) is(6)Xkx=XxBk,k=0,1,…,mwhere [18](7)B=010⋯0002⋯0⋮⋮⋮⋱⋮000⋯N000⋯0N+1xN+1,B0=10⋯001⋯0⋮⋮⋱⋮00⋯000⋯1N+1xN+1.From (4) and (6), we have the matrix relation(8)yxk=XxkA=XxBkA,k=0,1,…,m.In addition, we can obtain the following matrix forms of the expressions y0(x)2,  y1(x)y0(x),  y1(x)2,  y2(x)y1(x),  y2(x)y0(x) and y2(x)2, by similar computations as (4)–(8) [15, 17]:(9a)y0x2=XxX-xA-(9b)y1xy0x=XxBX-xA-(9c)y1x2=XxBX-xB-A-(9d)y2xy1x=XxB2X-xB-A-(9e)y2xy0x=XxB2X-xA-(9f)y2x2=XxB2X-xB2-where(10)y0x=yx,y1x=y′x,y2x=y′′xX-x=diag⁡XxXx⋯Xx;X-xN+1xN+12B-=diag⁡BB⋯B;B-N+12xN+12B2-=diag⁡B2B2⋯B2;B2-N+12xN+12A=a0a1⋮aNN+1x1,A-=a0Aa1A⋮aNAN+12x1.Now we can define the collocation points as(11)xi=a+b-aNi,i=0,1,…,N;a≤x0<x1<⋯<xN=b.By substituting the collocation points (11) into (1), we get the system of matrix equations, for i=0,1,…,N(12)∑k=0mPkxiykxi+∑p=02∑q=0pQpqxiypxiyqxi=gxior briefly(13)∑k=0mPkYk+∑p=02∑q=0pQpqYp,q=Gwhere(14)Pk=diag⁡Pkx0Pkx1⋯PkxNQpq=diag⁡Qpqx0Qpqx1⋯QpqxNYk=ykx0ykx1⋮ykxN,Yp,q=ypx0yqx0ypx1yqx1⋮ypxNyqxN,G=gx0gx1⋮gxN.We now put the collocation points (11) into the relation (8) and obtain the matrix equation, for i=0,1,…,N,(15)ykxi=XxiBkA⇒Yk=XBkAwhere(16)Yk=ykx0ykx1⋮ykxN,X=Xx0Xx1⋮XxN=1x0…x0N1x1…x1N⋮⋮⋱⋮1xN…xNN.On the other hand, we can write the nonlinear part of (13) as(17)∑p=02∑q=0pQpqYp,q=Q00Y0,0+Q10Y1,0+Q11Y1,1+Q20Y2,0+Q21Y2,1+Q22Y2,2.By substituting the collocation points (11) into the relations (9a), (9b), (9c), (9d), (9e), and (9f), respectively, the matrices Y0,0,Y1,0,Y1,1,Y2,0,Y2,1, and Y2,2 are obtained from the following matrix forms:(18)Y0,0=X0,0∗A-Y1,0=X1,0∗A-Y1,1=X1,1∗A-Y2,0=X2,0∗A-Y2,1=X2,1∗A-Y2,2=X2,2∗A-where(19)X0,0∗=Xx0X-x0Xx1X-x1⋮XxNX-xNX1,0∗=Xx0BX-x0Xx1BX-x1⋮XxNBX-xNX1,1∗=Xx0BX-x0B-Xx1BX-x1B-⋮XxNBX-xNB-X2,0∗=Xx0B2X-x0Xx1B2X-x1⋮XxNB2X-xNX2,1∗=Xx0B2X-x0B-Xx1B2X-x1B-⋮XxNB2X-xNB-X2,2∗=Xx0B2X-x0B2-Xx1B2X-x1B2-⋮XxNB2X-xNB2-A-=a0Aa1A⋯aNAT.By substituting the matrix relations (15) and (18) into (13), we obtain the fundamental matrix equation(20)∑k=0mPkXBkA+∑p=02∑q=0pQpqXp,q∗A-=Gor briefly(21)WA+VA-=Gwhere(22)W=∑k=0mPkXBk=wij;i,j=0,1,…,NV=∑p=02∑q=0pQpqXp,q∗=vmn;m=0,1,…,Nn=0,1,…,N+12Also we can write the fundamental matrix equation (21) in the augmented matrix form(23)W;V:Gor clearly(24)W;V:G=w00w01⋯w0N;v00v01⋯v0N+12:gx0w10w11⋯w1N;v10v11⋯v1N+12:gx1⋮⋮⋱⋮⋮⋮⋱⋮⋮wN0wN1⋯wNN;vN0vN1⋯vNN+12:gxN.Besides, by using the matrix relation (8), the fundamental matrix equation of the mixed conditions (2) becomes(25)∑k=0m-1akjXa+bkjXb+ckjXcBkA=λjor briefly(26)UA+O∗A-=λ⇒U;O∗:λor clearly(27)U;O∗:λ=u00u01⋯u0N;00⋯0:λ0u10u11⋯u1N;00⋯0:λ1⋮⋮⋱⋮⋮⋮⋱⋮⋮um-1,0um-1,1⋯um-1,N;00⋯0:λm-1so that(28)U=uj0uj1⋯ujN;j=0,1,…,m-1=∑k=0m-1akjXa+bkjXb+ckjXcBkλ=λ0λ1⋯λm-1TO∗=00⋯0zeromatrix;O∗mxN+12.Consequently, to find the Taylor coefficients, an(n=0,1,…,N), related to the approximate solution (3) of problem (1)-(2), by replacing the m row matrices (26) by the last m rows (or any m rows) of the augmented matrix (23), we obtain new augmented matrix.(29)W~:V~:G~=w00w01⋯w0N;v00v01⋯v0N+12;gx0w10w11⋯w1N;v10v11⋯v1N+12;gx1⋮⋮⋱⋮⋮⋮⋱⋮⋮wN-m,0wN-m,1⋯wN-m,N;vN-m,0vN-m,1⋯vN-m,N+12;gxNu00u01⋯u0N00⋯0λ0u10u11⋯u1N00⋯0λ1⋮⋮⋱⋮⋮⋮⋱⋮⋮um-1,0um-1,0⋯um-1,N00⋯0λm-1.From this nonlinear system, that is, the matrix equation W~A+V~A-=G, the unknown coefficients, an(n=0,1,…,N), are determined; therefore the Taylor polynomial solution is obtained as(30)yNx=∑n=0Nanxn.


## 3. Accuracy of Solutions and Residual Error Estimation

We can check the accuracy of the obtained solutions as follows [18]. Since the truncated Taylor series (3) is an approximate solution of Eq. (1), then the solution yNx and its derivatives are substituted in Eq. (1), the resulting equation must be satisfied approximately; that is, for(31)x=xl∈a,b,l=0,1,2,…RNxl=∑k=0mPkxlykxl+∑p=02∑q=0pQpqxlypxlyqxl-gxl≅0or(32)RNxl≤10-kl;klis  any  positive  numbers.If max10-kl=10-k is prescribed, then the truncation limit N is increased until the difference RNxl at each of the points becomes smaller than the prescribed 10-k. Therefore, if RNxl→0 when N is sufficiently large enough, then the error decreases.On the other hand, by means of the residual function,RN(x), and the mean value of the function, RN(x), on the interval a,b, the accuracy of the solution can be controlled and the error can be estimated [19].Thus, the upper bound of the mean error,R-n, is as follows:(33)∫abRNxdx≤∫abRNxdxand(34)∫abRNxdx=b-aRNc,a≤c≤b∫abRNxdx=b-aRNc≤∫abRNxdxRNc≤∫abRNxdxb-a=R-nMoreover we use different error norms for measuring errors. These are defined as follows:(1)
L2=∑i=0nei21/2(2)
L∞=max⁡ei,0≤i≤n(3)
RMS=∑i=1nei2/n+1whereei=yxi-yNxi; also y and yN are the exact and approximate solutions of the problem, respectively [20].


## 4. Numerical Examples

The method of this study is useful in finding the solutions of a class of nonlinear equations in terms of Taylor polynomials and the accuracy. We illustrate it by the several numerical examples and perform all of them on the computer using a program written separately in MATLAB R2017b.Example 1.
Consider the second-order nonlinear differential equation(35)y′′x-2y′x+yx+y2x-y′′xyx=gxwith the initial and boundary conditions(36)y0=3,y′0=2(37)0≤x≤1where gx=2+2ex.
While the exact solution isyx=1+2ex, the proposed method is applied and the approximate solutions of (35) under conditions (36) are obtained as y2x=3+2x+x2, y3x=3+2x+x2+0.3942x3, y4x=3+2x+x2+0.3251x3+0.1072x4, y5x=3+2x+x2+0.334149x3+0.079586x4+0.022571x5 for N = 2,3,4,5, respectively.
In Tables1 and 2, we see that absolute errors and L2, L∞, and RMS errors are calculated for N=2,3,4,5, respectively.Table 1
Absolute error of Example1 for N = 2, 3, 4, 5 and h = 0.1.
x l Absolute errors e 2 e 3 e 4 e 5 0 0 0 0 0 0.1 3.41836e-04 5.23638e-05 6.01615e-06 4.98080e-07 0.2 2.80552e-03 3.48083e-04 3.31963e-05 2.24337e-06 0.3 9.71762e-03 9.25784e-04 7.15951e-05 3.92686e-06 0.4 2.36494e-02 1.57940e-03 9.86752e-05 4.72833e-06 0.5 4.74425e-02 1.83245e-03 1.05041e-04 5.66754e-06 0.6 8.42376e-02 9.09599e-04 1.22880e-04 8.24884e-06 0.7 1.37505e-01 2.29481e-03 2.57394e-04 1.01147e-05 0.8 2.11082e-01 9.25145e-03 7.21536e-04 2.60625e-06 0.9 3.09206e-01 2.18344e-02 1.87440e-03 6.66051e-05 1 4.36564e-01 4.23636e-02 4.26365e-03 2.56735e-04Table 2
L 2 , L ∞, and RMS errors of Example 1 for N=3,4,5.
N L 2-error L ∞-error R M S-error 3 4.86820e-02 4.23636e-02 2.43410e-02 4 4.72452e-03 4.26365e-03 2.11287e-03 5 2.65709e-04 2.56735e-04 1.08475e-04Example 2.
Consider the second-order nonlinear differential equation of the form(38)d2yd2x+dydx2+y2x=1-sin⁡xwith the initial conditions(39)y0=0,y′0=1.See [21].
Similarly, while the exact solution isyx=sin⁡x, the approximate solutions of (38) under conditions (39) are obtained as y2x=x, y3x=x-0.16398x3, y4x=x-0.16971x3+0.00987x4, y5x=x-0.16673x3+0.00030x4+0.00793x5 for N = 2,3,4,5, respectively.
ConsideringN=2,3,4,5, the obtained approximate solutions are compared with the exact solution in Figure 1 and the absolute errors are demonstrated in Table 3.Table 3
Absolute errors of Example2 for N=2,3,4,5 and h=0.1.
xn Absolute errors e 2 e 3 e 4 e 5 0 0 0 0 0 0.1 1.66583e-04 2.60335e-06 2.13964e-06 3.73468e-08 0.2 0.00133 1.88292e-05 1.12187e-05 1.53195e-07 0.3 0.00448 5.23333e-05 2.24296e-05 2.16761e-07 0.4 0.01058 8.69376e-05 2.71103e-05 1.79108e-07 0.5 0.02057 7.69613e-05 2.24136e-05 2.26104e-07 0.6 0.03536 6.21533e-05 2.06813e-05 6.36595e-07 0.7 0.05578 4.62827e-04 5.84302e-05 1.25213e-06 0.8 0.08264 1.31385e-03 2.04858e-04 4.68499e-07 0.9 0.11667 2.86832e-03 5.69792e-04 6.33607e-06 1 0.15853 5.45098e-03 1.31098e-03 2.90151e-05Figure 1
Comparison of the Taylor polynomial solutions and exact solution of Example2 for N=2-5.The upper bound of the mean error,R-n, of Example 2 can be calculated as R-2=7.930306e-01,R-3=4.02186e-02,R-4=1.706537e-02,R-5=4.145337e-04 as in the method given and R-n of Example 2 is illustrated in Figure 2 for N=2,3,4,5.Figure 2
The residual error functions of Example2 for N=2-5Example 3.
Consider Volterra’s Population Model, a model for population growth of a closed system introduced by Volterra.
The model is defined by the nonlinear integro-differential equation [22](40)κdudx=ux-u2x-ux∫0tuxdx,u0=0.15where κ is a prescribed parameter. In order to examine the mathematical behavior of the individual's scale population u (t), we first set(41)yt=∫0tuxdx,y′t=ut,y′′t=u′tSubstituting (41) into (40) and using κ=0.25 for numerical purpose, we find(42)y′′t-4y′t+4y′t2+4yty′t=0with the initial conditions(43)y0=0,y′0=0.15See [23].Because the problem has not exact solution, we follow the same procedure and obtain the approximate solutions of (42) under the conditions (43). For N=5, the obtained solution using our method is y5x=0.15x+0.2550x2+0.27107596x3-0.12335678x4-0.108288366x5. In Table 4, the results of the approximate solutions of the problem that is obtained with Runge-Kutta method [24], Adomian decomposition method [23], and the proposed method on x∈[0,1] are presented in Table 4.Table 4
The results of the approximate solutions of Example3 and h=0.1.
xn Runge-Kutta method Adomian decomposition method Proposed method 0 0 0 0 0.1 1.77737e-02 1.77742e-02 1.78076e-02 0.2 4.19657e-02 4.19670e-02 4.21365e-02 0.3 7.36730e-02 7.36813e-02 7.40067e-02 0.4 1.13439e-01 1.13495e-01 1.13882e-01 0.5 1.60989e-01 1.61261e-01 1.61540e-01 0.6 2.15187e-01 2.16202e-01 2.15944e-01 0.7 2.74266e-01 2.77198e-01 2.75111e-01 0.8 3.36210e-01 3.42288e-01 3.35980e-01 0.9 3.99101e-01 4.04491e-01 3.95286e-01 1 4.61346e-01 4.37632e-01 4.44430e-01Example 4.
Consider the following differential equation:(44)y′x-2yx+y2x=gx(45)gx=x-12,0≤x≤12x2-2,12<x≤1.It is solved as follows with the proposed method under boundary conditions. (46)y0=0.(a)
For0≤x≤1/2, P0x=-2, P1x=1, Q00x=1.Forg1x=x-12, N=2 and the ordering points are x0=0,x1=1/4,x2=1/2 for the interval 0≤x≤1/2.The solution is investigated in the following form:(47)yx≅∑n=01anxn.Using our proposed method, the fundamental matrix equation is obtained as(48)P1XB+P0XA+Q00X00∗A-=G.Here,(49)P1=100010001X=X0X14X12=10011411611214Q00=100010001B=010002000P0=-2000-2000-2A=a0a1a2X00∗=X0X-0X14X-14X12X-12=100000000114116141161641161641256112141214181418116X-x=diag⁡XxXxXxA=a0a1a2T,A-=a0Aa1Aa2ATG1x=g10g114g112T=191614T.The following fundamental matrix equation is obtained as(50)-210-2038-2012A+100000000114116141161641161641256112141214181418116A-=191614.100A+000000000A-=0 is obtained for the boundary condition y0=0 and using the proposed method, the desired matrix equation(51)-210100-2012A+100000000000000000112141214181418116A-=1014is constructed.
Thus, Taylor’s coefficients are found asa0=0,a1=1,a2=0 and for 0≤x≤1/2, the solution(52)yx=x,0≤x≤12is obtained.(b)
Forg1x=x2-2,  N=2 and the ordering points are x0=1/2,x1=3/4,x2=1 for the interval 1/2 <x ≤ 1.y(1/2) = 1/2 is obtained for x0=1/2 by considering yx=x. Thus, the following matrix equation is written as follows:(53)-201211214-2-10A+112141214181418116000000000111111111A-=-7412-1.Thus, by calculating Taylor’s coefficients a0=1,a1=-1,a2=0 for 1/2<x≤1 the solution(54)yx=1-x,12<x≤1is obtained.
The solution of (44) under condition (46)(55)yx=x,0≤x≤121-x,12<x≤1overlaps with the exact solution of the problem.


## 5. Conclusion

Nonlinear differential equations used in engineering, physics, mathematics or in many modelling problems are usually difficult to solve analytically. In this study, a Taylor matrix method has been given to solve a class of nonlinear differential equations in order to find a consistent approximate solution. A significant advantage of the proposed method, the Taylor coefficients of the solution can be found obviously by developing computer programs. The results related to examples have been shown in Tables 1–4 and Figures1 and 2. As seen from tables and figures, when the value ofN is increased, the numerical results rapidly improve. Residual error analysis has been also developed for the accuracy of solutions. The proposed method can be improved with new strategies to solve other high order nonlinear differential equations.



---

*Source: 8256932-2018-09-13.xml*