Page 210 - Numerical methods for chemical engineering

196 4 Initial value problems

,τ m h −1 ,..., τ m h −k ]. We then generate a corrector polynomial that agrees
where a k = x[τ m h
with the predictor at uniformly spaced times in the recent past,
(p)
(c)
π (− j) = π (− j) j = 0, 1, 2,..., m h (4.192)
and that also matches the time derivative at t k+1 = t k + t,
(c)
dπ dx
= ( t) (4.193)
dτ τ=1 dt t k+1
(c)
With x [k+1] = π (1) and M [k+1] = M(x [k+1] ), this yields the equation
(c)

M [k+1] ˙ π (1) = ( t) f x [k+1] (4.194)
that we solve for the state x [k+1] at the new time t k+1 = t k + t.
(c)
Rather than construct π (τ), we find it easier to form the polynomial
(c)
(p)
δπ(τ) ≡ π (τ) − π (τ) δπ(− j) = 0 j = 0, 1, 2,..., m h (4.195)
Using Newton interpolation,
δπ(τ) = δπ[1] + δπ[1, 0](τ − 1) + δπ[1, 0, −1](τ − 1)(τ)
+ δπ[1, 0, −1, −2](τ − 1)(τ)(τ + 1) + ···
+ δπ[1, 0, −1, −2,. . . , −m h ](τ − 1)(τ)(τ + 1) ... (τ + m h − 1) (4.196)

(c)
Starting with δπ[1] = δπ(1) = π (1) − π (p) (1), from δπ(0) = 0 we have
δπ[0] − δπ[1] 0 − δπ(1)
δπ[1, 0] = = = δπ(1) (4.197)
(0 − 1) (−1)
Next, we use δπ(0) = δπ(−1) = 0 to compute
δπ[0, −1] − δπ[1, 0] 0 − δπ(1) 1
δπ[1, 0, −1] = = = δπ(1) (4.198)
(−1 − 1) (−2) 2
Similarly, as δπ(0) = δπ(−1) = δπ(−2) = 0,
1
δπ[0, −1, −2] − δπ[1, 0, −1] 0 − δπ(1) δπ(1)
2
δπ[1, 0, −1, −2] = = =
(−2 − 1) (−3) 3!
(4.199)
Continuing this approach, we find the general result
δπ(1)
δπ[1, 0, −1, −2,..., − j] = (4.200)
j!
The difference polynomial is then
m h +1 j−2

(
1 0
δπ(τ) = δπ(1) 1 + (τ + k) (4.201)
j!
j=1 k=−1
Taking the derivative of this polynomial at τ = 1,
 

m h +1 j−2 j−2 

1 0 
δ ˙π(1) = δπ(1)   (4.202)
j!  (1 + p)

 j=1 k=−1 p=−1 

p =k