Abstract: A family of sixth order and fourth order P-stable methods for solving second order initial value problems is considered. The nonlinear algebraic systems, which results on applying one of the methods in this family to a nonlinear differential system, may be solved by using a modified Newton method. The local error estimation technique based on the derivation of suitable formula pairs is used in which we compute two approximations of the solution, one with a fourth order method and the other with a sixth order method. The error estimate is then obtained by subtracting our two approximations. The methods in each pair are chosen to have certain features in common. They have the same iteration matrix and some of the function evaluations are common to both methods. Finally numerical results are presented of these methods.
INTRODUCTION
We consider the class of direct hybrid methods proposed by Cash (1981) then extended by Khiyal (1991) for solving the second order initial value problem.
(1) |
The basic methods has the form
(2) |
The off-step values are defined by
(3) |
(4) |
and
When the method (2-4) is applied to a nonlinear differential system (1), a nonlinear algebraic system must be solved at each step. This may be solved by using a modified Newton iteration scheme. The resulting iteration matrix involves J, J2 and J3, where J is an approximation for the Jacobian matrix of f with respect to y. Since matrix products are expensive, especially for large systems and any sparsity in J will be weakened in J2 and J3 and any ill-conditioning in J will be magnified in its powers, we wish to avoid the calculation of J2 and J3. Cash (1981), Thomas (1988), Khiyal (1991) and Khiyal (2007b) have derived sixth order P-stable (Lambert and Watson, 1976) methods of the form (2-4) for which the iteration matrix is a true perfect cube.
By taking β2 = 0 in (2) and choosing the remaining parameters appropriately, several authors (Cash, 1981; Chawla, 1981; Chawla and Neta, 1986; Costabile and Costabile, 1982; Khiyal, 1991; Thomas, 1988; Khiyal, 2007a) have derived fourth order methods of the form (2-4) which are P-stable.
The present consider the fourth order method of Thomas (1987), Khiyal (1991) and Khiyal (2007a) and sixth order methods of Thomas (1988), Khiyal (1991) and Khiyal (2007b).
We give necessary and sufficient conditions for there to exit
| • | Fourth order, P-stable, two-evaluations, two-step methods with iteration matrix (I-rh2J)2 |
| • | Sixth order, P-stable, three-evaluations, two-step methods with iteration matrix (I-rh2J)3 |
The method may be chosen so that the value of r is the same for both the fourth and sixth order methods. We derive some formula pairs. (The formula pair technique discussed by Khiyal (1991) and Khiyal and Thomas (2006).
MATERIALS AND METHODS
The modified Newton iteration scheme for methods of the form (2-4) is given by:
(5) |
For sixth order methods as defined by Khiyal (2007b)
(6) |
Where:
and
(7) |
To avoid the calculation of J2 and J3 Thomas (1988) and Khiyal (1991) has shown that the iteration matrix (5) may be factorized as a true perfect cube
Where:
The necessary and sufficient conditions for this are
(8a) |
(8b) |
The resulting methods are P-stable if and only if
(9) |
hold for all H. These conditions are satisfied when r is greater than or equal to the largest real root of the polynomial equation
(10) |
Note that the largest root of (10) is r = 0.6564 to four significant figure. Then the iteration scheme (5) has the form
(11) |
The Newton iteration scheme for fourth order methods of the form (2-4) with β2 = 0 is given by
(12) |
Where:
and
(13) |
A true perfect square iteration matrix is
Where:
For the fourth order methods, this is the case if and only if
(14) |
and the resulting methods are then P-stable if and only if
(15) |
Derivation of formula pairs: The three evaluations P-stable, sixth order
methods of the type given by Thomas (1988), Khiyal (1991), Khiyal and Thomas
(2006) and Khiyal (2007b) may be combined in a variable step code with appropriate
two evaluations P-stable fourth order methods of the type given by Thomas (1987),
Khiyal (1991) and Khiyal (2007a). The iteration matrix for the sixth order methods
is
for P-stability. For common iteration matrix we choose r >=R. The local error estimate is given by
Where,
SHK64A: The sixth order method is the method given by Khiyal (1991)
and Khiyal (2007b) for which
Sixth order method:
(16) |
Where:
β2 is free parameter (we choose β2 = 1).
Fourth order method
(17) |
Where:
β1 is free parameter (we choose β1 = 1).
r is greater than or equal to the largest root of the cubic Eq.
10 for both sixth and fourth order methods.
SHK64B: The sixth order method is the method given by Khiyal (1991)
and Khiyal (2007b) for which
Sixth order method
(18) |
Where:
β1 are free parameters.
Fourth order method
(19) |
Where:
β1 are free parameters.
r is greater than or equal to the largest root of the cubic Eq.
10 for both the sixth and fourth order methods.
SHK64C: Consider the sixth order method given by Khiyal (1991) and Khiyal
(2007b) for which the points
Sixth order method
(20) |
Where:
β1 are free parameters.
Fourth order method
(21) |
Where:
β1 are free parameters.
r is greater than or equal to the largest root of the cubic Eq. 10 for both the sixth and fourth order methods.
SHK64D: Consider the sixth order given by Khiyal (1991) and Khiyal (2007b)
for which
Sixth order method
(22) |
Where:
β2 are free parameters.
Fourth order method
(23) |
Where:
β1 is free parameter.
r is greater than or equal to the largest root of the cubic Eq. 10 for both sixth and fourth order methods.
SHK64E: Consider the` sixth order method given by Khiyal (1991) and
Khiyal (2007b) for which
Sixth order method
(24) |
Where:
β1 are free parameter.
Fourth order method
(25) |
Where:
β1 is free parameter.
r is greater than or equal to the largest root of the cubic Eq. 10 for both sixth and fourth order methods.
SHK64F: Consider the sixth order method given by Khiyal (1991) and Khiyal
(2007b) for which is
Sixth order method
(26) |
Where:
β2 is free parameter.
Fourth order method
(27) |
Where:
β1 is free parameter.
r is greater than or equal to the largest root of the cubic Eq. 10 for both sixth and fourth order methods.
RESULTS
We are mainly concerned with solving oscillatory stiff initial value problem. We have tried a number of explicit scalar (nonstiff) test problems of the form (1). They give similar results and so we restrict our attention to one oscillatory example.
Example 1:
This is a pure oscillation problem whose solution has maximum amplitude unity and period approximately six. We have calculated error as |Error at t = 6|
To verify that our techniques works for systems, we use a test problem a moderately stiff system of two equations
Example 2:
For this example we have deliberately introduced coupling from the stiff (linear)
equation to the nonstiff (nonlinear) equation. For this example again we have
calculated error as
In Table 1-6, we present the following statistics:
| • | No. of evaluation of the differential equation right hand side f, FCN |
| • | No. of evaluation of the Jacobian |
| • | No. of iteration overall, NIT |
| • | No. of iteration on steps where the iteration converges, NSIT |
| • | No. of steps overall, NST |
| • | No. of successful steps to complete the integration, NSST |
| • | No. of steps where the stepsize is changes, NCST |
| • | No. of failed steps, NFST |
| • | No. of LU factorization of the iteration matrix, NFAC |
| • | No. of function evaluation required on a per step basis, rather than on each iteration, NFPS |
The results are given in Table 1-6. The cost of the starting technique is not included in the tables. Before we can employ this error estimation technique, we need a special starting technique to obtain the values y1 and y2. Thus we take two steps with fourth order, two stages Implicit Runge Kutta method. The technique used is given by Khiyal and Rashid (2005). After two steps we use the formula pair technique. After word at the end of each step, we form the error estimate Len+1. To estimate the local error we use predictor corrector approach. For predictor we use the polynomial of degree five given by Khiyal (1991).
If the local error is satisfied, we use the sixth order direct hybrid method. At the each step, we form the error estimate;
Where, is the approximation for y(tn+1) obtained by the method of order 2m, m = 2 and 3. The step size for the next step is given by:
Where
| Table 1: | Methods are compared for example 1 with Tol = 10-2 |
| Table 2: | Methods are compared for example 1 with Tol = 10-4 |
| Table 3: | Methods are compared for example 1 with Tol = 10-6 |
| Table 4: | Methods are compared for example 2 with Tol = 10-2 |
| Table 5: | Methods are compared for example 2 with Tol = 10-4 |
| Table 6: | Methods are compared for example 2 with Tol = 10-6 |
CONCLUSIONS
In the study, we have derived some formula pairs, consisting of a fourth order and a sixth order direct hybrid method. The two methods in each pair have been chosen to have some features in common, so that the computational cost of using the formula pair is reduced. The formula pairs provide an estimate of the local error and this allows the stepsize to be varied so that the size of the local error is controlled.