Chapter 7.  Numerical Methods for Initial Value Problems

Table of Contents
7.1. Numerical Solution of Initial Value Problems
7.1.1. Explicit Euler Method (Forward Euler Method)
7.1.2. Numerical Stability and Stability Domain
7.1.3. Modified Euler Method (Trapezoidal Rule)
7.1.4. Implicit Euler method
7.1.5. Summary of Euler’s Method
7.1.6. General One-Step Procedures
7.1.6.1. Classical Runge-Kutta Methods
7.1.6.2. General Runga-Kutta Methods
7.1.6.3. Stability Area of Runge-Kutta Methods of Order 1≤p≤4
7.1.7. Step Size Control
7.1.7.1. 3. Order Runge-Kutta, Two Calculations of an Integration Step
7.1.7.2. Embedded Methods
7.1.8. Linear Multi-Step Methods
7.1.9. Activation of Linear Multi-Step Procedures
7.1.10. System of Differential Equations
7.1.11. BDF Methods
7.1.12. Remarks on Stiff Differential Equations
7.1.13. Implicit Runge-Kutta Methods
7.1.14. Comparison of Methods for Numerical Solution of Initial Value Problems (IVP)

7.1.  Numerical Solution of Initial Value Problems

The solutions for the ordinary differential equations one has been dealing with till now must, at least for the non-linear problems, be solved by a numerical method solution using a digital computer.

In order to derive appropriate methods, one shall consider the following numerical method used to solve a scalar differential equation of this kind

 

x ˙ = f ( x , t )

(7.1)

with the initial condition

 

x ( 0 ) = x 0 .

(7.2)

The methods derived in the following sections can be extended to a differential equation of higher order by means of the previously introduced substitution method. The application of this method on systems of differential equation is also unproblematic.

The solution of initial value problems, in numerical methods, allow for the determination of solutions x ( t n ) for a series of discrete points in time (grid points) t n with

 

t n = t n 1 + h n .

(7.3)

where h n is the time increment which can generally change for every step. For the sake of simplicity , one shall refer to it as constant time increment h from now on.

7.1.1.  Explicit Euler Method (Forward Euler Method)

One obtains the Forward Euler Method by substituting the differential term on the left-hand side of (7.1) with a difference term:

 

x ˙ n x n + 1 x n h

(7.4)

Forward Euler Method.
Figure 7.1. Forward Euler Method.


Solving for x n + 1 , this yields an approximate formula for the value of x at position t n

 

x n + 1 = x n + h f ( x n , t n )

(7.5)

Thus, the solutions can be recursively calculated as follows:

 

x 1 = x 0 + h f ( x 0 , t 0 ) x 2 = x 1 + h f ( x 1 , t 1 ) x n = x n 1 + h f ( x n 1 , t n 1 )

(7.6)

Hence, an easily applicable and understandable method is made available. But the practical application of this method faces various difficulties:

  1. The error which emerges within the process of discretization (substitution of the differential by the difference quotient) strongly depends on the time increments h .

  2. In some cases one can observe a trend of numerical instability. An analysis of this effect is based on the so-called test equation

 

x ˙ = α x with x ( 0 ) = x 0 < 0, α and α < 0.

(7.7)

The exact solution

 

x = x 0 e α t ,

(7.8)

is always positive and converges with t to zero. The least prerequisite for a useful approximation method is that it also converges to zero and always remains positive.

If we apply the forward Euler Law on equation (7.7), we will obtain a solution at position t n

 

x n + 1 = ( 1 α ) x n

(7.9)

which means that the approximation solution remains positive and it only disappears where t converges to infinity, if it satisfies that α h < 2 .

If it yields 2 < α h < 1 , the approximation solution indeed disappears where t converges to infinity. The approximation solution oscillates around the zero point value. The solution increases with each integration section and changes the algebraic sign at the same time where α h < 2 . One only obtains a meaningful numerical solution when h < 1 α . This characteristic of the Euler method is referred to as numerical instability.

Example 7.1: Numerical instability of the Euler method

The solution of the following equation

x ˙ = 10 x

should numerically be calculated.

The application of the Euler method in the interval [ 0, 16 ] with various increments reveals the results represented in [Figure 7.2]. Obviously, the error drastically increases if we double the increment.

Furthermore, the solution oscillates for h = 0.16 and it finally diverges for h = 0.32 [Figure 7.3].

Application of forward Euler method.
Figure 7.2. Application of forward Euler method.


Forward Euler method; instability.
Figure 7.3. Forward Euler method; instability.


In order to analyze the characteristics of oscillating solutions, one frequently allows α to be specified as a complex number, see the following section.

7.1.2.  Numerical Stability and Stability Domain

For the stability analysis of a numerical integration method the linear test equation

 

x ˙ = λ x ,   λ ,   x ( t = 0 ) = x 0 < 0

(7.10)

is applied. The exact solution of this initial value problem

 

x = x 0 e λ t

(7.11)

always remains positive and converges towards zero for t . If one tries to solve the test equation by an integration method, the least requirement to obtain a useful approximate solution is that it also converges towards zero and remains ever positive.

Definition 7.1

A numerical integration method is called (numerically) stable in a domain A = { α h : α  komplex mit  ( α ) , h } of the complex number field, if the series { x n } of the approximate solutions decreases absolutely at the points of time t n for t n (according to the behaviour of the exact solution).

Let’s analyze the explicit Euler’s method. By inserting the test equation (7.10) in (7.5) we obtain

 

x n + 1 = ( 1 + λ h ) x n = ( 1 + λ h ) n + 1 x 0 .

(7.12)

According to the stability condition this method is numerically stable if

 

lim n | x n + 1 | = 0.

(7.13)

This is only valid for

 

| 1 + λ h | R ( z ) = R ( λ h ) < 1.

(7.14)

I.e., ( 1 + z ) must be inside a circle of radius 1.0 (the so-called unit circle) around the origin. Thus, z must be inside the circle of radius 1.0 around the point ( 1, 0 ) , such that the explicit Euler’s method remains numerically stable. The stability domain of the explicit Euler’s method is depicted in [Figure 7.4]. The explicit Euler’s method is neither A-stable nor F-stable.

Stability domain of the explicit Euler’s method.
Figure 7.4. Stability domain of the explicit Euler’s method.


Definition 7.2

A numerical integration method for differential equations is A-stable or Absolute-stable if and only if its stability domain covers at least the complete left z-half plane.

Definition 7.3

A numerical integration method for differential equations is F-stable or faithfully-stable if and only if its stability domain is exclusively the left z-half plane.

7.1.3.  Modified Euler Method (Trapezoidal Rule)

The modified Euler Method results from the application of the Trapezoidal Rule on the integration of the function f ( x , t ) .

Derivation of the Trapezoidal Rule.
Figure 7.5. Derivation of the Trapezoidal Rule.


Trapezoidal Rule:

 

a b f ( x , t ) d t b a 2 ( f ( a ) + f ( b ) )

(7.15)

applied on f ( x , t ) , we obtain with t n = a , t n = b and h = a b :

 

x n + 1 x n = t n t n + 1 f ( t , x ) d t h 2 ( f ( x n + 1 , t n + 1 ) + f ( x n , t n ) )

(7.16)

and therefore as a calculation rule:

 

x n + 1 = x n + h 2 ( f ( x n + 1 , t n + 1 ) + f ( x n , t n ) )

(7.17)

In contrast to the Forward Euler Method, in this case the resulting calculation rule is not explicitly solved for the searched quantity x n + 1 . Thus this method is said to be an implicit integration method. The question emerges, how to solve (7.17) for x n + 1 in practice.

Case 1: f is linear in x

Obviously, (7.17) can easily be solved for x n + 1 .

Example 7.2:

x ˙ = a x + cos ( t )

We obtain a calculation rule:

x n + 1 = x n + h 2 ( a x n + 1 + cos ( t n + 1 ) + a x n + cos ( t n ) )

This special case provides the option to explicitly solve for x n + 1 . We obtain the following formula

x n + 1 ( k ) = x n + h 2 [ f ( x n + 1 ( k 1 ) , t n + 1 ) + f ( x n , t n ) ] .

Case 2: f is a non-linear function of x.

It is only possible in special cases to explicitly solve for x n + 1 . Instead, we must generally resort to numerical methods for the solution of non-linear algebraic equation systems. In this case, we have to have a good initial value for x n for an iterative solution. Thus, we apply an iterative method. One possibility is the method of successive substitution and fixed-point iteration respectively (Section [Section 9.3]):

 

x n + 1 ( k ) = x n + h 2 [ f ( x n + 1 ( k 1 ) , t n + 1 ) + f ( x n , t n ) ]

(7.18)

where x n + 1 ( k ) identifies the value of k t h iteration where x n + 1 and x n + 1 ( 0 ) represents an appropriate initial value for iteration. We chose x n , as an initial value. Thus, the first iteration increment is identical with the formula of the explicit Euler Rule. As soon as the difference of the successive two iterations is | x n + 1 ( k ) x n + 1 ( k 1 ) | ´ < ε , we cancel the iteration. Here, ε is a given accuracy bound which should be chosen near to the machine accuracy. However, the stability domain decreases by applying the fixed-point iteration, as depicted in [Figure 7.6]. For this reason, the Newton-Raphson iteration method (Section [Section 9.3]) is applied alternatively as a rule.

Stability domain after applying the fixed-point iteration.
Figure 7.6. Stability domain after applying the fixed-point iteration.


To sum up, we obtain the following calculation rule:

 

x n + 1 ( 0 ) = x n x n + 1 ( 1 ) = x n + h 2 [ f ( x n + 1 ( 0 ) , t n + 1 ) + f ( x n , t n ) ] (forward Euler method) x n + 1 ( 2 ) = x n + h 2 [ f ( x n + 1 ( 1 ) , t n + 1 ) + f ( x n , t n ) ] ( Runge-Kutta 2 nd order )

(7.19)

etc. until convergence arrives.

The question is why the modified Euler method has a higher accuracy and an improved stability behavior among numerical integration methods.

In order to answer that question, we have to consult the test equation x ˙ = α x . In this case, we obtain the following calculation rule

 

x n + 1 = x n α h 2 ( x n + 1 + x n ) x n + 1 = 1 α h 2 1 + α h 2 x n

(7.20)

If we develop the coefficient of the factor in front of x n into a Taylor series, we obtain

 

( 1 α h 2 ) 1 1 + α h 2 = ( 1 α h 2 ) ( 1 α h 2 + ( α h ) 2 4 geometric series ) = 1 α h + 1 2 ( α h ) 2 1 4 ( α h ) 3 +

(7.21)

If we apply the same analysis for the forward Euler-method, we obtain

 

1 α h

(7.22)

as represented above. If we compare the same results (7.21) and (7.22) with the exact result

 

e α h = 1 α h + 1 2 ( α h ) 2 1 6 ( α h ) 3 +

(7.23)

we obtain the following cancellation error for the modified Euler formula

 

1 12 ( α h ) 3 +

(7.24)

and for the explicit Euler formula as

 

1 2 ( α h ) 2 + ...

(7.25)

Thus, the Euler formula precisely describes the exact solution locally up to two elements of second order (method of second order). In contrast to that, the forward Euler method only represents an approximation up to the first element of first order (method of first order).

A further important advantage of the Trapezoidal Rule is its stability behavior. One obtains the following equation for n t h integration section:

 

x n + 1 = ( 1 α h 2 1 + α h 2 ) n + 1 x 0

(7.26)

and the absolute value of the approximation solution exactly coincides where Re ( α ) < 0 , if it yields

 

| 1 α h 2 | < | 1 + α h 2 |

(7.27)

But this is the case where all Re ( α ) < 0 , i.e. the trapezoidal rule is stable for arbitrary step sizes. The trapezoidal rule is A -stable as well as F -stable.

7.1.4.  Implicit Euler method

We obtain the implicit Euler method by substituting the forward difference quotient by the backward quotient in the explicit Euler’s process.

The implicit Euler method.
Figure 7.7. The implicit Euler method.


The application of this method too needs an iteration to calculate x n + 1 .

7.1.5.  Summary of Euler’s Method

  1. The explicit Euler’s method represents a favored method which is very easily applicable and easy to program, especially in real-time operation. The local discretization error behaves like h 2 ( O ( h 2 ) ) where h 0 , the global error behaves like h . A great disadvantage is the instability, which can, as the case may be, be eliminated by the choice of extreme short increments.

  2. The modified Euler method (Trapezoidal Rule) is A -stable, the local discretization error behaves like O ( h 3 ) , the global error like O ( h 2 ) . A disadvantage (which is similar to all implicit methods) is the necessity to solve a non-linear equation system within each integration section.

  3. In respect to the discretization error, the implicit Euler method behaves like the explicit method, but in contrast to the explicit version, the implicit one has the advantage of being A -stable.

7.1.6. General One-Step Procedures

The Euler methods discussed above are the simplest methods of the great class of one-step procedures. Here, approximations x n + 1 will only be calculated at t n + 1 = t n + h n solely from the approximation x n at the points t n and the increments h n .

One-step procedures can generally be written on the form

 

x n + 1 = x n + h n Φ ( x n , t n , h n )

(7.28)

with the so-called process regulation Φ .

Example 7.3

For the explicit Euler method it is true:

Φ ( x n , t n , h n ) = f ( x , t ) .

To measure the quality of one-step procedures, we utilize the following terms:

Initial value problem

x ˙ = f ( x , t ) , x ( t 0 ) = x 0 .

This means that

ε ( h ) : = x ( t n + h ) ( x n + h Φ ( x n , t n , h n ) )

is the local error which emerges in this step of the method defined in (7.28).

This procedure is consistent, if it yields

ε ( h ) = O ( h ) ,

it is of p t h order where

ε ( h ) = O ( h p + 1 ) .

Annotations 7.1:

If the approximation x n of x n ( t n ) with the one-step procedure

x n + 1 = x n + h n Φ ( x n , t n , h n ) , n = 0, , N

is calculated, Φ satisfies a Lipschitz condition in respect to x [ a , b ] ×

| Φ ( x , t , h ) Φ ( y , t , h ) | L | x y |

and the one-step procedure is consistent of p + 1 order

| x ( t + h ) ( x ( t ) + h Φ ( x n , t n , h n ) ) | c h p + 1 ,

the global error is true for

| δ n + 1 | = | x n + 1 x ( t n ) | c ( b a ) e h ( b a ) h p

with h : = max j = 0,1,..., n + 1 h j .

Expressed in words this means that the consistency order of the global error is always one order lower than the local error.

7.1.6.1. Classical Runge-Kutta Methods

A basic disadvantage of Euler’s method is the low accuracy achieved. This demands very short integration increments h and leads to high computing times and an accumulation of round-off errors during the calculation.

Very early endeavors have been made to increase the degree of accuracy. An option is to calculate the right-hand side of differential equations at additional interpolation points.

We consider again the differential equation (7.1). If x n is given, we can integrate (7.1) in the interval [ t n , t n + 1 ] in order to calculate the function value x n + 1 where t n + 1 = t n + h :

 

x n + 1 = x n + t n t n + 1 f ( x , t ) d t

(7.29)

We obtain the Runga-Kutta-method if we approximately integrate the right-hand side of (7.29)

Runga-Kutta Method of Second Order

The following statement results from the application of the Trapezoidal Rule in the approximation integration of the integral in (7.29)

 

t n t n + 1 f ( x , t ) d t h 2 [ f ( x n , t n ) + f ( x n + 1 , t n + 1 ) ]

(7.30)

The value for x n + 1 is unknown, thus the term f ( x n + 1 , t n + 1 ) will be approximated by the explicit Euler method. Thus, we obtain the following formula

 

x n + 1 = x n + h 2 [ f ( x n , t n ) + f ( x n + h f ( x n , t n ) ) ]

(7.31)

which will usually be formulated in the following calculation formula:

 

k 1 = h f ( x n , t n ) , k 2 = h f ( x n + k 1 , t n + 1 ) , x n + 1 = x n + 1 2 [ k 1 + k 2 ] .

(7.32)

The Runga-Kutta method at hand is also referred to as a predictor-corrector method on the basis of the Euler method. In this context, the explicit Euler method plays the role of the predictor whereas the Trapezoidal Rule inherits the role of the corrector.

To determine the accuracy, with which this method discretizes the differential equation, we develop x in the surrounding of t n with the Taylor series

 

x n + 1 = x n + h f + h 2 2 [ f t + f x f ] + h 3 6 [ f t t + 2 f t x f + f t t f 2 + f t f x + f x 2 f ] + O ( h 4 )

(7.33)

at the same time, the partial derivatives of f will be shortened by:

f t = [ f t ] t = t n x = x n , f x = [ f x ] t = t n x = x n , etc.

For means of comparison, we develop (7.31) in an appropriate Taylor series

 

x n + 1 = x n + h f + h 2 2 [ f t + f x f ] + h 3 4 [ f t t + 2 f t x f + f x x f 2 ] + O ( h 4 )

(7.34)

By comparing (7.33) with (7.34), it becomes clear that the error which occurs with every integration section is proportional to h 3 .

Hereafter, the remaining Runga-Kutta methods will only be specified by the error order and the calculation formula without derivation.

Runge-Kutta methods of third order

Calculation formula:

 

k 1 = h f ( x n , t n ) k 2 = h f ( x n + k 1 2 , t n + h 2 ) k 3 = h f ( x n k 1 + 2 k 2 , t n + h ) x n + 1 = x n + 1 6 ( k 1 + 4 k 2 + k 3 )

(7.35)

Local discretization error: O ( h 4 ) .

Runge-Kutta Method of Fourth Order

This method is based on Simpson’s 1 / 8 -Rule and yields

 

k 1 = h f ( x n , t n ) k 2 = h f ( x n + k 1 2 , t n + h 2 ) k 3 = h f ( x n + k 2 2 , t n + h 2 ) k 4 = h f ( x n + k 3 , t n + h ) x n + 1 = x n + 1 6 ( k 1 + 2 k 2 + 2 k 3 + k 4 )

(7.36)

The local error order of this method is O ( h 5 ) .

This version of the Runge-Kutta method is also referred to as classical Runge-Kutta method:

Stability Observation

In order to analyze the Runga-Kutta method, we consider again the test equation with a real α at first

 

x ˙ = α x with α < 0.

(7.37)

For a given value x n we obtain an approximation value where x n + 1

 

x n + 1 = e α h x n

(7.38)

By applying the Runge-Kutta method of fourth order in equation (7.38), we obtain the following statement

 

x n + 1 = [ 1 α h + 1 2 ( α h ) 2 1 6 ( α h ) 3 + 1 24 ( α h ) 4 ] γ x n

(7.39)

It is obvious that the factor γ just contains the first five summands of the expansion power series where e α h . A comparison with the real value of the exponential function reveals that the deviation from the true value increases as α < 0 and h grows and that instability is observable when α h < 2.785 , [Figure 7.8]. This is due to the numerical solution which increases with each integration section ( | γ | < 1 ) while the true solution decreases ( e α h < 1 ) .

We receive similar results when we make the same analysis with further Runge-Kutta methods.

Similar to Euler method, we can more generally consider the case that α is a complex number. Here, we also analyze of the behavior of oscillating solutions. In this case, we do not obtain an interval of the real axis as a stable or instable area but a region in the complex plane, [Figure 7.10]

Instability area in the Runga-Kutta method of 4th order.
Figure 7.8. Instability area in the Runga-Kutta method of 4th order.


7.1.6.2. General Runga-Kutta Methods

The calculation formula of an ODE discussed in section [Section 7.1.6.1]

 

x ˙ = f ( x , t )

(7.40)

can be generally represented for a method with m function evaluations ( m -stepped method)

 

x n + 1 = x n + h j = 1 m β j k j

(7.41)

 

k j = f ( t n + ς j h , x n + h l = 1 m γ j l k l ) , j = 1,.., m

(7.42)

The values x n j = x n + h l = 1 m γ j l k l can be interpreted as an approximation of the solution with one integration section in the following points in time: t n + ζ j h

Transition points in the Runge-Kutta method.
Figure 7.9. Transition points in the Runge-Kutta method.


We choose the coefficients to approximate x n as much as possible. A clear representation of the coefficients β i , ζ i γ i j 1 i , j m is frequently based on the following schema (Butcher schema)

 

ς 1 γ 11 γ 12 γ 1 m ς 2 γ 21 γ 22 γ 2 m ς m γ m 1 γ m 2 γ m m β 1 β 2 β m

(7.43)

It is real for

 

ς i = j = 1 m γ i j , i = 1,..., m .

(7.44)

If

 

γ i j = 0 for j i .

(7.45)

the variables x i can be directly calculated from already known quantities, i.e. we refer to an explicit method (e.g. the classical Runge-Kutta methods already discussed in section [Section 7.1.6.1]).

Example 7.4: Butcher Schema for Classical Runge-Kutta Methods

Euler’s method:

0 0 1

Runge-Kutta method of 4th order:

0 0 0 0 0 1 2 1 2 0 0 0 1 2 0 1 2 0 0 1 0 0 1 0 1 6 1 3 1 3 1 6

Annotations 7.2:

  1. From section [Section 7.1.6.1] we infer that that there is at least one Runge-Kutta method of p = m order where m 4 . The question emerges whether there also exists a method where p < m and whether there is always, at least, one method where p = m . One must negate this in both cases.

  2. As a matter of fact, we can prove the following connection ([Table 7.1]):

Table 7.1.  Accumulation of the amount of transition values and achieved order.

Number of interim values

1

2

3

4

5

6

7

8

9

10

Achievable order

1

2

3

4

4

5

6

6

7

7


i.e. the method of 4 t h order represents a kind of optimum in respect to this observation.

7.1.6.3. Stability Area of Runge-Kutta Methods of Order 1≤p≤4

In this case we can prove that the area of absolute stability of a method of order p results from

 

| 1 + h λ + ( h λ ) 2 2 ! + ... + ( h λ ) p p ! R ( z ) = R ( h λ ) | 1

(7.46)

This particularly means that all methods of order p have the same stability area.

The stability border results from the complex solution z of the statement

R ( z ) = e i θ

with arbitrary θ of the interval ( 0,2 π ) .

The stability areas of the Runge-Kutta method up to p = 4 are represented in [Figure 7.10]

Stability area of the Runge-Kutta method of order 1-4.
Figure 7.10. Stability area of the Runge-Kutta method of order 1-4.


7.1.7.  Step Size Control

Errors are inevitable if iterative methods are applied. The emerging error depends on the used step size. Practically however, mostly a desired accuracy is merely known. Thus, one would like to provide the accuracy instead of the step size as input of the numerical integration method.

For this purpose, the algorithm is supposed to choose a step size automatically in every iteration step, such that the local discretization error lies below a desired value ζ . Thus, the error which is made in an iteration step must be estimated in advance.

Therefore, two different possibilities exist. On the one hand, an iteration step can be performed twice by using different step sizes. From the difference of the calculated solutions the local error can be estimated. On the other hand, two different solutions can also be obtained by using two different integration algorithms (embedded methods).

7.1.7.1. 3. Order Runge-Kutta, Two Calculations of an Integration Step

For the local discretization error of the Runge-Kutta method of third order with an integration interval of length h it yields

 

E h = B h 4 .

(7.47)

where the constant B depends from the envisioned differential equation. If one integrates the same interval in two integration steps of length h 2 , the following approximation statement of the discretization error at the end of the interval is true

 

2 E h 2 = 2 B ( h 2 ) 4 = 1 8 B h 4

(7.48)

If one subtracts the equation (7.47) from (7.48), we obtain

 

E h 2 E h 2 = B h 4 1 8 B h 4 = 7 8 B h 4

(7.49)

The left-hand side of the equation (7.49) can be calculated by first accomplishing an integration step of length h and subsequently repeating the calculation with the integration increment h 2 . If the specific results are referred to as x h and x h 2 respectively, we obtain in this case

 

E h 2 E h 2 = x h x h 2

(7.50)

By replacing equation (7.49) with (7.50) and after solving for B , it yields

 

B = 8 7 ( x h x h 2 ) h 4

(7.51)

If B is known, we can calculate an estimated value by means of equation (7.47)

 

h ( ζ B ) 1 4 .

(7.52)

7.1.7.2. Embedded Methods

It is not efficient to calculate the same integration step twice in succession in order to estimate the local integration error. For this reason, embedded methods are more advantageous. Hereby, two different algorithms are applied which use the greatest possible number of steps. The Runge-Kutta-Fehlberg 4 / 5 algorithm (RKF4/5) is a possible choice for an embedded algorithm. The Butcher scheme has the following form:

 

0 0 0 0 0 0 0 1 2 1 4 0 0 0 0 0 3 8 3 32 9 32 0 0 0 0 12 13 1932 2197 7200 2197 7296 2197 0 0 0 1 439 216 8 36801 513 845 4104 0 0 1 2 8 27 2 3544 2565 1859 4104 11 40 0 x 1 25 216 0 1408 2565 2197 4104 1 5 0 x 2 16 135 0 6656 12825 28561 56430 9 50 2 55

(7.53)

Thereby, x 1 is a Runge-Kutta solution of fourth order and x 2 is a Runge-Kutta solution of order five.

Therefore,

 

ε ~ h 5 h ~ ε 5 .

(7.54)

The relative error can be obtained by

 

ε r e l = | x 1 x 2 | max ( | x 1 | , | x 2 | , δ )

(7.55)

with δ = 10 10 .

The objective is to choose a new step size, such that the relative error lies near the relative tolerance.

We thus want the following equation to be fulfilled:

 

ς = | x 1 x 2 | max ( | x 1 | , | x 2 | , δ ) .

(7.56)

Therefore, the proposal for the choice of the step size is

 

h n e w = ς · max ( | x 1 | , | x 2 | , δ ) | x 1 x 2 | 5 · h o l d .

(7.57)

If the error is too great, the step size is decreased, if the error is too small, the step is increased. Notice, that hereby, steps are never repeated, even if the error is much too great. Therefore, this type of step size control is called optimistic. In contrast, a conservative step size control repeats a step with a new step size if the estimated error is greater than the tolerance.

The idea of step size control can also be understood as a (feedback) control problem. In this case, the above given formula equals a P-controller ([Figure 7.11]).

Step size control depicted as a (feedback) control problem
Figure 7.11. Step size control depicted as a (feedback) control problem


Kjell Gustafsson developed a PI-controller in order to control the step size. The corresponding formula is:

 

h n e w = ( 0.8 · t o l r e l ε r e l n e w ) 0.3 n · ( ε r e l o l d ε r e l n e w ) 0.4 n h o l d .

(7.58)

7.1.8.  Linear Multi-Step Methods

The previous integration methods had all in common that the approximation value x n + 1 directly results from its predecessor value x n . Within one integration step, we exclusively use the information of the last step. Methods with these characteristics are referred to as one-step procedures.

The question emerges whether one can achieve an enhanced accuracy by consulting the calculated values of the previous steps x n k , x n k + 1 , , x n 1 . These methods are called multi-step procedures.

The derivation of such procedures includes the formal integration of the initial value problem (7.1), (7.2) (compare the derivation of the Runge-Kutta method in section [Section 7.1.6])

 

x ( t n + 1 ) = x ( t n k + 1 ) + t n k + 1 t n + 1 f ( x ( τ ) , τ ) d τ .

(7.59)

The integrand will subsequently be integrated approximately by an interpolative quadrature formula with grid points t n k + 1 ,..., t n , t n + 1 . If k = 2 and the grids equidistant and the increment h , we then obtain for e.g The Simpson’s Rule:

 

x n + 1 x n 1 = h [ 1 3 f ( x n + 1 , t n + 1 ) + 4 3 f ( x n , t n ) + 1 3 f ( x n 1 , t n 1 ) ] .

(7.60)

If , on the other hand, we act on the original differential equation

 

x ˙ ( t ) = f ( x ( t ) , t ) | t = t n + 1

(7.61)

by approximating the derivation directly on the basis of a numerical differential formula, as e.g.

 

x ˙ ( t n + 1 ) 1 2 h [ 3 x ( t n + 1 ) 4 x ( t n ) + x ( t n 1 ) ]

(7.62)

we obtain the following procedure

 

3 2 x n + 1 2 x n + 1 2 x n 1 = h f ( x n + 1 , t n + 1 )

(7.63)

Both methods are examples of linear multi-step procedures. This will lead us to the following definition.

Definition 7.4

A linear multi-step procedure with n increments (also: linear n -step procedure), which aims at the determination of the approximations x n for the solution x ( t ) , with the initial value problem (7.1), (7.2) is defined by the specification of n initial values

 

x ( t j ) = x j , j = 0,1,..., n 1

(7.64)

and the calculation rule (difference equation)

 

j = 0 n α j x n j + 1 = h j = 0 n β j f ( x n j + 1 , t n j + 1 )

(7.65)

with

 

α j , β j , α 0 0 and | α n | + | β n | < 0.

(7.66)

Annotations 7.3:

  1. The linear multi-step procedure given in (7.64) and (7.65) is referred to as linear because the increment function (method function) depends linearly from the function values f ( x n j + 1 , t n j + 1 )

 

Φ ( x n k + 1 ,..., x n + 1 , t i ; h ) j = 0 n β j f ( x n j + 1 , t n j + 1 )

(7.67)

  1. The condition α 0 0 guarantees that the implicit differential equation (7.65) holds an exact solution, at least for sufficient small increments h .

  2. By means of the condition | α n | + | β n | < 0 , the step number n is exactly determined.

  3. β 0 = 0 is true for explicit linear multi-step procedures.

Example 7.5:

  1. By inserting α 0 = 1 and α 1 = 1 in (7.65), we obtain an implicit procedure ( β 0 0 )

 

x n + 1 = x n + h j = 0 n β j f ( x n j + 1 , t n j + 1 ) ,

(7.68)

which is referred to as Adams-Moulton formula.

  1. If we proceed like in 1. , but insert β = 0 , we obtain the following procedure

 

x n + 1 = x n + h j = 1 n β j f ( x n j + 1 , t n j + 1 )

(7.69)

This procedure constitutes the class of Adams-Bashford formulae.

  1. If we approximate the derivative on the basis of backward differences in the differential equation (7.61), we obtain a category of implicit integration procedures

 

j = 0 n α j x n j + 1 = h β 0 f ( x n + 1 , t n + 1 )

(7.70)

with representatives which are referred to as Backward Difference Formulae or, abbreviated, BDF formulae, cp. Section [Section 7.1.11]. This class of procedures plays an important role in the solution of stiff initial value problems and in the solution of DAEs.

An example of a favored multi-step procedure is the method of Adams-Bashford, which is based on the following formula:

 

x n + 1 = x n + h 24 ( 55 f ( x n , t n ) 59 f ( x n 1 , t n 1 ) + 37 f ( x n 2 , t n 2 ) 9 f ( x n 3 , t n 3 ) )

(7.71)

The local method error yields

 

251 720 h 5 x ( 5 ) ( η n ) = O ( h 5 ) with η n [ t n , t n + 1 ] .

(7.72)

In addition to explicit linear multi-step methods, one often utilizes implicit procedures in practice. Reasons are:

  1. they are more accurate when compared to explicit procedures,

  2. they have considerably better stability characteristics and

  3. they have easy strategies to estimate errors and to control step sizes.

In order to calculate a good starting value, which is necessary for the solution of the non-linear conditional equation where x n + 1 , we use a multi-step procedure with the form of a predictor-corrector method. One example is the procedure by Adams-Moulton which is defined by the predictor (7.38) and the corrector

 

x n + 1 = x n 1 + h 720 ( 251 f ( x n + 1 , t n + 1 ) + 646 f ( x n , t n ) 264 f ( x n 1 , t n 1 ) + 106 f ( x n 2 , t n 2 ) 19 f ( x n 3 , t n 3 )

(7.73)

Both formulae can also be used independently from each other.

In this case, the first formula is an explicit method (right hand side does not depend on x n + 1 ), whereas the second formula is an implicit method (right hand side depends on x n + 1 ). The first formula predominantly aims at achieving a good starting value for the iterative solution of the second. The second formula describes an implicit method which must be solved iteratively. In practice, we often use one or two iteration steps. The method error for the Adam-Moulton procedure is as follows

 

3 160 h 6 x ( 6 ) ( η n ) = O ( h 6 ) , η n [ t n , t n + 1 ]

(7.74)

In most cases, a predictor-corrector procedure consists of an explicit method (the predictor) and an implicit method (the corrector) with an error order which is at least equal to the predictor.

An error estimation (e.g. required for the increment control) results from additional computing with a double increment, as described within the scope of the Runge-Kutta methods:

 

e n , h 1 31 ( f h ( x n , t n ) f 2 h ( x n , t n ) )

(7.75)

In (7.71) and (7.74) we had encountered a problem which is always involved in the application of multi-step procedures: In order to start the calculation we already need solution values at the four sampling points x 1 , x 2 , x 3 , x 4 . These points must be calculated by means of another method (e.g. Runge-Kutta procedure). A further disadvantage is the complicated way of increment variations in the integration process (in contrast to one-step procedures).

An advantage is that the calculation of a new solution value merely requires one analysis of the differential equation. In contrast to that, Runge-Kutta methods of similar error order basically require analysis of differential equation at several points.

Stability analysis can also be made for multi-step procedures. But to deal with it here would go beyond the scope of this lecture.

7.1.9.  Activation of Linear Multi-Step Procedures

Because only the starting value x 0 = x ( 0 ) is given in a initial value problem, we must determine further n 1 initial values with a sufficient accuracy in order to start the linear multi-step procedure. As a general rule, there are to basic strategies:

  1. We use an one-step procedure with a automatic step size control and calculate approximates with given accuracy for the starting values at points t 1 , , t n 1 . Here, the following problems might occur:

    1. The grid points which are automatically determined by the step size control are not necessarily equidistant. In this case, we must calculate the required values at equidistant sampling points by interpolating the calculated sampling points.

    2. Even if the increment h is constant, it does not have to be identical to the increments which are needed by the linear multi-step procedure to comply with the accuracy requirements. Therefore, we advantageously use a one-step procedure which has the same consistency order like the multi-step procedure.

  2. We use a family of multi-step procedures with ascending consistency order. We start with a method of lower order and successively increase the order. As discussed in 1 , the grid points might also not be equidistant.

7.1.10.  System of Differential Equations

The numerical integration of a system of differential equation ensues analogically to the procedure in scalar differential equation. This will be exemplified by the Runge-Kutta procedure. In this case, the application of (7.59) to the equation system yields

 

x ˙ = f ( t , x )

(7.76)

the calculation rule

 

k 1 = h f ( x n , t n ) k 2 = h f ( x n + k 1 , t n + 1 ) x n + 1 = x n + 1 2 [ k 1 + k 2 ]

(7.77)

i.e. the scalar quantities must be substituted by vectors or, expressed differently, the numerical integration procedure will simply be applied equation wise.

7.1.11.  BDF Methods

[10]

The multi-step procedures discussed before are based on numerical solutions of the integral equations (7.59).The class of so-called “Backward Difference Formulae“ (BDF methods) will be constructed with the help of numerical differentiation.

In order to determine an approximate value of x n + 1 for x ( t n + 1 ) , we define an interpolation polynomial q by the points

 

( t n k + 1 , x n k + 1 ) ( t n k + 2 , x n k + 2 ) ( t n , x n ) ( t n + 1 , x n + 1 )

(7.78)

The polynomial q ( ζ ) can be written as follows

 

q ( ς ) = x ( t n + ς h ) = j = 0 k ( 1 ) j ( ς + 1 j ) j x n + 1

(7.79)

with the so-called backward differences

 

0 x n : = x n , j + 1 x n : = j x n j x n 1 .

(7.80)

The unknown value xn+1will be determined in such a way that the polynomial satisfies the differential equation

 

q ˙ ( t n + l ) = f ( t n + l , x n + l ) , l { 0, 1, 2, ... }

(7.81)

at one point of equation (7.78).

When l = 0 we obtain explicit formulae, namely when k = 1 the explicit Euler method and when k = 2 the Midpoint Method (cp. exercise). The formulae when k = 3 are instable.

When l = 1 , we obtain implicit formulae, the so-called BDF methods.

 

j = 0 k α j j x n + 1 = h f n + 1

(7.82)

with the coefficient

 

α j = ( 1 ) j d d ς ( ς + 1 j ) | ς = 1

(7.83)

With

 

( 1 ) j ( ς + 1 j ) = 1 j ! ( ς 1 ) ς ( ς + 1 ) ... ( ς + j 2 )

(7.84)

we obtain

 

α 0 = 0, α j = 1 j ! for  j 1

(7.85)

thus,

 

j = 1 k 1 j j x n + 1 = h f n + 1

(7.86)

Where k = 4 , we obtain for e.g.

 

25 x n + 1 48 x n + 36 x n 1 16 x n 2 + 3 x n 3 = 12 h f n + 1

(7.87)

The formulae (7.86) are stable when k 6 ; when k < 6 , there are instable.

7.1.12.  Remarks on Stiff Differential Equations

The following example will help us to understand the phenomenon which is referred to as stiff differential equations in literature. Let us first of all envision the following example:

Example 7.6:

We have the following equations

 

x ˙ 1 = λ 1 + λ 2 2 x 1 + λ 1 λ 2 2 x 2 , x ˙ 2 = λ 1 λ 2 2 x 1 + λ 1 + λ 2 2 x 2

(7.88)

with the constants λ 1 , λ 2 where λ i < 0 .

The exact general solution of this equation is

 

x 1 = c 1 e λ 1 t + c 2 e λ 2 t , x 2 = c 1 e λ 1 t c 2 e λ 2 t .

(7.89)

Both solutions converge to zero as t . Using the explicit Euler method, we obtain the numerical solution

 

x 1 n = c 1 ( 1 + h λ 1 ) i + c 2 ( 1 + h λ 2 ) i x 2 n = c 1 ( 1 + h λ 1 ) i c 2 ( 1 + h λ 2 ) i

(7.90)

These approximate solutions evidently only converge to zero, if

 

| 1 + h λ 1 | < 1 and | 1 + h λ 2 | < 1.

(7.91)

Here, h must be decreased in such way that

 

h < min ( 2 | λ 1 | , 2 | λ 2 | )

(7.92)

In this special technical process, which is described by this equation, suppose that | λ 2 | | λ 1 | .

In the analytical approach, the contribution of the term e λ 2 t is negligibly small when compared to the contribution of e λ 1 t . But in the numerical approach, the solution component of e λ 2 t becomes relevant for the choice of the minimal values of h due to eq (7.92).

Suppose that λ 1 = 1 and λ 2 = 1000 , it would be h 0.005 . If the second solution component does not exist, h 2 would be sufficient. Hence, although e 1000 t is more or less relevant to solve the problem, the factor 1000 determines the choice of the integration increment. This characteristic of the differential equation system in numerical integration is referred to as stiff. If the system at hand is stiff, not only the A -stability of the integration method has to be considered but the damping behavior as well. For this purpose the term L-stability is introduced.

Definition 7.5

A numerical integration method for differential equations is called L-stable if and only if it is A -stable and the following additionally holds:

 

| x n | 0 and Re ( h λ ) .

(7.93)

Annotation 7.4:

  1. All F -stable methods are A -stable, but never L -stable.

  2. While the solution portions including high damping are only weakly numerically damped when using F -stable methods, these portions are strongly numerically damped when using L -stable integration methods.

A principle difference between the characteristic of implicit and explicit integration methods will be discussed in the following.

Example 7.7: Integration of the Equation of Motion for Single Mass Pendulum [9]

The motion of a linear not stimulated single mass pendulum will be described by the state equation

 

x ˙ 1 = x 2 x ˙ 2 = c m x 1 d m x 2

(7.94)

For simplification, we write

x : = x 1

v : = x 2

Hence, the equations (7.94) yields

 

x ˙ = v v ˙ = c m x d m v

(7.95)

Application of the explicit Euler method on (7.95) results in

 

x k + 1 = x k + h v k v k + 1 = v k h ( c m x k + d m v k )

(7.96)

In the damped case ( d = 0 ) , the pendulum oscillates with constant amplitudes and the frequency

f 0 = ω 0 2 π with ω 0 = c m

The amplitude results from the initial condition.

The parameters

m = 1, c = 50, d = 0, x ( 0 ) = 1, v ( 0 ) = 0

with the increments

h = 1 m s , 5 m s , 10 m s

result in the outcomes represented in [Figure 7.12].

Integration of one-mass pendulum by means of the explicit Euler method.
Figure 7.12. Integration of one-mass pendulum by means of the explicit Euler method.


It becomes obvious that the oscillating amplitude does not remain constant (in contrast to reality), but that it grows with increasing increment.

In order to describe this characteristic, we calculate the total energy of a pendulum:

E ( x , v ) = 1 2 c x 2 + 1 2 m v 2

By means of the solutions by the Euler method, we obtain

 

E k + 1 = 1 2 c x k + 1 2 + 1 2 m v k + 1 2 = 1 2 c ( x k + h v k ) 2 + 1 2 m ( v k h c m x k ) 2 = 1 2 c x k 2 + c x k h v k + 1 2 c h 2 v k 2 + 1 2 m v k 2 m h c m x k + 1 2 m c 2 m 2 x k 2 = 1 2 c x k 2 + 1 2 m v k 2 + c h 2 m ( 1 2 c x k 2 + 1 2 m v k 2 ) = ( 1 + c h 2 m ) E k

(7.97)

(7.97) shows obviously that the total energy of the numerical solution increases with the factor

1 + c h 2 m = 1 + ( h ω 0 ) 2

in every step.

With the time period

T = 1 f 0 = 2 π ω 0

and the number of integration steps per period,

n = T h = 2 π h ω 0

the total energy increases per period with the factor of

( 1 + ( h ω 0 ) 2 ) n = ( 1 + ( h ω 0 ) 2 ) 2 π h ω 0 < 1

Since E = 1 2 c x max 2 at v = 0 (reversal point), the amplitude increases for each oscillation with the factor

( 1 + ( h ω 0 ) 2 ) 2 π h ω 0 .

Suppose this consideration would be true for the implicit Euler method, we would first obtain

 

x k + 1 = x k + h v k + 1

(7.98)

 

v k + 1 = v k h ( c m x k + 1 + d m v k + 1 )

(7.99)

If we insert equation (7.98) in (7.99), we obtain

v k + 1 = v k h c m x k ( h 2 c m + h d m ) v k + 1

or after transformation

( 1 + h 2 c m + h d m ) v k + 1 = v k h c m x k

and, subsequently, after calculation

v k + 1 = ( 1 + h 2 c m + h d m ) ( 1 + h 2 c m + h d m ) v k ( h 2 c m + h d m ) ( 1 + h 2 c m + h d m ) v k h 2 c m + h d m ( 1 + h 2 c m + h d m ) x k = v k h ( c m + h d + h 2 c x k + c h + d m + h d + h 2 c v k ) = v k h ( c m h x k + d h m h v k )

with the modified damping

d h = d + h c

and the modified mass

m h = m + h d + h 2 c = m + h d h

Suppose that the same parameter values of the explicit Euler method would also be real in this case. Hence, we will obtain the following results represented in [Figure 7.13].

Integration of one-mass oscillator by means of the implicit Euler method.
Figure 7.13. Integration of one-mass oscillator by means of the implicit Euler method.


Obviously, the oscillating amplitude decreases. After energy observation, analogue to the explicit method, we now obtain a decrease in energy for each oscillating period with the factor

( 1 + ( h ω 0 ) 2 ) 2 π h ω 0 < 1

The behavior of the total energy is represented in [Figure 7.14] in both cases.

Energy behaviour of a plane oscillator.
Figure 7.14. Energy behaviour of a plane oscillator.


This example suggests an essential general characteristic of integration methods, which will be discussed later:

  1. Explicit methods insert additional (only numerical qualified) excitement into the system.

  2. Implicit methods lead to an additional (numerical qualified) damping of the system.

In fact, this behavior also reoccurs in complex applications in practice.

A-stable methods are best suited to calculate stable problems (cp. [Section 7.1.2]). But in the past, we have only been dealing with the implicit Euler methods and the Trapezoidal Rule, with a local consistency order of 1 and 2 respectively, as an integrated part of A -stable methods. According to the theorem by Dahlquist, there is however no “better” method.

Theorem 7.1 (Dahlquist):

  1. Explicit multi-step methods are never A -stable.

  2. The order of an A -stable implicit multi-step method is at the most 2 .

  3. The Trapezoidal Rule is an A -stable method of second order with the lowest error constants.

Therefore, we introduce a further definition, which even though weakens the A -stability approach, can be expediently be utilized in many applications.

Definition 7.6:

A method is A ( α ) -stable if a stability area comprises a sector of the following form

 

A ( z ) : = { z C : | arg ( z ) | α } .

(7.100)

It is referred to as A -stable, if A ( 0 ) -stable where α < 0 (see [Figure 7.15]).

A(0)-stability.
Figure 7.15. A(0)-stability.


The BDF-methods are e.g. A ( α ) -stable with the opening angle α ([Table 7.2]).

Table 7.2.  Opening angle α for different k

k

1

2

3

4

5

6

α [ ]

90

90

86

73,3

51,8

17,8


7.1.13.  Implicit Runge-Kutta Methods

In comparison to the explicit Runge-Kutta methods we dealt with in section [Section 7.1.6], the matrix Γ in equation (7.43) can be completely filled. This means that the values k i in equation (7.42) cannot be sequentially calculated, but that with each integration step we have to solve a non-linear equation system.

Example 7.8: Implicit Euler Method

The implicit Euler method

x n + 1 = x n + h f ( x n + 1 , t n + 1 )

is obviously a single-step implicit Runge-Kutta method.

The Trapezoidal Rule

x n + 1 = x n + h 2 ( f ( x n , t n ) + f ( x n + 1 , t n + 1 ) )

can be regarded as a two-step implicit Runge-Kutta method.

k 1 = f ( x n , t n )

k 2 = f ( x n + h ( 1 2 k 1 + 1 2 k 2 ) , t n + h )

x n + 1 = x n + h ( 1 2 k 1 + 1 2 k 2 )

Butcher’s scheme is according to (7.43) as follows

Implicit Euler method

1 1 1

Trapezoidal Rule

0 0 0 1 0,5 0,5 0,5 0,5

A great disadvantage of implicit Runge-Kutta methods is the necessity to solve a non-linear equation system for the m · N variables h i (with N representing the dimensions of the equation system which must be solved). On the other hand, the additional constants in (7.43) can thereunto be used in order to

  1. achieve a higher consistency order with identical number of steps m .

  2. seriously improve the stability of the numerical method, which is “almost” A -stable.

Example 7.9: The Midpoint Method

The “Midpoint Method”

k 1 = f ( x n + h 2 h k 1 , t n + h 2 )

x n + 1 = x n + h ( 1 2 k 1 + 1 2 k 2 )

is obviously a single-step method which, however has a consistency order of 2.

Within the framework of this lecture, we will not cater to the construction and the application of implicit Runge-Kutta methods in detail. Detailled data is given by Jumann, M. (2004) for example.

7.1.14.  Comparison of Methods for Numerical Solution of Initial Value Problems (IVP)

The described methods can be divided into the following three classes

  1. Euler method,

  2. Runge-Kutta method,

  3. Multi-step procedure,

  4. BDF method

Furthermore, there exist other classes (e.g. extrapolation technique) which will not be discussed in this lecture.

Runge-Kutta methods and multi-step procedures allow for error estimation and increment step control. Additionally, multi-step procedures are open to variations of the method order. Commercially distributed program libraries (e.g. IMSL, NAG) offer sophisticated computer programs for all classes. Furthermore, there are programs with published source code which are freely available (e.g. Shampine and Gordon, 1983).

Advantages and Disadvantages of this Application:

In respect to the analyses of the differential equation, multi-step procedures call for the least effort. An explicit multi-step procedure often requires only a single analysis per integration step of the so-called right hand side. Implicit methods, on the other hand, further require the appropriate evaluations needed for the (“few”) iterations. Compared to this is the effort required for the increment control of the multi-step procedure. Multi-step procedures include relatively high efforts in terms of calculation control (overhead). To sum up, most of the advantages result from a complex structure of the differential equation as a result of which each analysis implies comparatively high calculating efforts.

Runge-Kutta methods of lower order prove to be more advantageous in cases of simply built differential equations or a discontinuity on the right-hand side. On the other hand, they become more disadvantageous in the case of complex differential equations. Furthermore, Runge-Kutta methods allow for very simple step increment regulations.

In the last years, Euler methods went through a kind of renaissance in the field of real time applications. Nevertheless, the potentially negative stability characteristics of the explicit Euler method should be considered. Thus, Euler methods can be recommended in case of stability problems, which can also cause problems because of the necessary iterations in real time applications. Furthermore, Euler methods require for very short integration steps (increments) because of their low consistency order.