Fabien B. Analytical System Dynamics: Modeling and Simulation

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

228 5 Numerical Solution of ODEs and DAEs

−h

i=1

− hˆγf(ˆy

(k+1)

, t

(k+1)

) = 0.

Taking y

(k+1)

as the estimate of the solution, the simpliﬁed Newton’s method

gives

ˆy

(k+1)

= y

(k+1)

+ h [I − hˆγJ]

−1

(

i=1

(

− b

f(y

(k)

, t

(k)

) + ˆγf(y

(k+1)

, t

(k+1)

)

= y

(k+1)

+ [I − hˆγJ]

−1

(

i=1

+ h

f(y

(k)

, t

(k)

)

+ hˆγf(y

(k+1)

, t

(k+1)

)

, (5.29)

where [e

, e

, ···, e

] = [(

−b

), (

−b

), ···, (

−b

)]A

−1

. Also, we use

the approximation ∂f(ˆy

(k+1)

, t

(k+1)

)/∂y ≈ J = ∂f(y

(k)

, t

(k)

)/∂y. Thus, the

local error estimate is

(k+1)

= y

(k+1)

− ˆy

(k+1)

= −[I − hˆγJ]

−1

(

i=1

+ h

f(y

(k)

, t

(k)

)

+ hˆγf(y

(k+1)

, t

(k+1)

)

. (5.30)

In the case of the Radau IIA, s = 3 method η

(k+1)

= O(h

s+1

) = O(h

Using (5.30) the step size can be adjusted as

σ =







i=1

|η

(k+1)

atol

+ rtol

(k+1)







−1/4

h = h min(fac

, max(fac

, βσ)). (5.31)

where atol ∈ R

is the absolute error tolerance, rtol ∈ R

is the relative

error tolerance, and the step size adjustment parameters 0 < β < 1, fac

fac

> 0. As discussed in the previous section, if the local error estimate

satisﬁes the desired tolerance, the new step size at t

(k+1)

h ≥ h. Otherwise,

the integration from t

(k)

is repeated with a new step size

h < h.

With these ingredients an implicit Runge-Kutta method for the solution

of the ODE (5.1)-(5.2) can be implemented as outlined in Algorithm 5.2.3.

5.3 Numerical Solution of DAEs 229

Starting at t = t

the algorithm executes the loop deﬁned by steps 2 through

20 until t = t

. In step 3 the algorithm solves the equations (5.22) via the

simpliﬁed Newton’s method. For these Newton iterations we use Z

(0)

= 0

however, other initial estimates are possible (see Section 5.3.5 below, and

Olsson and S¨oderlind, (1998)). If the simpliﬁed Newton’s iterations fail to

converge the step size is reduced by a factor 0.5 and the Newton’s iterations

are repeated. If the Newton’s iterations converge then the local error estimate

is computed along with a new step size

h. If the local error is suﬃciently

small, i.e., σ ≤ 1, then the solution advances to time t + h with y(t + h) = ¯y.

Otherwise, if σ > 1 the equations (5.22) are solved again at time t, with a

new step size

h < h.

Algorithm 5.2.3 Implicit Runge-Kutta Method

Input: An initial time t

, a ﬁnal time t

, the initial condition y(t

), an initial step size

, absolute error tolerance atol ∈ R

, relative error tolerance rtol ∈ R

, step size

adjustment parameters 0 < β < 1, fac

> fac

> 0, the minimum allowable step size

min

, and the maximum allowable number of iterations MAX ITER > 0.

Output: y(t

)

1: y = y(t

), t = t

, and h = h

2: for ITER = 0, 1, · · · , MAX ITER do

3: if the simpliﬁed Newton’s method, (5.24), converges then

4: ¯y = y +

i=1

5: η = − [I − hˆγJ]

−1



i=1

+ h

f(y, t) + hˆγf (¯y, t + h)



6: σ =



i=1



|η

atol

+rtol

|¯y





−1/4

h = h min(fac

, max(fac

, βσ))

8: if σ ≤ 1 then

9: if t = t

, then y(t

) = ¯y, STOP end if

10: t = t + h

11: y = ¯y

12: if t +

h > t

, then

h = t

− t end if

13: else

14: h =

15: end if

16: else

17: h = h/2

18: end if

19: if h < h

min

, then STOP end if

20: end for

21: STOP, too many iterations

230 5 Numerical Solution of ODEs and DAEs

5.3 Numerical Solution of DAEs

We now consider the extension of the methods described above to systems of

diﬀerential-algebraic equations (DAEs) of the form

Φ(y, ˙y, t) = 0, (5.32)

where t is the time, y(t) ∈ R

is the state vector, and ˙y = dy/dt is the state

derivative. It is assumed that Φ(y, ˙y, t) ∈ R

is continuously diﬀerentiable

with respect to all of its arguments. The system (5.32) is sometimes called

implicit diﬀerential equations (IDEs).

We note that the equation (5.32) includes ordinary diﬀerential equations

(5.1), since we can write

Φ(y, ˙y, t) = ˙y − f(y, t) = 0.

The description (5.32) also includes the Lagrangian DAEs (4.8). To show

this we simply take y = [q

, f

, s

, e

, λ

, µ

]

. Thus, (5.32) can include a

combination of diﬀerential equations and algebraic equations.

An important structural property of the system (5.32) is the diﬀerentiation

index. To deﬁne this quantity consider the system of equations

Φ(y, ˙y, t) = 0

Φ(y, ˙y, t) = 0. (5.33)

Then the diﬀerentiation index is the smallest integer p ≥ 0 for which the

system (5.33) can be used to determine an explicit system of ordinary diﬀer-

ential equations of the form ˙y =

f(y, t). In the following example we consider

the diﬀerentiation index for some typical dynamic systems described using

the Lagrangian DAEs.

Example 5.4.

(i) Consider the mass-spring-damper system shown here.

5.3 Numerical Solution of DAEs 231

Using y

as the generalized displacement, and y

as the generalized ﬂow

(velocity), the diﬀerential equations of motion for this system can be writ-

ten as the implicit diﬀerential equations

= ˙y

− y

= 0

= ˙y

− (k/m)y

− (b/m)y

− F/m = 0.

Therefore, Φ(y, ˙y, t) = [Φ

, Φ

]

, where y = [y

, y

]

. We can see that

˙y = [ ˙y

, ˙y

]

can be determined explicitly from the equations above. Hence,

no derivatives of Φ are required, and system has diﬀerentiation index p = 0.

In fact it is easy to see that a system of ordinary diﬀerential equations has

diﬀerentiation index 0.

(ii)Consider the dynamic system described by a damper and an applied force,

as shown below.

Using y

as the generalized displacement, and y

as the generalized ﬂow,

the diﬀerential equations of motion for this system can be written as the

implicit diﬀerential equations

= ˙y

− y

= 0,

= by

+ F = 0.

We note that ˙y

does not appear explicitly in this system of equations.

Taking one time derivative of the system however, gives the equation

b ˙y

F = 0,

which can be used to determine ˙y

. Hence, this system has diﬀerentiation

index 1.

(iii)Consider the dynamic system described by a spring and an applied force,

as shown below.

Using y

as the generalized displacement, and y

as the generalized ﬂow,

the diﬀerential equations of motion for this system can be written as the

implicit diﬀerential equations

232 5 Numerical Solution of ODEs and DAEs

= ˙y

− y

= 0,

= ky

+ F = 0.

We note that ˙y

does not appear explicitly in this system of equations.

Taking two time derivatives of the system however, gives the equation

k ˙y

F = 0,

which can be used to determine ˙y

. Hence, this system has diﬀerentiation

index 2.

(iv)Consider the simple pendulum shown below. In this model we will use y

and y

as the displacements, and y

= ˙y

, and y

= ˙y

as the corresponding

ﬂows.

The Lagrangian DAEs for this system can be written as

= ˙y

− y

= 0,

= ˙y

− y

= 0,

= m ˙y

+ 2y

= 0,

= m ˙y

+ 2y

− mg = 0,

= y

+ y

− l

= 0,

where y

is the Lagrange multiplier associated with the displacement con-

straint y

+ y

− l

= 0.

In these diﬀerential-algebraic equations ˙y

does not appear explicitly. How-

ever, three time derivatives of the equation Φ

= 0 produces

˙y

(2y

˙y

+ 2y

˙y

− ˙y

g).

Therefore, this system has diﬀerentiation index 3.

5.3.1 Hessenberg form for the DAE

In many cases it is possible to separate the variables in y so that we can write

(5.32) explicitly as coupled diﬀerential equations and algebraic equations. For

5.3 Numerical Solution of DAEs 233

example, it can be seen the Lagrangian DAEs have such a structure. Some

important structural arrangements of the DAEs are known as the Hessenberg

forms, which are deﬁned below. For the sake of simplicity we give the au-

tonomous form of these systems. The non-autonomous form can be deduced

by appending the diﬀerential equation ˙τ = 1, where τ represents the time

variable.

• Hessenberg Index-1 Form.

The DAEs are said to be in Hessenberg index-1 form if they can be written

˙x = f(x, z),

0 = g(x, z), (5.34)

where x(t) ∈ R

, z(t) ∈ R

, and the Jacobian, g

= ∂g/∂z, is nonsin-

gular along the solution to the DAEs. In this formulation x is called the

diﬀerential variable and z is called the algebraic variable.

These index-1 DAEs can be reduced to a system of ordinary diﬀerential

equations by taking one time derivative of the algebraic equation g(x, z) =

0. In fact,

g(x, z) = g

˙x + g

˙z = 0, where g

= ∂g/∂x. Since g

nonsingular we have ˙z = −g

−1

˙x = −g

−1

f. Hence, only one time

derivative of the algebraic equation is required to determine ˙z in terms of

x and z.

Clearly we can write (5.34) in the form of (5.32) by using y = [x

, z

]

Also note that the DAEs in Example 5.4-(ii) can be put in Hessenberg

index-1 form if we select x = y

and z = y

• Hessenberg Index-2 Form.

The DAEs are said to be in Hessenberg index-2 form if they can be written

˙x = f(x, z),

0 = g(x), (5.35)

where x(t) ∈ R

, z(t) ∈ R

, and the matrix g

is nonsingular along

the solution to the DAEs, where f

= ∂f/∂z.

Two time derivatives of the algebraic equation g(x) = 0 gives ˙z =

−(g

)

−1

f, where (g

= ∂(g

f)/∂x. Note that we can write

(5.35) in the form of (5.32) by using y = [x

, z

]

, and the DAEs in Ex-

ample 5.4-(iii) can be put in Hessenberg index-2 form if we select x = y

and z = y

• Hessenberg Index-3 Form.

The DAEs are said to be in Hessenberg index-3 form if they can be written

˙x = f(x, z),

˙z = h(x, z, u),

234 5 Numerical Solution of ODEs and DAEs

0 = g(x), (5.36)

where x(t) ∈ R

, z(t) ∈ R

, u(t) ∈ R

and the matrix g

nonsingular along the solution to the DAEs, where h

= ∂h/∂u. Taking

three time derivatives of the algebraic equation g(x) will allow us to reduce

the DAEs to a system of ordinary diﬀerential equations (see Problem 15).

We can write (5.36) in the form of (5.32) by using y = [x

, z

, u

]

, and

the DAEs in Example 5.4-(iv) can be put in Hessenberg index-3 form if

we select x = [y

, y

]

, z = [y

, y

]

and u = y

5.3.2 Implicit Runge-Kutta methods for DAEs

Implicit Runge-Kutta methods can be applied directly to the diﬀerential

equations (5.32), and the Hessenberg form DAEs (5.34), (5.35) and (5.36).

In the case of the diﬀerential equations (5.32) the k-th step of an s-stage

implicit Runge-Kutta method has the form

(k+1)

= y

(k)

+ h

i=1

0 = Φ(Y

, Y

, τ

), i = 1, 2, ···, s,

= y

(k)

+ h

j=1

= t

(k)

+ c

h, (5.37)

where Y

∈ R

is the stage derivative, and Y

∈ R

is the stage value.

Therefore, each step of the method requires that we solve the equations

Φ(Y

, Y

, τ

) = 0, i = 1, 2, ···, s, for Y

∈ R

and Y

∈ R

In the case of the Hessenberg index-1 DAE (5.34) the k-th step of an

s-stage implicit Runge-Kutta method can be written as

(k+1)

= x

(k)

+ h

i=1

(k+1)

= z

(k)

+ h

i=1

= f (X

, Z

), i = 1, 2, ···, s,

0 = g(X

, Z

), i = 1, 2, ···, s,

= x

(k)

+ h

j=1

5.3 Numerical Solution of DAEs 235

= z

(k)

+ h

j=1

, (5.38)

where X

∈ R

and X

∈ R

are the stage derivative and stage value for

the x variables, respectively. Similarly, Z

∈ R

and Z

∈ R

are the stage

derivative and stage value for the z variables, respectively.

We can develop similar formulas for the Hessenberg index-2 form and the

Hessenberg index-3 form. (See Problem 16.)

Hairer, Lubich and Roche (1989), and Hairer and Wanner (1996) have

investigated the properties of implicit Runge-Kutta methods applied to DAEs

in the Hessenberg form. An important set of results they develop show the

relationship between the local error and the structure of the DAEs for a given

set of Runge-Kutta coeﬃcients. Here, we summarize the results they obtain

for the local error behavior of the Radau IIA methods.

• Hessenberg Index-1 Form.

Suppose an s-stage Radau IIA Runge-Kutta method is applied to the

system (5.34) with consistent initial conditions (x

(0)

, z

(0)

), then the local

error in the x and z variables satisfy

(1)

− x(t

(1)

) = O(h

2s−1

(1)

− z(t

(1)

) = O(h

2s−1

where t

(1)

= t

(0)

+ h, and x(t

(1)

), z(t

(1)

) is the exact solution to the DAEs

at time t

(1)

. The initial conditions are consistent if g(x

(0)

, z

(0)

) = 0. It is

interesting to note that the local error is the same as that obtained when

the s-stage Radau IIA method is applied to ordinary diﬀerential equations.

• Hessenberg Index-2 Form.

Suppose an s-stage Radau IIA Runge-Kutta method is applied to the

system (5.35) with consistent initial conditions (x

(0)

, z

(0)

), then the local

error in the x and z variables satisfy

(1)

− x(t

(1)

) = O(h

2s−1

(1)

− z(t

(1)

) = O(h

where t

(1)

= t

(0)

+ h, and x(t

(1)

), z(t

(1)

) is the exact solution to the DAEs

at time t

(1)

. In this case the initial conditions are consistent if g(x

(0)

) =

0 and g

(0)

)f(x

(0)

, z

(0)

) = 0. Note here the reduction in the order of

the local error of the algebraic variable, z, when compared to Hessenberg

index-1 systems.

• Hessenberg Index-3 Form.

Suppose an s-stage Radau IIA Runge-Kutta method is applied to the

system (5.36) with consistent initial conditions (x

(0)

, z

(0)

, u

(0)

), then the

local error in the x, z and u variables satisfy

(1)

− x(t

(1)

) = O(h

2s−2

236 5 Numerical Solution of ODEs and DAEs

(1)

− z(t

(1)

) = O(h

(1)

− u(t

(1)

) = O(h

s−1

where t

(1)

= t

(0)

+ h, and x(t

(1)

), z(t

(1)

), u(t

(1)

) is the exact solution to

the DAEs at time t

(1)

. In this case the initial conditions are consistent

if g = 0, g

f = 0 and (g

f + g

h = 0, when evaluated at

the initial condition x

(0)

, z

(0)

, u

(0)

. Again we note that the order of the

local error, for the algebraic variables, is signiﬁcantly reduced here, when

compared to the order attained for index-1 systems.

5.3.3 Index reduction

From the results in the previous section we observe that, for the algebraic

variables (z and u), the order of the local error decreases as the diﬀerentia-

tion index increases. Therefore, we can assert that the computational expense

for solving diﬀerential-algebraic equations increases with increasing diﬀeren-

tiation index. That is, the higher the diﬀerentiation index, the greater the

computational eﬀort required to achieve a desired solution accuracy for the

algebraic variables since, smaller step sizes will be required.

Due to this order reduction in the local error it is desirable to reformulate

a system of diﬀerential-algebraic equations with high diﬀerentiation index

into a problem with lower diﬀerentiation index. One approach to reducing

the diﬀerentiation index of DAEs in Hessenberg form is to diﬀerentiate the

algebraic equations. Each such diﬀerentiation will reduce the index of the

original problem by one, as illustrated in the following example.

Example 5.5.

Consider the index-3 DAEs

˙x = f(x, z),

˙z = h(x, z, u),

0 = g(x).

If we replace the algebraic equation g(x) = 0 with

g(x) = 0 we obtain the

system

˙x = f(x, z),

˙z = h(x, z, u),

0 = g

(x)f(x, z),

which are index-2 DAEs.

5.3 Numerical Solution of DAEs 237

If we replace the equation g(x) = 0 with

g(x) = 0, in the original index-3

DAEs, we obtain the system

˙x = f(x, z),

˙z = h(x, z, u),

0 = (g

(x)f(x, z))

f(x, z)

+ g

(x)f

(x, z)f(x, z)

+ g

(x)f

(x, z)h(x, z, u),

which are index-1 DAEs.

A third time derivative of the algebraic equation g(x) = 0 will lead to system

of ﬁrst-order diﬀerential equations in x, z and u.

Unfortunately, the index reduction technique described in Example 5.5

does not always lead to satisfactory numerical results. The main reason for the

poor performance of these reduced index DAEs is that they do not explicitly

satisfy the constraint g(x) = 0. As a result, the numerical solution tends to

drift away from this constraint as the time variable increases (see Problem

9).

Nevertheless, there is an index reduction technique for index-3 DAEs that

works well in practice. This method can be described as follows.

Consider the index-3 LDAEs

˙q − f = 0,

M(q, f)

f + φ

(q)

λ + Υ (q, f) = 0,

φ(q) = 0. (5.39)

We assume that the displacement constraints are linearly independent, and

the inertia matrix, M, is nonsingular at all times. Taking two time derivatives

of the displacement constraint φ = 0 we can reduce this LDAEs to a system

of ordinary diﬀerential equations

˙x = F (x),

where x = [q

, f

]

and

F (x) =

−M

−1

(q)



−1



−1

(φ

f −φ

−1

− M

−1

To ensure that the solution to this ODEs does not drift from the constraints

φ = 0 and

φ = φ

f = 0 we solve the system

˙x = F (x),

0 = H(x), (5.40)