# Introduction to Matlab by liaoqinmei

VIEWS: 6 PAGES: 16

• pg 1
```									    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
» 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
   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)
   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)

```
To top