Perturbation-Iteration Method for First-Order Differential Equations and Systems

Abstract

The previously developed new perturbation-iteration algorithm has been applied to differential equation systems for the first time. The iteration algorithm for systems is developed first. The algorithm is tested for a single equation, coupled two equations, and coupled three equations. Solutions are compared with those of variational iteration method and numerical solutions, and a good agreement is found. The method can be applied to differential equation systems with success.

1. Introduction

Perturbation methods are classical methods which have been used over a century to obtain approximate analytical solutions. The method has been successfully applied to differential equations, integro differential equations, and algebraic equations. Many different perturbation techniques such as the method of multiple scales, the method of averaging, the renormalization method, the Lindstedt-Poincare method, the method of matched asymptotic expansions, and their variants were developed [1, 2].

The major limitation of the perturbation methods is the requirement of a small parameter. Sometimes the small parameter may also be artificially introduced into the equations. The solutions therefore have a limited range of validity. Although the solutions are valid for weakly nonlinear problems, they are not admissible usually for strongly nonlinear problems.

Several techniques have been proposed in the literature recently to obtain admissible solutions which do not require small parameter assumption. While a complete review of the attempts is beyond the scope of this work, linearized perturbation method, parameter expanding method, new time transformations as modifications of Lindstedt-Poincare method, and iteration methods can be mentioned as examples [3–14].

Recently, a class of alternative perturbation-iteration algorithms has been proposed. The fundamentals of the algorithms were outlined for the first-order differential equations by Pakdemirli et al. [15]. Several iteration algorithms can be derived by taking different number of terms in the perturbation expansions and different order of correction terms in the Taylor series expansions. The perturbation iteration algorithm is called PIA(𝑛, 𝑚) where 𝑛 represents the correction terms in the perturbation expansion and 𝑚 represents the highest order derivative term in the Taylor series. This new method has been successfully implemented to Bratu-Type equations [16]. Solutions obtained by this new method and those obtained by variational iteration method (VIM) were contrasted, and it is shown that while PIA(1, 1) algorithm usually produces identical results with VIM, higher order algorithms PIA(𝑛, 𝑚) produce better results. The new perturbation-iteration technique was applied to nonlinear heat transfer equations by Aksoy et al. [17] very recently. One of the main advantages is that the new method does not require initial assumptions or transformation of the equations to another form. Actually, the techniques were developed first for algebraic equations [18–20] and then adopted to ordinary differential equations [15–17].

In this study, the iteration algorithms for single equations are generalized to arbitrary number of first-order coupled equations. An application of the algorithm to a single equation which is a degenerate case is treated first. Then, coupled systems with two and three equations are solved. Solutions are contrasted with available other approximate solutions or numerical solutions. It is found that the method can be effectively applied to differential equation systems as well.

2. Perturbation-Iteration Algorithm PIA(1, 𝑚)

In this section, a perturbation-iteration algorithm PIA(1, 𝑚) is constructed by taking one correction term in the perturbation expansion and correction terms of 𝑚th-order derivatives in the Taylor series expansion:

Consider the following system of first-order differential equations. 𝐹_𝑘𝑢̇_𝑘,𝑢_𝑗,𝜀,𝑡=0,𝑘=1,2,…,𝐾,𝑗=1,2,…,𝐾, (1) where 𝐾 represents the number of differential equations in the system and the number of dependent variables. 𝐾=1 for a single equation. In the open form, the system of equations is 𝐹₁=𝐹₁𝑢̇₁,𝑢₁,𝑢₂,…𝑢_𝐾𝐹,𝜀,𝑡=0,₂=𝐹₂𝑢̇₂,𝑢₁,𝑢₂,…𝑢_𝐾⋮𝐹,𝜀,𝑡=0,_𝐾=𝐹_𝐾𝑢̇_𝐾,𝑢₁,𝑢₂,…𝑢_𝐾,𝜀,𝑡=0. (2) Assume an approximate solution of the system 𝑢_𝑘,𝑛+1=𝑢_𝑘,𝑛+𝜀𝑢^𝑐_𝑘,𝑛 (3) with one correction term in the perturbation expansion. The subscript 𝑛 represents the 𝑛th iteration over this approximate solution. The system can be approximated with a Taylor series expansion in the neighborhood of 𝜀=0 as 𝐹_𝑘=^𝑀_𝑚=01𝑑𝑚!𝑑𝜀^𝑚𝐹_𝑘𝜀=0𝜀^𝑚,𝑘=1,2,…,𝐾, (4) where 𝑑=𝜕𝑢̇𝑑𝜀^𝑘,𝑛+1𝜕𝜕𝜀𝜕𝑢̇_𝑘,𝑛+1+^𝐾_𝑗=1𝜕𝑢_𝑗,𝑛+1𝜕𝜕𝜀𝜕𝑢_𝑗,𝑛+1+𝜕𝜕𝜀 (5) is defined for the 𝑛+1th iterative equation 𝐹_𝑘𝑢̇_𝑘,𝑛+1,𝑢_𝑗,𝑛+1,𝜀,𝑡=0. (6) Substituting (5) into (4), one obtains an iteration equation 𝐹_𝑘=^𝑀_𝑚=01𝑢̇𝑚!^{𝑐𝑘,𝑛}𝜕𝜕𝑢̇_𝑘,𝑛+1+^𝐾_𝑗=1𝑢^𝑐_𝑗,𝑛𝜕𝜕𝑢_𝑗,𝑛+1+𝜕𝜕𝜀^𝑚𝐹_𝑘𝜀=0×𝜀^𝑚=0,𝑘=1,2…𝐾 (7) which is a first-order differential equation and can be solved for the correction terms 𝑢^𝑐_𝑘,𝑛. Then, using (3), the 𝑛+1th iteration solution can be found. Iterations are terminated after a successful approximation is obtained.

Note that for a more general algorithm, 𝑛 correction terms instead of one can be taken in expansion (3) which would then be a PIA(𝑛, 𝑚) algorithm. The algorithm can also be generalized to a differential equation system having arbitrary order of derivatives.

3. Applications

Applications of the theory developed will be outlined in this section. A first-order single equation and coupled systems with two and three equations will be treated.

Example 1. The following first-order differential equation arising in the cooling problem of a lumped system [21] (1+𝜀𝑢)𝑑𝑢(0)𝑑𝑡+𝑢=0,𝑢=1 (8) will be treated with PIA(1, 1) and PIA(1, 2) algorithms. The specific heat is assumed to be a linear function of temperature and the equation is cast into nondimensional form as outlined in [21].

(i) PIA(1, 1) Algorithm. For the equation considered, taking 𝑀=1 and 𝐾=1, (7) reduces to 𝑢̇^𝑐_1,𝑛+𝑢^𝑐_1,𝑛𝑢̇=−^1,𝑛+𝑢_1,𝑛𝜀−𝑢_1,𝑛𝑢̇_1,𝑛 (9) which is the determining iteration equation for the perturbation correction term. Assuming an initial solution, using then (9) and (3), successive iteration functions can be determined.

An initial trial function 𝑢_1,0=𝑒^−𝑡 (10) which satisfies the boundary condition is selected. Substituting this trial function into (9), solving for the correction term, and using (3), one has 𝑢_1,1=(1+𝜀)𝑒^−𝑡−𝜀𝑒^−2𝑡. (11) The successive iterations are 𝑢_1,2=𝜀1+𝜀+²2+𝜀³6𝑒^−𝑡−𝜀1+2𝜀+𝜀²𝑒^−2𝑡+32𝜀²(1+𝜀)𝑒^−3𝑡−23𝜀³𝑒^−4𝑡,𝑢_1,3=𝑒^−𝑡840840+840𝜀+420𝜀²+140𝜀³+35𝜀⁴+21𝜀⁵+7𝜀⁶+𝜀⁷−𝑒^−2𝑡𝜀366+6𝜀+3𝜀²+𝜀³²+𝑒^−3𝑡4𝜀²(1+𝜀)²6+6𝜀+3𝜀²+𝜀³−𝑒^−4𝑡3𝜀³8+20𝜀+21𝜀²+12𝜀³+3𝜀⁴+5𝑒^−5𝑡𝜀72⁴39+93𝜀+87𝜀²+29𝜀³−43𝑒^−6𝑡𝜀20⁵(1+𝜀)²+7𝑒^−7𝑡6𝜀⁶(1+𝜀)−16𝑒^−8𝑡𝜀63⁷. (12) These results are the same with the results of the variational iteration method given in [21].

(ii) PIA(1, 2) Algorithm. A higher iteration algorithm can be constructed by taking 𝑀=2. For this choice, (7) reduces to 1+𝜀𝑢_1,𝑛𝑢̇^𝑐_1,𝑛+1+𝜀𝑢̇_1,𝑛𝑢^𝑐_1,𝑛𝑢̇=−^1,𝑛+𝑢_1,𝑛𝜀−𝑢_1,𝑛𝑢̇_1,𝑛. (13) Taking the same initial trial function as given in (10), the successive iterations are 𝑢_1,1=𝑒^−𝑡+2𝜀𝑒^𝑡−1+𝜀²𝑒^𝑡−𝑒^−𝑡2(𝑒^𝑡+𝜀)²,𝑢_1,2=𝑒^−𝑡+2𝜀𝑒^𝑡−1+𝜀²𝑒^𝑡−𝑒^−𝑡2(𝑒^𝑡+𝜀)²+𝜀³24(𝑒^𝑡+𝜀)⁵𝑒^−𝑡−1+𝑒^𝑡2×−16𝑒^2𝑡+4𝑒^3𝑡−15𝑒^𝑡𝜀−10𝑒^2𝑡𝜀+𝑒^3𝑡𝜀−3𝜀²−6𝑒^𝑡𝜀²−3𝑒^2𝑡𝜀². (14) In (13), during the iterations, the resulting equation comes out to be a variable coefficient system. For obtaining the last iteration, due to complexity, 𝑢_1,0 is taken instead of 𝑢_1,1 in the coefficients of the left-hand side. The third iteration result is not given here for brevity.

The error in Table 1 is defined as %error=numericalsolution−approximatesolutionnumericalsolution×100. (15) As can be seen from Table 1, PIA(1, 2) performs better than VIM and PIA(1, 1).

Table 1 

Comparison of percentage errors of PIA(1, 2) with PIA(1, 1) and VIM for 𝜀=1.

Example 2. Two coupled stiff system will now be considered. Solutions will be obtained by PIA(1, 1) algorithm. The coupled system is [22] 𝑢̇₁=−1002𝑢₁+1000𝑢²₂,𝑢̇₂=𝑢₁−𝑢₂−𝑢²₂, (16) with the initial conditions 𝑢₁(0)=1,𝑢₂(0)=1, (17) for which exact solutions are available as 𝑢₁=𝑒^−2𝑡,𝑢₂=𝑒^−𝑡. (18) An artificial perturbation parameter is inserted as follows: 𝐹₁=𝑢̇_1,𝑛+1+1002𝑢_1,𝑛+1−𝜀1000𝑢²_2,𝑛+1𝐹=0,₂=𝑢̇_2,𝑛+1+𝑢_2,𝑛+1−𝑢_1,𝑛+1+𝜀𝑢²_2,𝑛+1=0. (19) For (19), (7) reduces to 𝜀𝑢̇^𝑐_1,𝑛+1002𝜀𝑢^𝑐_1,𝑛=−𝑢̇_1,𝑛−1002𝑢_1,𝑛+𝜀1000𝑢²_2,𝑛,𝜀𝑢̇^𝑐_2,𝑛+𝜀𝑢^𝑐_2,𝑛=−𝑢̇_2,𝑛−𝑢_2,𝑛−𝜀𝑢²_2,𝑛+𝑢_1,𝑛+1. (20) If the initial trial functions are taken as 𝑢_1,0𝑢=1,_2,0=1, (21) the successive iterations are 𝑢_1,1=500+𝑒501^−1002𝑡,𝑢501_2,11=−−𝑒501^−1002𝑡[phone]𝑒^−𝑡,𝑢1001_1,2=500−125751501500𝑒^−2004𝑡+1260031297535012006000𝑒^−1003𝑡+5020025011006009𝑒^−2𝑡−10020012006000𝑒^−𝑡502002501+𝑒^−1002𝑡−1007012008+2515032530012000𝑡,𝑢251252001_2,21=−+𝑒125751501^−2004𝑡−2521321367595031003𝑒^−1003𝑡251252001+𝑒^−1002𝑡335670670−839182520846672000𝑡251503253001+𝑒^−𝑡504264776770742984+5042647767767640032006𝑡.502002501 (22) These results are identical with the results of the variational iteration method given in [22].

Example 3. The problem of spreading of a nonfatal disease in a population which is assumed to have constant size over the period of the epidemic is considered in [23]. The following system determines the progress of the disease: 𝑢̇₁=−𝛽𝑢₁𝑢₂,𝑢̇₂=𝛽𝑢₁𝑢₂−𝛾𝑢₂,𝑢̇₃=𝛾𝑢₂,𝑢₁(0)=20,𝑢₂(0)=15,𝑢₃(0)=10. (23) The system is solved using PIA(1, 1). Perturbation parameter is artificially introduced as 𝐹₁=𝑢̇_1,𝑛+1+𝛽𝜀𝑢_1,𝑛+1𝑢_2,𝑛+1𝐹=0,₂=𝑢̇_2,𝑛+1−𝛽𝜀𝑢_1,𝑛+1𝑢_2,𝑛+1+𝛾𝜀𝑢_2,𝑛+1𝐹=0,₃=𝑢̇_3,𝑛+1−𝛾𝜀𝑢_2,𝑛+1=0. (24) For the previous equations, (7) reduces to 𝑢̇_1,𝑛+𝜀𝑢̇^𝑐_1,𝑛+𝛽𝜀𝑢_1,𝑛𝑢_2,𝑛𝑢̇=0,_2,𝑛+𝜀𝑢̇^𝑐_2,𝑛+𝜀𝛾−𝛽𝑢_1,𝑛𝑢_2,𝑛𝑢̇=0,_3,𝑛+𝜀𝑢̇^𝑐_3,𝑛−𝛾𝜀𝑢_2,𝑛=0. (25) The initial trial functions are 𝑢_1,0𝑢=20,_2,0𝑢=15,_3,0=10. (26) The following iteration results are obtained for 𝛽=1/100 and 𝛾=1/50: 𝑢_1,1𝑢=20−3𝑡,_2,1𝑢=15+2.7𝑡,_3,1𝑢=10+0.3𝑡,_1,2=20−3𝑡−0.045𝑡²+0.027𝑡³,𝑢_2,2=15+2.7𝑡+0.018𝑡²−0.027𝑡³,𝑢_3,2=10+0.3𝑡+0.027𝑡²,𝑢_1,39=20−3𝑡−2×10⁻²𝑡²+5612×10⁻⁴𝑡³−6218×10⁻⁵𝑡⁴−153095×10⁻⁷𝑡⁵−5672×10⁻⁸𝑡⁶+7297×10⁻⁸𝑡⁷,𝑢_2,3=15+27×10⁻²9𝑡+5×10⁻²𝑡²−2817×10⁻⁵𝑡³−5138×10⁻⁵𝑡⁴+153095×10⁻⁷𝑡⁵+5672×10⁻⁸𝑡⁶−7297×10⁻⁸𝑡⁷,𝑢_3,3=10+3×10⁻¹𝑡+27×10⁻³𝑡²+325×10⁻³𝑡³−272×10⁻⁵𝑡⁴,𝑢_1,49=20−3𝑡−2×10⁻²𝑡²+5612×10⁻⁴𝑡³−63638×10⁻⁶𝑡⁴−1266034×10⁻⁸𝑡⁵−976412×10⁻⁹𝑡⁶[phone]×10⁻¹¹𝑡⁷[phone]×10⁻¹³𝑡⁸−139833288956×10⁻¹⁵𝑡⁹−131800772×10⁻¹⁶𝑡¹⁰[phone]×10⁻¹⁸𝑡¹¹[phone]×10⁻¹⁵𝑡¹²−2518767952×10⁻¹⁸𝑡¹³+5904914×10⁻¹⁸𝑡¹⁴[phone]×10⁻¹⁸𝑡¹⁵,𝑢_2,4=15+27910𝑡+𝑡500²−2817×10⁻⁵𝑡³−261814×10⁻⁷𝑡⁴+1276294×10⁻⁸𝑡⁵+4473814×10⁻¹⁰𝑡⁶−993624328×10⁻¹¹𝑡⁷−102798011112×10⁻¹³𝑡⁸[phone]×10⁻¹⁵𝑡⁹[phone]×10⁻¹⁶𝑡¹⁰−9435789637×10⁻¹⁸𝑡¹¹−334611125×10⁻¹⁵𝑡¹²[phone]×10⁻¹⁸𝑡¹³+5904914×10⁻¹⁸𝑡¹⁴−177147245×10⁻¹⁸𝑡¹⁵,𝑢_3,43=10+10𝑡+27𝑡1000²+3𝑡25000³−28172×10⁻⁷𝑡⁴−5135×10⁻⁸𝑡⁵+51035×10⁻⁹𝑡⁶+81×10⁻¹⁰𝑡⁷−72928×10⁻¹⁰𝑡⁸. (27) The previous solutions are the same as those obtained from variational iteration method [23].
Functions 𝑢₁₋₃ are plotted in Figures 1, 2, and 3. The higher iterations, that is, 𝑛=4,5,6, calculated by symbolic programs are compared with the numerical solutions. As the number of iterations increase, the approximate analytical solutions converge to the numerical solutions.

Figure 1 

Comparison of iterations (𝑛=4,5,6) of PIA(1, 1) and numerical solutions for 𝑢1.

Figure 2 

Comparison of iterations (𝑛=4,5,6) of PIA(1, 1) and numerical solutions for 𝑢2.

Figure 3 

Comparison of iterations (𝑛=4,5,6) of PIA(1, 1) and numerical solutions for 𝑢3.

4. Concluding Remarks

The newly developed perturbation-iteration algorithm is applied to systems of equations for the first time. The theory is developed first and then applied to three different problems. Based on this study and on the previous work [17], one can conclude that while PIA(1, 1) algorithm produces compatible results with the VIM method, PIA(1, 2) produces better results than the PIA(1, 1) and the VIM.