Research Article | Open Access

Volume 2014 |Article ID 565137 | https://doi.org/10.1155/2014/565137

Ahmad Fadly Nurullah Rasedee, Mohamed bin Suleiman, Zarina Bibi Ibrahim, "Solving Nonstiff Higher Order Odes Using Variable Order Step Size Backward Difference Directly", Mathematical Problems in Engineering, vol. 2014, Article ID 565137, 10 pages, 2014. https://doi.org/10.1155/2014/565137

# Solving Nonstiff Higher Order Odes Using Variable Order Step Size Backward Difference Directly

Revised03 Jul 2014
Accepted22 Jul 2014
Published19 Aug 2014

#### Abstract

The current numerical techniques for solving a system of higher order ordinary differential equations (ODEs) directly calculate the integration coefficients at every step. Here, we propose a method to solve higher order ODEs directly by calculating the integration coefficients only once at the beginning of the integration and if required once more at the end. The formulae will be derived in terms of backward difference in a constant step size formulation. The method developed will be validated by solving some higher order ODEs directly using variable order step size. To simplify the evaluations of the integration coefficients, we find the relationship between various orders. The results presented confirmed our hypothesis.

#### 1. Introduction

In this paper, we will focus only on nonstiff ODEs of the form given where in the interval , is the order of the ODE. We impose the condition that is continuous. This implies that solution for exists.

Many science and engineering problems are in the form of higher order initial value problems (IVPs) ODEs. Authors such as Gear , Suleiman [2, 3], Hall and Suleiman , and Omar  suggested a new approach for solving higher order ODEs directly. An algorithm was designed by Suleiman  to solve stiff and nonstiff higher order ODEs directly without reducing the order of the problems to first order. He called it the direct integration method. The drawbacks to methods described by Omar , Suleiman , Gear , and Lambert  are the tedious calculations of the divided differences and recurrence relations in computing the integration coefficients.

Here, a new and efficient algorithm for solving IVPs of higher order ODEs directly using variable order step size method in its backward difference formulation (MSBD method) is developed. In contrast to the integration coefficients used in Suleiman’s  direct integration (DI) method, the code MSBD calculates the integration coefficients once at the start and once when calculating the last point. The simpler nature of the backward difference compared to the divided difference gives an added advantage to the MSBD over the DI method when formulating and coding the method. This produces more elegant error formulae. The DI method in Suleiman  has been proven to converge. We note that the DI method is formulated based on divided difference while the proposed MSBD technique is based on backward difference. Therefore, in a similar way, we can prove the convergence of the MSBD. A brief explanation of the DI code is given below.

Let be the integrating polynomial with degree ; interpolating points is denoted by and obtained through divided difference.

We define to be the -fold integral and the integration coefficients obtained are denoted as follows: For the case , where .

Let be given explicitly by which are the predicted derivatives.

Using the coefficients, the predictor and corrector are given by

#### 2. Derivation of th Order Explicit and Implicit Backward Differences Method

In this section, the integration coefficients of the backward difference formulae are derived in order to obtain a relationship between the explicit and implicit formulae. This will make the computation easier.

Integrating (1) once yields Let be the interpolating polynomial which interpolates the values ; then Next, approximating in (8) with and letting give us where Let the generating function for the coefficients be defined as follows: Substituting (12) in gives which leads to Hence, the coefficients of the backward difference formulation, , are given by To derive the th order generating function, we integrate (1) number of times and mathematical induction gives This gives the generalized explicit coefficients which are denoted by Similarly, deriving the generalized implicit formulae gives the following relationship for the generating functions: where the coefficients are given by and is the Lagrange coefficient.

#### 3. The Relationship between the Explicit and Implicit Integration Coefficients

Calculating the integration coefficients directly is time consuming if large numbers of integrations are involved. By obtaining a recursive relationship between the coefficients, we are able to obtain the implicit integration coefficient more efficiently. The relationship between the explicit and implicit coefficients is discussed below.

For first order coefficients, It can be written as By substituting into (22), we have Solving the equation above gives the recursive relationship For second order coefficient, It can be written as Substituting (24) into the equation above gives or which can be simplified as Solving the equation above gives the recursive relationship For th order coefficient, using mathematical induction, we have In similar manner to the case of the second order coefficients, we obtain the general th order coefficients as follows:

#### 4. Order and Step Size Selection

An important decision to be made is whether to accept the result of an integration step. Although the efficiency of an algorithm is affected by the strategy for selecting the order and step size, its reliability is affected by the acceptance criteria used.

The literature on variable order step size codes show numerous strategies of order and step size selection. Varying the order in a multistep method is very easy. The order depends directly on the back values stored. After the completion of any step, the order can be increased by one if none of the back values used in the previous step are discarded. The order can also be decreased to any value desired by discarding the appropriate number of back values. Through experience, it is suggested that the order strategies which are not biased to selecting a lower order when implemented in an Adams code are efficient for nonstiff problems. After considering certain criteria, we adopted the order strategies similar to those used by Shampine and Gordon .

Let be the calculated step size and the final step size. See Suleiman . In practice, we multiply by a safety factor such that , in order to give a more conservative estimate for and hence reduce the numbers of steps rejected. For our code, we take the safety factor as 0.8. Shampine and Gordon  discuss a practical approach to the result of the convergence and stability when using a variable step size. Their experience suggests that their ratio of successive step sizes will need to be restricted in order to ensure stability. For this reason and reasons of economy discussed in the section of generating the algorithm, it would appear that the use of a constant step size is desirable. Opposing this view, however, is the fact that the use of the maximum possible step size minimizes the number of steps to be taken and hence the number of derivative evaluations.

#### 5. Error Estimation

Adopting the same strategies from Hall and Watt  for higher order ODEs, the estimations for the local errors on each step of the integration are shown below.

The predictor takes the form of where the constant is independent of . In algorithm, the corrector can be written as follows: where denotes the th backward difference using for . It can be shown that where (36), can be simplified for computation to

The Milne error estimate mentions that the local truncation error (LTE) is simply the difference between the two possibilities of and giving The selection of an appropriate , for , to control the order and step size is as mentioned in Suleiman .

When reducing the order, the estimated error would be Since a different predictor, of order one less, would have been used, from the mean value theorem, we are able to write In practice, to avoid the extra derivative evaluation required for (40), the decision to reduce the order is based on the estimate which is immediately available.

To consider the increase in order, the estimated error is where the predicted difference is based on a predictor with an extra term included. In estimating , we obtain leading to which establishes the asymptotic validity of using

#### 6. Changing the Step Size

We implement the step size changing technique as mentioned in Lambert . Implementing the Adams-Bashforth methods as predictors and Adams-Moulton methods as correctors or which is commonly known as ABM methods in PECE mode in backwards difference form, we then adopt the algorithm for doubling or halving the step size from Krogh  which was derived for the ABM method.

#### 7. Last Step for Solving Higher Order ODEs

As we mentioned previously, the method we proposed calculates the integrating coefficients only once. In the case of the last step, generally, the final step size will not always fit into our current step size strategy of , and . Here, we need to calculate the coefficients once more for the last step when the step size is in the form of , . The following are the generating functions for the explicit coefficients: and implicit coefficients: The following yields the relationship between the explicit and implicit coefficients: Tables 1 and 2 contain values for the coefficients when .

 0 1 2 3 4 5 6
 0 1 2 3 4 5 6

#### 8. The MSBD Algorithm

The previous write-up was devoted to the calculation of the integration coefficients for the explicit formulae which can form the basis of the predicted values and also the implicit coefficients which are obtained from the explicit formulae.

For the sake of clarity, we present the algorithm for the MSBD.

Step 1. The integration coefficients are calculated from the algorithm in (16), (18), and (34).

Step 2. Use the back values to obtain the predictor in (30).

Step 3. Calculate the corrected values in terms of the predicted values given in (35).

Step 4. The errors , and are obtained in (39), (42), and (46).

Step 5. Determine whether satisfies the local accuracy requirements which we take to be TOL, where with give the absolute error test; which gives a relative error test, while give a mixed error test.

Step 6. The order is lowered by one if, for , and . If we have available, we lower the order if and . The order is raised by one only after successful step at constant step size such that, for and, for . We restrict in the range .

Step  7. Step size selection is also determined based on the local accuracy requirements. The step size algorithm for doubling and halving is derived by Krogh .

Step 8. If , where is the current step size and denotes the end of the interval, repeat Steps 2–7. If not, determine where . We then calculate the new integration coefficients as given in (47), (48), and (50). Using the new coefficients, repeat Steps 2–7 and exit the program.

#### 9. Test Problems and Numerical Results

##### 9.1. Numerical Results

Tables 3, 4, 5, and 6 show the numerical results for Problems 14, when solved with direct integration, backwards difference, and also reducing to first order systems. Problems 1 and 2 are well-behaved problems whereas Problems 3 and 4 are without exact solutions. For the first two problems, we evaluate the maximum and average values of the error in the computed solution . The definition of the error is as defined above. The following notation will indicateTIME: the execution time taken in microseconds,FS: the number of failed steps,STEPS: total steps,MTD: the name of the method,SUCSTEP: the number of successful step,MAXE: the maximum error,AVER: the average error,TOL: the tolerance used,DI: direct integration,D1: direct integration with the reduction to first order system,1PBDVSO: backward difference,B1: backward difference with the reduction to first order system.The errors calculated are defined as with as the th component of . correspond to the absolute error test, correspond to the mixed error test, and correspond to the relative error test. The mixed error test is used for all test problems. The maximum and average errors are given by where is the number of equations in the systems.

 TOL MTD STEPS FS MAXE AVER TIME (total) 10−2 D1B1DI1PBDVSO 113897671 5310 8.80309 (−2)2.07488 (−1)8.78700 (−2)1.17774 (−1) 1.19400 (−1)6.87734 (−2)2.56159 (−2)1.70700 (−2) 2078 1628 969 969 10−3 D1B1DI1PBDVSO 1461118579 7021 2.68984 (−2)3.05138 (−2)2.12907 (−2)1.81903 (−2) 5.25240 (−2)6.28779 (−3)9.00022 (−3)6.05764 (−3) 2560 2071 1135 1063 10−4 D1B1DI1PBDVSO 15117094149 2110 4.30130 (−3)5.37539 (−4)4.25614 (−3)1.97342 (−5) 5.45887 (−3)1.06165 (−4)1.53472 (−3)5.05592 (−6) 3024 3169 1284 1951 10−5 D1B1DI1PBDVSO 218200161160 1700 1.84329 (−5)4.40121 (−4)6.68588 (−4)1.83864 (−5) 1.33438 (−5)1.37204 (−4)2.70629 (−4)6.47215 (−6) 3860 3652 2094 2080 10−6 D1B1DI1PBDVSO 275121179176 3010 1.80685 (−6)1.38184 (−5)3.15166 (−4)5.88468 (−6) 1.81711 (−6)2.65555 (−6)1.29162 (−4)2.26518 (−6) 52274077 2294 2272 10−7 D1B1DI1PBDVSO 219330193197 9010 4.47927 (−6)1.06409 (−7)1.44784 (−4)2.58941 (−7) 7.17250 (−6)6.01321 (−8)5.87351 (−5)9.00261 (−8) 5413 6257 2502 2513 10−8 D1B1DI1PBDVSO 348287205209 2200 3.49871 (−8)4.59227 (−7)6.69118 (−5)6.96499 (−8) 2.75009 (−8)1.19307 (−7)2.76364 (−5)2.35911 (−8) 6310 5367 8665 2664 10−9 D1B1DI1PBDVSO 260424226225 7000 2.68246 (−8)4.48526 (−10)3.11980 (−5)1.87994 (−9) 1.79939 (−8)1.44218 (−10)1.25997 (−5)4.25589 (−10) 8633 8271 29442851 10−10 D1B1DI1PBDVSO 526475370248 91000 3.45300 (−9)7.88269 (−9)1.44822 (−7)4.13390 (−9) 2.61398 (−9)1.39352 (−9)5.81042 (−8)8.63144 (−10) 9866 8822 4639 3109
 TOL MTD STEPS FS MAXE AVER TIME (total) 10−2 D1B1DI1PBDVSO 13513496173 0000 4.83764 (−2)2.63945 (−2)4.99681 (−3)1.01436 (−4) 3.21705 (−2)2.14665 (−2)2.47698 (−3)8.58874 (−5) 6621 4358 1961 3218 10−3 D1B1DI1PBDVSO 175121157162 1000 8.81849 (−3)3.90401 (−2)6.91170 (−5)2.40484 (−3) 5.48179 (−3)1.65274 (−2)4.87625 (−5)1.71894 (−3) 8294 3896 2890 3005 10−4 D1B1DI1PBDVSO 233173141142 2000 2.85719 (−3)1.44401 (−3)1.79264 (−4)2.60404 (−4) 1.39202 (−3)9.91992 (−4)1.19952 (−4)1.89056 (−4) 11684 5520 2622 2635 10−5 D1B1DI1PBDVSO 226263238237 1200 3.05429 (−4)1.97450 (−4)5.14212 (−6)1.03862 (−5) 1.22401 (−4)5.02997 (−5)3.89969 (−6)8.18350 (−6) 11095 8492 4137 4268 10−6 D1B1DI1PBDVSO 337262214218 1000 1.12897 (−5)6.52399 (−5)1.18826 (−5)1.56415 (−6) 5.53620 (−6)2.34931 (−5)9.01852 (−6)1.09680 (−6) 160608429 3766 3941 10−7 D1B1DI1PBDVSO 421337368370 2000 4.00671 (−7)1.70966 (−6)1.79221 (−7)2.75078 (−7) 4.00671 (−7)1.05119 (−6)1.35202 (−7)1.91969 (−7) 19735 10899 6213 6584 10−8 D1B1DI1PBDVSO 610523335336 4310 3.05239 (−7)4.90321 (−8)3.11508 (−8)3.64286 (−8) 2.22380 (−7)3.52515 (−8)2.58590 (−8)1.21667 (−8) 29060 16840 5641 5996 10−9 D1B1DI1PBDVSO 534555578575 2510 4.96195 (−8)9.16995 (−8)1.20539 (−8)3.80099 (−9) 1.71145 (−8)5.35093 (−8)9.25523 (−9)2.95718 (−9) 25234 17610 9347 10114 10−10 D1B1DI1PBDVSO 800671516517 2400 2.12904 (−10)8.63683 (−10)8.37294 (−10)3.30834 (−9) 1.21736 (−10)3.89440 (−10)6.51148 (−10)2.52921 (−9) 38235 21617 8541 9011
 TOL MTD (1) 10−2 1PBDVSODIB1D1 1.9788520539 + 0001.9163232164 + 0001.8780198255 + 000−1.#IND000000 + 000 10−4 1PBDVSODIB1D1 1.8706114418 + 0001.8694497531 + 0001.8694030590 + 000−1.#IND000000 + 000 10−6 1PBDVSODIB1D1 1.8694735377 + 0001.8694497531 + 0001.8694368915 + 000−1.#IND000000 + 000 10−8 1PBDVSODIB1D1 1.8694389431 + 0001.8694388843 + 0001.8694402726 + 000−1.#IND000000 + 000 10−10 1PBDVSODIB1D1 1.8694388662 + 0001.8694388834 + 000Fail to compute−1.#IND000000 + 000
 TOL MTD (1) 10−4 1PBDVSODIB1D1 −1.5171246226  −  009−2.5152953310  −  0089.9810188187  −  009Fail to compute 10−6 1PBDVSODIB1D1 1.0035458019  −  0081.0609214996  −  0089.9998786959  −  009Fail to compute 10−8 1PBDVSODIB1D1 1.0004172282  −  0081.0004031492  −  0081.0000000079  −  008Fail to compute 10−10 1PBDVSODIB1D1 9.9999571459  −  0099.9999855296  −  0099.9999999555  −  009Fail to compute

The result of the backwards difference method (1PBDVSO) is measured against D1, B1, and DI method for solving the same problems. The variable step size order codes are implemented for solving the problems. We begin to compute the solution using a higher tolerance ().

Problem 1. Consider Solution is as follows: First order system is as follows: Solution is as follows: For source, see Shampine and Gordon .

Problem 2. Consider Solution is as follows: First order system is as follows: Solution is as follows: For source, see Suleiman .

Problem 3. Van der Pol’s equation is as follows: First order system is as follows: For source, see Suleiman .

Problem 4. Control theory is as follows: