Multistep variable methods for exact integration of perturbed stiff linear systems

1 Introduction

Perturbed harmonic oscillators are of particular interest in many areas of physics and engineering. They are also of great interest in Astrodynamics, as newtonian equations of motion can be reduced to harmonic oscillators by means of the Kustaanheimo–Stiefel [19] and Burdet–Ferrándiz [2] transformations. In addition, all natural phenomena, which can be modelled using perturbed oscillators, permit a model using perturbed differential linear equation systems.

The numerical methods used to solve these kinds of systems should preferably have the property that, if the terms of perturbation disappear from the independent variable t at an arbitrary moment, the numerical method then integrates the non-perturbed system without any discretisation error [19]. Methods with this desirable property are found in [3, 5, 10, 11, 15, 16].

In [4] the Φ-function series method is presented, for the precise integration of perturbed linear differential equations. Compared with the G-function series [16], this method has the advantage of integrating the perturbed problem without any discretisation error, under certain conditions and with only the first two terms of the series. It also maintains the desirable properties of the Scheifele methods.

The Φ-function series method is extremely precise. However, it is difficult to calculate coefficients for each specific case it is applied to, which makes it difficult to implement on a computer.

In order to solve this problem, this article describes a variable-step, variable-order multistep method (VSVO), which maintains the good properties of the Φ-function series method and incorporates an algebraic procedure to compute the coefficients of the algorithm.

This multistep method is obtained by approximating the perturbation function derivatives, which appear in the series method, using divided differences. This allows us to establish a recurrent calculation procedure to compute the coefficients of the algorithm. This recurrent procedure is based on the relationship existing between the divided differences of the perturbation function and the elementary and complete symmetric functions.

Both explicit and implicit multistep VSVO algorithms are constructed and, from these, a predictorcorrector algorithm. A computational algorithm is also designed to implement the method on a computer.

The excellent behaviour of the method can be seen when it is applied to stiff and highly oscillatory problems, comparing the numerical results obtained with those calculated by other well-known integrators and even with methods based on Scheifele’s G-functions.

2 Basic ideas and formulations

Let us consider the following IVP:

where x : R ⟶ R^m, A ∈ M(m, R) and f : R^m × R ⟶ R^m. The components of the vector perturbation field f(x(t), t) are f_i(x(t), t) with i = 1, …, m and the field is continuous, with continuous derivatives until a certain order that satisfies the conditions for existence and uniqueness of solution. This type of system is called a perturbed linear system.

Assuming that g(t) = f(x(t), t) is analytical in I with regard to t, where it is sufficient that f is analytical in its arguments. In terms of the linear operator derivation D, with respect to the variable t, (2.1) can be written as follows:

for which it is supposed that x(t) will be the only solution, in I, which can be developed in a power series.

Applying the operator (D + B) to (2.2), where B ∈ M(m, R), and noting L₂ = D² + (A + B)D + BA, the new IVP is obtained:

whose exact solution x(t) is the same as that of (2.1) and (2.2).

The idea that leads us to consider this ‘enlarged’ IVP, is that of cancelling the perturbation with the operator (D + B).

Given that g(t) is analytical in its arguments, we can write

with [4]

The solution of the IVP (2.5) is obtained by adding a specific unperturbed IVP solution with null initial conditions to the general solution of the perturbed IVP with given initial conditions. The former can be obtained by resolving the following specific IVPs:

where X_j is a real function with values in the ring M(m, R) of the squared matrices of order m, with I_m and 0 being, respectively, the unit and neutral elements of said ring. The solutions of (2.6) are the so-called Φ-functions.

3 Multistep methods

The series method described in [4] is very precise, however there is a difficulty adapting it to each specific problem. To solve this difficulty, we will proceed to describe the conversion of the series method to a multistep method similar to the SMF [[10], [11]] and EIpPC [15], VSVO [26] methods, which has the advantage of precisely integrating the perturbed problem under certain conditions.

Denoting by t_n = t_n−1 + h_n, with n = 1, 2, …, an approach to the solution x(t), in the point t = t₁, ie x₁ = x(t₁) is given by:

Let us suppose that we have calculated an approximation to the solution x(t) and its derivative x′(t) in the point t = t_n, we shall call these approximations x_n and x′n, respectively.

As

to calculate an approximation to the solution at the point t_{n + 1}, the change was made to the independent variable t = τ + t_n, becoming (3.2) in the system:

Calculating the expansion coefficients

where

The approximation to the solution in point t_n+1 = (n + 1)h, is given by

3.1 Explicit multistep method MDFSpE, for perturbed systems

In order to obtain an explicit method the derivatives of the perturbation function are substituted by divided differences.

To make a variable step explicit multistep method of p-step, series

is truncated, such that the higher order derivative is (p − 1), i.e, q_p−1 = g^(p−1)(t_n), so that:

Properly rearranging (3.7), we obtain:

Therefore it would be necessary to use divided differences (p − 1)-th order of each of the component fields of the function g, in the values t_n, …, t_n−p+1, that is g_gi[t_n, …, t_n−p+1] function with i = 1, & m.

As the divided differences of an arbitrary function g(t), satisfy the identity

with P_k(t)=t^k/k! and H_i=t_n –t_n-i [15,25].

Denoting by

the next matrix, with m × p order and H = max {H₁, …, H_p−1}, it verify the identity

where O_p×m is a matrix of order p × m, whose i-th row is

and A_p is the non singular matrix of order p:

Using a more compact notation

Truncating the expansion obtained and solving Z_p×m it results in:

replacing in (3.8) the derivatives of the components fields of the perturbed function, we can write

being p_j with j = 1, …, p the j-th. column of the matrix Dp,n×Ap−t.

Let (dij)j=1,...,pi=1,...,p be the matrix Ap−t, then

which substituted in ∑k=0p−1(Φk+1(hn+1)+Φk+2(hn+1)B)  pk+1 allow us to write

Defining

we obtain

Replacing (3.19) in (3.15), we obtain the next formula for a explicit multistep method

We introduce the next notation, we obtain the next definition.

3.2 Implicit multistep method MDFSpI, for perturbed systems

In order to obtain an implicit method of p steps, we will use the same idea as in the previous section, the ∑k=0∞Φk(t)(qk+1+Bqk) series is truncated, such that the higher order derivative is p, i.e, q_p = g^(p)(t_n), to correspond with the latest g[t_n+1, …, t_n+1−p] divided difference, so that:

rearrange (3.21), we obtain:

It will be necessary to use the divided differences of order p-th of each of the fields of components of the function of g, in the values t_n+1, …, t_n+1−p.

Given that h_n+1 = t_n+1 − t_n, it is defined

Denoting by

the next matrix, with m × (p + 1) order, and H = max {h_n+1, H₁, …, H_p−1}, it verify the identity

where O_(p+1)×m is a matrix of order (p + 1) × m, whose i-th row is

and B_p is the non singular matrix of order p + 1

Using a more compact notation:

Truncating the expansion obtained and solving for Z_(p+1)×m results in:

replacing in (3.22) the derivatives of the components fields of the perturbed function, we can write

being p_j with j = 1, …, p + 1 the j-th column of the matrix Dp,n+1×Bp−t.

Let (dij)j=1,…,(p+1)i=1,…,(p+1) be the matrix Bp−t, then

which substituted in ∑k=0p(Φk+1(hn+1)+Φk+2(hn+1)B)pk+1 allow us to write:

Defining

we obtain

Replacing (3.33) in (3.29), we obtain the next formula for a implicit multistep method:

The next notation is introduced, giving rise to the following definition.

3.3 Predictor-corrector multistep method MDFSpPC, for perturbed systems

We define the predictor-corrector method, with variable step size of p step MDFSpPC for perturbed linear systems, which has as predictor to MDFpE and as corrector to MDFSpI.

The predictor-corrector method used is like P(EC)^μE^1−t with μ = t = 1.

4 Recurrent calculus of the matrices Ap-t, Bp-t and new definition the multistep methods

The next methods described in Definitions 3.1 and 3.2, present the difficulty that the coefficients of the matrices Ap−t=(dij)j=1,…,pi=1,…,p and Bp−t=(dij)j=1,…,p+1i=1,…,p+1 are not expressed in a recurrent way, which leads to difficulties in its codification to automatize its calculation.

Once this problem is resulted, the methods will be able to be codified and will enable one to chose the step size and the number of steps, that each execution requires.

The problem is then reduced to find a recurrent formula, that allows us to calculate the elements of matricesAp−t and Bp−t.