Fabien B. Analytical System Dynamics: Modeling and Simulation

Подождите немного. Документ загружается.

208 5 Numerical Solution of ODEs and DAEs

methods.

Local error

The local error at time t

(k+1)

is deﬁned as

(k+1)

= y(t

(k+1)

) − y

(k+1)

, (5.8)

where it is assumed that the numerical solution coincides with the exact

solution at t

(k)

, i.e., y(t

(k)

) = y

(k)

. Hence, η

(k+1)

is a measure of the error

made in the single step from t

(k)

to t

(k+1)

In the case of the explicit Euler method the local error can be determined

by writing Φ(y(t

(k+1)

)) in a Taylor series about the point y

(k+1)

. Doing so

gives

Φ(y(t

(k+1)

)) = Φ(y

(k+1)

) + DΦ(y

(k+1)

)[y(t

(k+1)

) − y

(k+1)

]

+ higher order terms

= δ

(k+1)

But, Φ(y

(k+1)

) = 0 because y

(k)

= y(t

(k)

). Also, DΦ(y

(k+1)

) = dΦ/dy

(k+1)

I. Hence, if we neglect the higher order terms we get

y(t

(k+1)

) − y

(k+1)

= η

(k+1)

= δ

(k+1)

. (b)

Therefore, in the case of the explicit Euler method, the local error is equal

to the local discretization error.

For the implicit Euler method the local error can be determined by writing

Ψ(y(t

(k+1)

)) in a Taylor series about the point y

(k+1)

. That is,

Ψ(y(t

(k+1)

)) = Ψ (y

(k+1)

) + DΨ(y

(k+1)

)[y(t

(k+1)

) − y

(k+1)

]

+ higher order terms

= δ

(k+1)

where DΨ (y

(k+1)

) = dΨ/dy

(k+1)

= I −h

f(y

(k+1)

, t

(k+1)

), and Ψ(y

(k+1)

) =

0 because y

(k)

= y(t

(k)

). By neglecting the higher order terms we can see that

the local error is

y(t

(k+1)

) − y

(k+1)

= η

(k+1)



I − h

f(y

(k+1)

, t

(k+1)

)



−1

(k+1)

. (c)

If the term h df/dy is small then the local error for the explicit and implicit

Euler methods will be the same. For stiﬀ diﬀerential equations however, the

term h df/dy can be very large, and as a result the local error must be de-

termined by scaling the discretization error as shown in (c).

5.2 Higher-order Methods for ODEs 209

Global Error

Consider the ODE (5.1) with initial condition y(t

(0)

) = y

(0)

then the global

error at time t

(k)

is deﬁned as

(k)

= y(t

(k)

) − y

(k)

. (5.9)

Using the explicit Euler method in (5.9) we get

(k+1)

= y(t

(k+1)

) − (y

(k)

+ hf (y

(k)

, t

(k)

))

= y(t

(k+1)

) − (y(t

(k)

) − E

(k)

+ hf (y(t

(k)

) − E

(k)

, t

(k)

))

= E

(k)

+ y(t

(k+1)

) − y(t

(k)

) − hf(y(t

(k)

), t

(k)

)

+hf

(k)

+ O(kE

(k)

)

= (I + hf

(k)

+ δ

(k+1)

+ O(kE

(k)

)

≈ (I + hf

(k)

+ δ

(k+1)

where f

f(y(t

(k)

), t

(k)

). This result shows that the global error at t

(k+1)

is a function of the global error at t

(k)

and the local discretization error δ

(k+1)

Thus, the global error is a propagation of the previous errors as well as the

current the local discretization error.

In practice controlling the global error when integrating initial value prob-

lems can be very expensive. This involves solving the problem once to estab-

lish an estimate of the global error. If the desired global error is not achieved

the problem is solved a second time with reduced step sizes.

Since controlling the global error is too expensive most numerical methods

regulate the local error instead. This is done by requiring that the local error

be less than some speciﬁed solution tolerance at each step. The details of how

such schemes are implemented will be discussed below.

5.2 Higher-order Methods for ODEs

The previous section introduced the explicit and implicit Euler methods

which are both ﬁrst-order, i.e., the leading term in the local discretization

error is O(h

). These methods tend to be ineﬃcient when the interval [t

, t

]

is large and the desired solution tolerance is small. Since, in this case small

step sizes, h, will be required to maintain a small local error. It is therefore

desirable to have integration methods that allow us to take large steps while

maintaining small errors. This can be accomplished if the local discretization

error of the method is of a high order. In this section we consider some single

step higher-order numerical integration methods.

210 5 Numerical Solution of ODEs and DAEs

5.2.1 The Taylor series method

To develop this method let t be some time where the solution to the ODE

is known. Let h be some small ﬁxed increment in the time. Then, expanding

the exact solution y(t + h) in a Taylor series about the time t, gives

y(t + h) = y(t) + h

y(t) +

y(t) + ···+

y(t) + R(ξ),

(5.10)

where R(ξ) =

p+1

(p+1)!

p+1

y(ξ) is called the remainder, with ξ being some

point between t and t + h, and p > 0 being some integer. Here we assume

that y(t) and f (y(t), t)) are smooth enough to assure that the representation

(5.10) is valid.

Now let y

(k)

approximate the exact solution at time t, and let y

(k+1)

ap-

proximate the exact solution at time t + h. Then by neglecting the remainder

in (5.10), and using the fact that f(y(t), t) = dy(t)/dt, we can obtain the

explicit discretization formula

(k+1)

= y

(k)

+ hT

(k)

, t

(k)

), (5.11)

where

(k)

, t

(k)

) = f(y

(k)

, t

(k)

) +

f(y

(k)

, t

(k)

) + ···+

p−1

f(y

(k)

, t

(k)

(5.12)

The evaluation of T

(k)

, t

(k)

) requires the total derivatives of f (y

(k)

, t

(k)

For example in the case of a scalar ODE

= f +

where

∂f

∂y

f +

∂f

∂t

∂

∂y



∂f

∂y



f + 2

∂

∂t∂y

f +

∂f

∂y

∂f

∂t

∂

∂t

Here, it is understood that all functions are evaluated at the point (y

(k)

, t

(k)

From (5.10) it can be seen that the local discretization error for the formula

(5.12) is

(k+1)

p+1

(p + 1)!

p+1

y(ξ)

p+1

(p + 1)!

f(y(ξ), ξ),

5.2 Higher-order Methods for ODEs 211

where t

(k)

< ξ < t

(k+1)

. Thus the formula (5.11) deﬁnes a p-th order explicit

method.

The formula (5.12) can be used to construct a Taylor series method for the

numerical integration of ordinary diﬀerential equations. Such a procedure is

given in Algorithm 5.2.1. Note that the explicit Euler method (Algorithm

Algorithm 5.2.1 Taylor Series Method

Input: An integer p > 0, and integer N > 0, h = (t

− t

)/N , t

(0)

= t

, y

(0)

= y(t

) = y

Output: y

(k)

, k = 1, 2, · · · , N.

1: for k = 0, 1, · · · , N − 1 do

2: y

(k+1)

= y

(k)

+ hT

(k)

, t

(k)

)

3: t

(k+1)

= t

(k)

+ h

4: end for

5.1.1) is a special case of the Taylor series method. In particular, when p = 1

the formula (5.12) deﬁnes the explicit Euler method.

5.2.2 Explicit Runge-Kutta methods

The implementation of higher order Taylor series methods require the total

derivatives of f(y, t) which can be be tedious to compute. The explicit Runge-

Kutta methods attempt to duplicate the higher order Taylor series methods

without computing the higher order derivatives of f (y, t).

To illustrate how these Runge-Kutta methods can be constructed we con-

sider the development of a method that has a local discretization error of

order h

. For the sake of simplicity we restrict our attention to scalar ODEs.

An example of an explicit Runge-Kutta method is

(k+1)

= y

(k)

+ h [b

+ b

] , (a)

where

= f (y

(k)

, t

(k)

)

= f (y

(k)

+ ha

, t

(k)

+ c

h).

Here, the method coeﬃcients b

, b

, a

and c

are selected to ensure that

the local discretization error of (a) is of order h

The formula (a) deﬁnes a ‘2-stage’ Runge-Kutta method, where k

and k

are called the stage derivatives. Comparing (a) to a Taylor series expansion

of y(t

(k+1)

) we see that the terms b

must coincide with f +(h/2)

for this method to have a local discretization error of order h

. Hence, the

combination b

+ b

is used to approximate the total derivative of f , and

this is to be accomplished using only function evaluations.

212 5 Numerical Solution of ODEs and DAEs

The local discretization error for the formula (a) is computed as

(k+1)

= y(t

(k+1)

) − y(t

(k)

) − h



+ b



, (b)

where

= f (y(t

(k)

), t

(k)

)

= f (y(t

(k)

+ ha

), t

(k)

+ c

h).

Now, write y(t

(k+1)

and

as Taylor polynomials centered at the point

(y(t

(k)

), t

(k)

) to get

y(t

(k+1)

) = y(t

(k)

) + hf +

+ O(h

)

= y(t

(k)

) + hf +



∂f

∂y

f +

∂f

∂t



∂

∂y



∂f

∂y



f + 2

∂

∂t∂y

f +

∂f

∂y

∂f

∂t

∂

∂t

+ O(h

)

= f

= f + hc

∂f

∂t

∂

∂t

+ ha

∂f

∂y

f +

∂

∂y

+ a

∂

∂y∂t

f + O(h

Using these in (b) yields

(k+1)

= [(1 − b

− b

)f] h +



− b



∂f

∂y

f +



− b



∂f

∂t





−



∂

∂y



−



∂

∂t



− b



∂

∂y∂t

f +



∂f

∂y



f +

∂f

∂y

∂f

∂t

+ O(h

From this result it can be seen that the local discretization error, δ

(k+1)

, will

be of order h

1 − b

− b

= 0. (c)

If (c) holds along with

− b

= 0, (d)

and

5.2 Higher-order Methods for ODEs 213

− b

= 0 (e)

then, the local discretization error will be of order h

. It is also evident from

(k+1)

that there is no way for us to select b

, b

, a

and c

so that the

coeﬃcient of h

is made to vanish, for all possible functions f(y(t), t). Thus,

the highest possible order that can be achieved by this 2-stage Runge-Kutta

method is order 2.

Equations (c), (d) and (e) are called the order conditions for the Runge-

Kutta method. These equations provide us with three equations for the four

coeﬃcients b

, b

, a

and c

. Therefore, one of the coeﬃcients can be selected

arbitrarily and as a result an inﬁnite number of solutions are possible.

Two well known solutions to these order conditions are;

• The explicit midpoint method

= 0, b

= 1, a

, c

• Huen’s method/Explicit trapezoidal method

, b

, a

= 1, c

= 1.

It is easy to verify that both of these methods yield local discretization error

of order h

, and are thus called second-order methods.

The general s-stage explicit Runge-Kutta method has the form

(k+1)

= y

(k)

+ h

i=1

= f (Y

, τ

= y

(k)

+ h

i−1

j=1

= t

(k)

+ c

h. (5.13)

The coeﬃcients of the method are given in A ∈ R

s×s

, b ∈ R

and c ∈ R

Here, the (i, j)-th element of the matrix A is denoted as a

, etc. For explicit

Runge-Kutta methods the matrix A is strictly lower triangular, i.e., a

= 0

if j ≥ i.

The coeﬃcients of a Runge-Kutta method can be written compactly as a

Butcher tableau,

c A

For example, the explicit Euler method, the explicit midpoint method and

the explicit trapezoidal method are written as

214 5 Numerical Solution of ODEs and DAEs

0 0

0 0 0

1/2 1/2 0

0 1

, and

0 0 0

1 1 0

1/2 1/2

respectively.

Using the approach outlined above it is possible to derive the order con-

ditions for higher order Runge-Kutta methods. As shown above, the 2-stage

Runge-Kutta formula is at best a second-order method. To achieve higher-

order methods we must use more stages in the Runge-Kutta formula. How-

ever, as the number of stages increases so does the complexity of the anal-

ysis required to determine the order conditions. In Problem 5 the reader is

asked to derive the order conditions for a 3-stage, third-order explicit Runge-

Kutta method. The order conditions for a 4-stage, fourth-order Runge-Kutta

method are;

i=1

= 1,

i=1

j=1

i=1

j=1

i=1

j=1

i=1

j=1

k=1

. (5.14)

These 8 nonlinear equations must be solved to determine the coeﬃcients of

the fourth-order method.

Perhaps the most famous Runge-Kutta method is the ‘classical’ fourth-

order formula given by the Butcher tableau

5.2 Higher-order Methods for ODEs 215

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

This 4-stage method is often called ‘the Runge-Kutta method’ but, as we

have seen there are many other Runge-Kutta methods. In fact, there are

several popular explicit fourth-order Runge-Kutta methods.

The development of explicit Runge-Kutta methods of order p ≥ 5 is very

complex. One of the diﬃculties is that the number of order conditions that

must be solved becomes very large. For example, a method of order p = 5

has 17 order conditions, a method of order p = 6 has 37 order conditions,

and a method of order p = 7 has 85 order conditions. Moreover, to realize a

method of order p = 5 requires at least 6 stages, a method of order p = 6

requires at least 7 stages, and a method of order 7 requires at least 9 stages.

The interested reader should consult Butcher (1987), and Hairer, Norsett and

Wanner (1993) for additional details.

Error estimation and step size adaptation

An important consideration for the eﬃcient integration of diﬀerential equa-

tions is the ability of adjust the step size so that the local error satisﬁes a

desired tolerance at each step. To see how such schemes are possible consider

the numerical integration of the scalar ODE ˙y = f (y(t), t), with initial condi-

tion y(t

) = y

. Here we will use two Runge-Kutta methods, the ﬁrst method

has order p, and the second method has order p−1. From the analysis in

the previous sections we know that at the k-th step the p-th order method

produces a local error

y(t

(k+1)

) − y

(k+1)

= Ch

p+1

+ O(h

p+2

), (a)

where y

(k+1)

is the result of the Runge-Kutta formula, y(t

(k+1)

) is the exact

solution to the ODE, and C is a constant that depends on the Runge-Kutta

formula and the function f(y(t), t). On the other hand, the method of order

p−1 produces a local error

y(t

(k+1)

) − ˆy

(k+1)

+ O(h

p+1

), (b)

where ˆy

(k+1)

is the result of the (p−1)-th order Runge-Kutta formula, and

C is a constant that depends on the Runge-Kutta formula and the function

f(y(t), t). We assume that both methods have the same initial condition for

the k-th step, i.e., y

(k)

= ˆy

(k)

= y(t

(k)

From (a) we see that y(t

(k+1)

) = y

(k+1)

+ Ch

p+1

+ O(h

p+2

). Using this

result in (b) we obtain a computable estimate of the local error, i.e.,

216 5 Numerical Solution of ODEs and DAEs

(k+1)

= y(t

(k+1)

) − ˆy

(k+1)

= y

(k+1)

− ˆy

(k+1)

+ O(h

p+1

). (c)

At each step we would like |η

(k+1)

| to be less than or equal to some speciﬁed

tolerance . If |η

(k+1)

| ≤  then the integration from t

(k)

to t

(k+1)

is considered

to be successful, and we accept y

(k+1)

as the numerical solution to the ODE

at t

(k+1)

. However, if |η

(k+1)

| >  then we need to repeat the integration with

a smaller step size so that the local error tolerance is not violated.

The computation of a new step size is accomplished using the local error

estimate given by (c). Suppose |η

(k+1)

| >  then

|η

(k+1)

| = |ˆy

(k+1)

− y

(k+1)

| ≈ |

| > .

Now, let

h be the step size that satisﬁes the condition |¯η

(k+1)

| = |

| = .

Then, dividing |¯η

(k+1)

| by |η

(k+1)

| gives







|η

(k+1)

which shows that

h and h are related via

h = h





|η

(k+1)



. (d)

Thus, if |η

(k+1)

| >  then (d) computes a new step size

h < h which is

used to repeat the integration of the ODE at time t

(k)

. If on the other hand

|η

(k+1)

| ≤  then (d) computes a new step size

h ≥ h which can be used to

advance the solution from time t

(k+1)

, i.e., t

(k+2)

= t

(k+1)

In practice the new step size

h is restricted to ensure that it does not vary

to greatly from the old step size h. This can be done using a formula of the

type

h = h min(fac

, max(fac

, β(|η

(k+1)

|/)

−1/p

)), (e)

where 0 < fac

< fac

and 0 < β < 1. Typical values for these parameters

are β = 0.9, fac

= 0.2 and fac

= 5, which ensures that 0.2h ≤

h ≤ 5h.

If y(t) is an n-dimensional vector then the error estimation and step size

computation must take account of the possible variability in the scaling of

the elements of y(t). In this case the step size computation formula (e) is

replaced with

σ =







i=1

|η

(k+1)

atol

+ rtol

(k+1)







−1/p

h = h min(fac

, max(fac

, βσ)). (5.15)

Here, η

k+1

is the i-th element of the vector η

k+1

= ˆy

(k+1)

−y

(k+1)

, atol ∈ R

is the absolute error tolerance, and rtol ∈ R

is the relative error tolerance.

5.2 Higher-order Methods for ODEs 217

Note that the elements of atol and rtol are nonnegative. In this formulation,

for the local error to be acceptable it is suﬃcient that |η

k+1

| ≤ atol

rtol

(k+1)

|. If rtol

= 0 then we require that |η

k+1

| not exceed the absolute

error measure atol

. This is called absolute error control. If atol

= 0 then

the local error will be acceptable if (|η

k+1

|/|y

(k+1)

|) ≤ rtol

. This is called

relative error control.

Embedded Runge-Kutta Formulas

The local error estimation and step size control scheme described in the pre-

vious section requires Runge-Kutta methods of order p and p−1. A practical

and eﬃcient approach to implementing such a strategy is to use so called

‘embedded’ Runge-Kutta methods. In this approach y

(k+1)

is computed by

(5.13), and ˆy

(k+1)

is computed from the formula

ˆy

(k+1)

= y

(k)

+ h

i=1

, (5.16)

where the coeﬃcients

, i = 1, 2, ···, s are selected to ensure that the local

error y(t

(k+1)

) − ˆy

(k+1)

= O(h

). The eﬃciency of this method is derived

from the fact that both y

(k+1)

and ˆy

(k+1)

are computed using the same stage

derivatives k

, i = 1, 2, ···, s.

The Butcher tableau for these embedded methods are of the form

c A

An example for such a method can be obtained by combining the explicit

Euler methods with the explicit trapezoidal method to get

0 0 0

1 1 0

1/2 1/2

1 0

For this embedded method we compute

(k+1)

= y

(k)

+ h





ˆy

(k+1)

= y

(k)

+ hk

= f (y

(k)

, t

(k)

)

= f (y

(k)

+ hk

, t

(k)

+ h).