Abstract: A family of sixth 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. We derive methods which require only three (new) function evaluations per iteration and for which the iteration matrix is a true real perfect cube. This means that at most one real matrix must be factorized at each step. Finally numerical results are presented of these methods.
INTRODUCTION
We consider an extension of the class of direct hybrid methods proposed by Cash (1981) 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
The methods are sixth order accurate if satisfy the following conditions (the order conditions are given in full in Khiyal (1991) and Thomas and Khiyal (1992)).
P-Stability (for definition, Lambert and Watson (1976))
To analyses the stability properties of the method, we apply the methods (2-4) to the scalar test equation
y″ + c2y = 0, c real. This yield
| (5) |
Where, H = ch. From condition Eq. 3, we have
| (6) |
while from Eq. 4, we have
| (7) |
Substituting the values of yn±α1 from Eq. 6 in 7 we have,
| (8) |
Substituting the values of yn±α1 and yn±α2 from Eq. 6 and 8, respectively in Eq. 5, we have
| (9) |
Where:
| (10) |
The general solution of this difference equation is
| B1 and B2 | = | Constants |
| λ1(H) and λ2(H) | = | The zeros of the stability polynomial |
| (11) |
Method (10) is unconditionally stable if |λ1(H)|≤1
and |λ2(H)|≤1 for all values of H. Making the Routh-Herwitz
transformation
Thus, the necessary and sufficient conditions for P-stability are
for all values of H. We find from (10) that τ2 = τ0 if
| (12) |
The resulting method are P-stable if and only if
| (13a) |
| (13b) |
In addition to the conditions (12), we impose the conditions as described by Khiyal (1991)
|
|
(14a) |
and
| (14b) |
Condition (14a) are assumed for simplicity while condition (14b) may be justified as follows:
We find that for a P-stable method, we must require β2(Z++Z–) ≠ 0. For suppose that β2(Z++Z–) = 0. Then, on imposing the symmetry condition Eq. 12 and condition Eq. 14a, we find fron the order condition that
Where, a1, a2, a7 and a8 are same as defined above.
On using the conditions (12), (14a) and the order conditions (a), (b), (f) and (g) above, the P-stability condition (13) become
These cannot be satisfied simultaneously. Thus for a P-stable method, β2(Z++Z–) ≠ 0. Then, from (g) and (h)
while from (e)
Now if A+–C+ = –α1 then, from (g) and (i), α1 = 0. If R+–T+ = –α2 then, from (d), α2 = 0. Thus we must have (14b). Note that if both α1 = 0. and α2 = 0, then condition (c) is violated.
With these assumption, the conditions for sixth order accuracy become
| (15) |
Where, a3, a4, a5, a6, a7, a9, a10, a11 and a12 are same as defined above.
These conditions were derived by Khiyal (1991).
Now from the order conditions (15), the necessary and sufficient conditions for P-stability become
| (16a) |
| (16b) |
for all real values of H = ch. These conditions derived by Khiyal (1991). Note that Cash (1981) considers a subset of these methods which is P-stable and does satisfy the conditions (16).
For methods of the form (2-4) applied to the scalar test equation
Where, c, v and w real, the numerical forced oscillation is in phase (Gladwell and Thomas, 1983) with its analytical counterpart whenever
| (17) |
Efficient sixth order p-stable formulae: Observe that 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. Thus on defining
| (18) |
the modified Newton iteration scheme is
| (19) |
Where:
| (20) |
and J is an approximation for the Jacobian matrix.
Since matrix products are expensive, especially for large systems and any sparsity in J will be lost in J2 and J3 and also any ill-conditioning in J will be magnified in its powers, we wish to avoid the calculation of J2 and J3.
To avoid the calculation of J2 and J3 Thomas (1987) and Khiyal (1991) has shown that the iteration matrix (20) may be factorized as a true perfect cube
The necessary and sufficient conditions for this are
| (21a) |
| (21b) |
The resulting methods are P-stable if and only if
| (22) |
hold for all H. These conditions are satisfied when r is greater than or equal to the largest real root of the polynomial equation
| (23) |
Note that the largest root of (23) is r = 0.6564 to four significant figure. Then the iteration scheme (19) has the form
| (24) |
Three-evaluation methods: In general, the sixth order direct hybrid methods require five function evaluation per iteration. Cash (1984) derived a four evaluation P-stable sixth order method but the iteration matrix is not a perfect cube. The method is
| (25) |
Voss and Serbin (1988) have changed the parameters of Cash’s (1984) method in an attempt to obtain a perfect cube iteration matrix.
| (26) |
Here, we show the number of function evaluation may be reduced further, to three per iteration. We consider the following three cases.
Case (i):
Case (ii): α1 = 0, α2 ≠ 0, yn-α2 is independent of yn+1 and in addition, either
Case (iii): α1 ≠ 0, α2 ≠ 0 and both and are independent of yn+1, so that the functions , i = 1,2, have to be evaluated once per step, not once per iteration.
We do not consider the case α1 = α2 = 0 because such methods violate the conditions for sixth order accuracy.
Case (i): Suppose that α1 ≠ 0, α2
= 0. To obtain three-evaluation methods, we suppose that the quantity
| (27) |
To avoid the evaluation of
| (a) | As α2 = 0, the evaluation of |
| (28) |
The P-stability conditions (22) are satisfied provided r is greater than or equal to the largest root of (23). The resulting three-evaluation, perfect cube, P-stable, sixth order method is;
| (29) |
β2 and α1 are free parameters and r is greater than or equal to the largest root of the cubic Eq. 23.
In general for this class of methods,
| (30) |
giving the method
| (31) |
β2 is a free parameter and r is greater than or equal to the largest real root of the cubic Eq. 23.
| (b) | The evaluation of f(tn-α2, yn-α2) can also be avoided if the points f(tn-α2, yn-α2), f(tn+α2, yn+α2) are coincident. Then the conditions for sixth order in phase with a true perfect iteration matrix are |
| (32) |
The P-stability conditions (22) are satisfied provided r is greater than or equal to the largest root of (23). The resulting three-evaluation, perfect cube, P-stable, sixth order method is
| (33) |
and r is greater than or equal to the largest root of the cubic Eq. 23.
A further saving may be made by requiring that
| (34) |
and the resulting method is
| (35) |
β2 is a free parameter and r is greater than or equal to the largest real root of the cubic Eq. 23.
Case (ii): Suppose that α1 = 0, α2 ≠ 0, the quantity yn-α2 is independent of yn+1 if and only if R– = 0, Y– = 0 and X–y″n-α1 + Z–y″n+α1 is independent of yn+1. By imposing these conditions the sixth order accuracy conditions (15) become
| (36) |
The P-stability conditions (22) are satisfied provided r is greater than or equal to the largest root of (23). To reduce the number of function evaluations required on each iteration to three, we consider two possibilities.
If y″n-α1 is identically equal to y″n, then α1 = 0, A– = 0, B– = 1, C– = 0, s– = 0, q– = 0, u– = 0, Z– = 0, (the condition Z– = 0 is required to ensure that X–y″n-α1+Z–y″n+α1 is independent of yn+1). We obtain the following sixth order, P-stable, in phase method:
| (37) |
α2 and β1 are free parameters while r is greater than or equal to the largest root of the cubic Eq. 23.
If the points (tn-α1, yn-α1), (tn+α1, yn+α1) are coincident, then α1 = 0, A+ = A–, B+ = B–, C+ = C–, s+ = s–, q+ = q–, u+ = u– and either Z–+X– = 0 or A+ = s+ = 0. However, if A+ = s+ = 0, we have q+ = q– = 0, from the order conditions (15) and the P-stability conditions are not satisfied. Thus when A+ = s+ = 0 there are no sixth order, P-stable, in Phase methods and so instead we must take Z–+X– = 0. We obtain the following sixth order, P-stable, in phase method:
| (38) |
α2 and β1 are free parameters and r is greater than or equal to the largest root of the cubic Eq. 23.
Case (iii): Suppose that α1 ≠ 0, α2 ≠ 0, To obtain three evaluation methods, we suppose that the quantities yn-α1 and yn-α2 are independent of yn+1. Then
| (39) |
The resulting methods are in phase if and only if
| (40) |
and are sixth order accurate if
| (41) |
The iteration matrix is a perfect cube if and only if
| (42) |
The P-stability conditions (22) are satisfied provided r is greater than or equal to the largest real root of (23).
Note that Chawla and Rao (1985) derived a method of this type, given by Eq. 39, 40 and 41 and by choosing
| (43) |
Where, α is a free parameter. They showed that this method is P-stable provided However, the iteration matrix for this method is not a perfect cube.
In general for this class (iii), f(tn-α1, yn-α1), f(tn-α2, yn-α2), have to be evaluated once per step. It is possible to avoid having to evaluate one of these two functions, by making particular choices of coefficients. We know that α1 ≠ 0, α2 ≠ 0. Also, note that α1 = 0, α2 = 0 cannot both equal to one, because otherwise conditions (15iii) for sixth order accuracy is not satisfied. However either α1 = 0 or α2 = 0 can be equal to 1 provided the other is not. Thus, we can avoid the evaluation of f(tn-α2, yn-α2) by taking α2 = 1 and yn-α2 ≡ yn-1. Alternatively we may avoid having to evaluate f(tn-α1, yn-α1) by taking α1 = 1 and yn-α1 ≡yn-1. Another possibility is that we may take α1 = 1/2 and require that y″n+α1 at one step is used for y″n-α1 on the next step.
We consider each of these cases separately below.
(a) Suppose that the quantity yn-α2 is identically equal to yn-1. Then
| (44) |
and from (39), (40), (41), (42) and (22), we obtain the following in phase, three-evaluation, sixth order, P-stable, perfect cube method:
| (45) |
α1 and β2 are free parameters and r is greater than or equal to the largest root of the cubic Eq. 23.
(b) Here we suppose that yn-α2 ≡ yn-1. Then
| (46) |
From Eq. 39, 40, 41, 42 and 22, we obtain the following three evaluation, sixth order, in phase, P-stable perfect cube method:
| (47) |
| Table 1: | Sixth order P-stable methods |
α2 and β1 are free parameters and r is greater than or equal to the largest root of the cubic Eq. 23.
(c) Here we suppose that y″n+α1 at one step is used for y″n-α1 on the next step. This can be achieved by taking
From (39), (40), (41),(42) and (22), we obtain the following three evaluation, sixth order, in phase, P-stable perfect cube method:
α2 and β2 are free parameters and r is greater than or equal to the largest root of the cubic Eq. 23.
We test several of the methods derived above (Table 1).
NUMERICAL 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:
y″ + sinh y = 0, y(0) = 1, y’(0) = 0.
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:
y″1 + sinh (y1 + y2)
= 0, y1(0) = 1, y’1(0) = 0
y″2 +104 y2 = 0, y2(0) = 10-8,
y’2(0) = 0
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 ||Error at t = 6||∞.
In Table 2-7, we present the following statistics:
| • | No. of evaluation of the differential Eq. right hand side f, FCN. |
| • | No. of evaluation of the Jacobian ∂f/∂y, JCB. |
| • | No. of iteration overall, NIT. |
| • | Np. 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 cost of the starting technique is not included in the tables. 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;
| Table 2: | Methods of Table 1 are compared for example 1 with Tol = 10–2 |
| Table 3: | Methods of Table 1 are compared for example 1 with Tol = 10–4 |
| Table 4: | Methods of Table 1 are compared for example 1 with Tol = 10–6 |
| Table 5: | Methods of Table 1 are compared for example 2 with Tol = 10–2 |
| Table 6: | Methods of Table 1 are compared for example 2 with Tol = 10–4 |
| Table 7: | Methods of Table 1 are compared for example 2 with Tol = 10–6 |
Where,
where Tol is the local error tolerance and h is the current step size and
The results obtained by each of the sixth order methods for example 1 and 2 with Tol = 10–2, 10–4, 10–6 are given in Table 2-7. Note that our local error estimate is really an estimate of the error in the predictor, although it is corrected value which is accepted as an approximation to y(tn+1). This predictor could not be used to advance the step as it not P-stable. Also, we would expect a P-stable method to perform in much the same way for stiff oscillatory system of example 2 as it does for the scalar problem of example 1. however, we do not achieve this for the sixth order methods when we use these predictor and interpolant to calculate the back values. More reliable interpolants and predictors are required and their development is the subject of current research.