Page 206 - Numerical Methods for Chemical Engineering

Differential-algebraic equation (DAE) systems 195

integrator, the velocity Verlet rule,
( t) 2
r α ← r α + v α ( t) + F α α = 1, 2,..., N
2m α
F α
v α ← v α + ( t) α = 1, 2,..., N
2m α
(4.185)
compute new forces F α α = 1, 2,..., N
F α
v α ← v α + ( t) α = 1, 2,..., N
2m α

Differential-algebraic equation (DAE) systems

In the previous sections, we have treated systems described by a set of ODEs. We now
consider the addition of algebraic equations to obtain a DAE system. We show how the
BDF method can be modified to accommodate a system of mixed differential and algebraic
equations. Consider the DAE system
M(x)˙x = f (x) (4.186)

M(x) is a state-dependent mass matrix.If M(x) is nonsingular, we can decouple the set of
equations into standard ODE form
˙ x = M −1 f (x) (4.187)

We consider here the case where M(x) is singular, as when the system is modeled by a
combination of differential and algebraic equations. As a particular example, consider the
system
˙ y = F(y, z)
(4.188)
0 = G(y, z)
System (4.188) is expressed in the DAE form (4.186) as

I 0 ˙ y F(y, z)
M ˙ x = = = f (x) (4.189)
0 0 ˙ z G(y, z)

BDF method for DAE systems of index one

We now show how the BDF method can be modified to simulate a DAE system when M(x)
is singular (Ascher & Petzold, 1998). The first step is to generate an explicit predictor
polynomial that extrapolates the state behavior at past times (not necessarily uniformly
spaced) into the future,
(p)
π (τ j ) = x(τ j ) = x [k− j] τ j = t k− j − t k j = 0, 1, 2,..., m h (4.190)
t
This polynomial is constructed easily with Newton interpolation,
(p)
π (τ) = a 0 + a 1 (τ − τ m h ) + a 2 (τ − τ m h )(τ − τ m h −1 )
(4.191)
)(τ − τ m h −1 ) ··· (τ − τ 1 )(τ)
+··· + a m h +1 (τ − τ m h