Page 199 - Numerical Methods for Chemical Engineering
P. 199

188     4 Initial value problems



                   and substitute these expansions into (4.159),

                       N             N                       N              N
                      	   [k+1]  [ j]  	  [k]  [ j]         	   [k]  [ j]  	   [k+1]  [ j]
                         c    w   =    c w   + ( t) (1 − θ)A   c w    + θ A   c   w
                          j             j                       j              j
                       j=1          j=1                      j=1           j=1
                                                                                     (4.161)
                                   [ j]
                   Using Aw [ j]  = λ j w , collecting terms yields
                            N
                                  [k+1]  [k]  [k]              [k+1]        [ j]

                       0 =     − c    + c  + c ( t)(1 − θ)λ j + c  ( t)θλ j w        (4.162)
                                  j      j    j                j
                           j=1
                   Setting to zero each component of (4.162) yields
                                            [k+1]
                                           c       1 + ( t)λ j (1 − θ)
                                            j
                                                =                                    (4.163)
                                             [k]
                                            c        1 − ( t)λ j θ
                                             j
                   For simplicity, let us assume that each λ j of interest for the determination of absolute stability
                   is real. Then, λ j < 0 and absolute stability requires
                               [k+1]
                              c
                              j        1 − ω j (1 − θ)
                                   =                ≤ 1   ω j =−( t)λ j > 0          (4.164)
                               [k]
                              c         1 + θω j

                               j
                   Let us consider the three cases of interest:
                                        [k+1]

                                        c
                                        j           1 − ω j
                        explicit Euler     [k]       =        ≤ 1  only when ω j ≤ 2  (4.165)
                                         c           1
                                         j   θ=0
                                         [k+1]
                                        c
                                        j            1
                        implicit Euler     [k]       =        ≤ 1 ∀ω j > 0           (4.166)
                                         c           1 + ω j
                                         j  θ=1
                                         [k+1]          1
                                        c
                                        j             1 − ω j
                                                        2
                        Crank–Nicholson    [k]     =    1    ≤ 1 ∀ω j > 0            (4.167)
                                         c  j    θ=1/2    1 + ω j
                                                        2
                   For both the implicit Euler and the Crank–Nicholson method, we find that for any ω j > 0,
                   (4.164) is satisfied and the methods are absolutely stable.
                   Definition A method is A-stable if it is absolutely stable for all positive time steps.
                     Thus, the implicit Euler and Crank–Nicholson methods, or more generally (4.150) when-
                   ever θ ≥ 1/2, are A-stable. Here, we have only shown this to be true when all stable eigen-
                   values are real. For four values of θ, Figure 4.11 plots the modulus of the growth coefficient
                        [k+1]  [k]
                   µ j = c  /c  as a function of ω j in the complex plane. The region of absolute stability,
                         j    j

                   |µ j |≤ 1, lies within the contour µ j = 1. When the rules are biased more towards the


                   future state than the old state, θ ≥ 1/2, the method is A-stable.
                     By contrast, the explicit Euler method loses absolute stability whenever ω j > 2 for
                   any stable eigenvalue. To see what happens in this case, Figure 4.12 plots the numerical
                                                                              −t
                   trajectories for ˙ x =−x, x(0) = 1, for which the true solution is x(t) = e . The explicit
   194   195   196   197   198   199   200   201   202   203   204