Docstoc

Introduction to Matlab

Document Sample
Introduction to Matlab Powered By Docstoc
					    Numerical Solution of Single ODEs
1. Euler methods
2. Runge-Kutta methods
3. Multistep methods
4. Matlab example
                        Euler Methods
   Initial value problem
       dy
           f ( x, y)    y( x0 )  y0
       dx
    » Often impossible to generate analytical solution y(x)
    » Instead compute approximate solution at equally spaced
      node points
        x1  x0  h     x2  x0  2h     x3  x0  3h 
        y ( x1 )  y1   y ( x2 )  y 2   y ( x3 )  y3   

    » h = step size
   Taylor series expansion
                            dy ( x) h 2 d 2 y( x)
      y( x  h)  y( x)  h                 2
                                                  
                             dx      2! dx
                            Forward Euler Method
           First-order Taylor series expansion
                                      dy( x)
             y ( x  h)  y ( x )  h         y( x)  hf ( x, y)
                                       dx
           Iterative calculation
 y1  y0  hf ( x0 , y0 )     y2  y1  hf ( x1 , y1 )
yn 1  yn  hf ( xn , yn )       n  0,1,2,

           Approximation error
                                   dy ( x) h 2 d 2 y( )
             y( x  h)  y( x)  h                            x   xh
                                    dx      2! dx 2
            » Local truncation error: O(h2)
            » Global truncation error: O(h)  first-order method
               Backward Euler Method
   Forward Euler
      dy( x)                 yn1  yn
              f ( x, y)                f ( xn , yn )     yn1  yn  hf ( xn , yn )
       dx                        h

    » Explicit formula for yn+1 (explicit method)
   Backward Euler
     dy( x)                yn1  yn
             f ( x, y)               f ( xn1 , yn1 )   yn1  yn  hf ( xn1 , yn1 )
      dx                       h

    » Implicit formula for yn+1 (implicit method)
    » Allows larger h values to be used with comparable
      errors  more stable
    » Generally preferred to forward Euler
                    Plug-Flow Reactor Example

qi, CAi                                                               qo, CAo
                                CA(z)                       z

           0                                                     L
                                                                         
                                                                      A k B
                z                                                      r  kC A
                    q dC A
                                kC A  0 C A (0)  C Ai
                    A dz
                 q C A ( z n 1 )  C A ( zn )
                                                kC A ( z n 1 )  0
                 A              z
                                  1
          C A ( zn 1 )                   C A ( z n ) C A ( z0 )  C Ai
                          1  kAz q
                                                            N
                                              1     
                                         1  kAz q  C Ai
               C A ( L)  C A ( z N )              
                                                    
       Plug-Flow Reactor Example cont.
   Analytical solution
                           kA 
                           q L
      C A ( L)  C Ai exp    
                              
   Numerical solution
                               N                        N
                       1                    1    
                  1  kAz q  C Ai   1  kAL Nq  C Ai
      C A ( L)                                  
                                                 
   Convergence formula
                       N
            1 
       lim            e a
       N  1  a N 
                    
   Convergence of numerical solution
                           N
                   1                       kA 
       lim  1  kA L Nq  C Ai  C Ai exp   q L 
                                                  
       N 
                                                 
              Improved Euler Method
   Predictor-corrector formulation
      Prediction              yn 1  yn  hf ( xn , yn )
                               *


      Correction   yn 1  yn  1 h[ f ( xn , yn )  f ( xn 1 , yn 1 )]
                                 2
                                                                  *



   Approximation error
    » Local truncation error: O(h3)
    » Global truncation error: O(h2)  second-order method
   Adaptive step size
    » Large h desirable if dy/dx changing slowly (speed)
    » Small h necessary if dy/dx changing rapidly (accuracy)
    » Adjust h to maintain the local error below a prespecified
      tolerance
    » Used in modern ODE solution software
                     Runge-Kutta Methods
   Basic Runge-Kutta (RK) method
              k1  hf ( xn , yn ) k 2  hf ( xn  1 h, yn  1 k1 )
                                                  2         2

         k3  hf ( xn  1 h, yn  1 k 2 ) k 4  hf ( xn  h, yn  k3 )
                        2         2

                     yn 1  yn  1 (k1  2k 2  2k3  k 4 )
                                  6


     » Global truncation error: O(h4)  fourth-order method
   Runge-Kutta-Fehlberg (RKF) method involves a
    combination of two Runge-Kutta formulas (see text)
   Comparison of single-step methods
    Method              Function Evaluations    Global Error    Local Error
    Euler                         1                  O(h)            O(h2)
    Improved Euler                3                 O(h2)            O(h3)
    RK                            4                 O(h4)            O(h5)
    RKF                           6                 O(h5)            O(h6)
                   Multistep Methods
   Definition
    » Single-step: calculation at current step only uses value at
      last step (yn  yn+1)
    » Multi-step: calculation at current step uses values at
      several previous steps
   General development
    » Integrate ODE
        dy( x)               xn1                             xn1
                f ( x, y)  dy( x)  y( xn1 )  y( xn )   f [ x, y( x)]dx
          dx                xn                               xn

    » Approximate f(x,y) with an interpolation polynomial
                           xn1
          yn 1  yn            p( x)dx
                        xn

    » Different polynomial produce different methods
                Adams-Bashforth Methods
   Utilize cubic polynomial that interpolates
            f n  f ( xn , y n )           f n 1  f ( xn 1 , yn 1 )
        f n  2  f ( xn  2 , y n  2 )   f n  3  f ( xn  3 , y n  3 )

   Iterative calculation

        yn1  yn  24 h(55 f n  59 f n1  37 f n2  9 f n3 )
                    1



    » Explicit method
    » Global truncation error: O(h4)  fourth-order method
   Initialization
        y4  y0  24 h(55 f 3  59 f 2  37 f1  9 f 0 )
                  1



    » Must compute y1, y2, y3 by another method of comparable
      accuracy (e.g. Runge-Kutta)
                Adams-Moulton Methods
   Utilize cubic polynomial that interpolates

        f n 1  f ( xn 1 , yn 1 )       f n  f ( xn , y n )
        f n 1  f ( xn 1 , yn 1 )   f n  2  f ( xn  2 , y n  2 )

   Iterative calculation
        yn1  yn  24 h(9 f n1  19 f n  5 f n1  f n2 )
                    1



    » Implicit method
    » Global truncation error: O(h4)  fourth-order method
   Predictor-corrector formulation
       yn 1  yn  24 h(55 f n  59 f n 1  37 f n  2  9 f n 3 )
        *            1


          yn 1  yn  24 h(9 f n*1  19 f n  5 f n 1  f n  2 )
                       1
                      Matlab Example
   Isothermal CSTR model
                      dC A q
         
      2A  B k
                           (C Af  C A )  kC A  f (C A )
                                               2

                       dt  V

   Model parameters: q = 2, V = 2, Caf = 2, k = 1
   Initial condition: CA(0) = 2
   Backward Euler formula
                         q                     2 
    C A,n 1  C A,n  h  (C Af  C A,n )  kC A,n   C A,n  hf (C A,n )
                         V                         
   Algorithm parameters: h = 0.01, N = 200
Matlab Implementation: iso_cstr_euler.m
h = 0.01;           for i=1:N
N = 200;              t(i+1) = t(i)+h;
Cao = 2;              f = q/V*(Caf-Ca(i))-
q = 2;                k*Ca(i)^2;

V = 2;                Ca(i+1)= Ca(i)+h*f;

Caf = 2;            end

k = 1;              plot(t,Ca)

t(1) = 0;           ylabel('Ca (g/L)')

Ca(1) = Cao;        xlabel('Time (min)')
                    axis([0,2,0.75,2.25])
                                   Euler Solution

>> iso_cstr_euler
            2.2


             2


            1.8


            1.6
 CA (g/L)




            1.4


            1.2


             1


            0.8
                  0   0.2   0.4   0.6   0.8       1      1.2   1.4   1.6   1.8   2
                                              Time (min)
             Solution with Matlab Function
function f = iso_cstr(x)      >> xss = fsolve(@iso_cstr,2)
Cao = 2;                      xss = 1.0000
q = 2;                        >> df = @(t,x) iso_cstr(x);
V = 2;                        >> [t,x] = ode23(df,[0,2],2);
Caf = 2;                      >> plot(t,x)
k = 1;                        >> ylabel('Ca (g/L)')
Ca = x(1);                    >> xlabel('Time (min)')
f(1) = q/V*(Caf-Ca)-k*Ca^2;   >> axis([0,2,0.75,2.25])
                           Matlab Function Solution

           2.2


            2                                                          Euler
                                                                       Matlab

           1.8


           1.6
Ca (g/L)




           1.4


           1.2


            1


           0.8
                 0   0.2    0.4   0.6   0.8       1      1.2   1.4   1.6   1.8   2
                                              Time (min)

				
DOCUMENT INFO
Shared By:
Categories:
Tags:
Stats:
views:6
posted:10/6/2011
language:English
pages:16