Lecture 34 - Ordinary Differential Equations - BVP

Document Sample
Lecture 34 - Ordinary Differential Equations - BVP Powered By Docstoc
					Lecture 34 - Ordinary Differential
        Equations - BVP
            CVEN 302
         November 28, 2001
  Systems of Ordinary Differential
         Equations - BVP
• Shooting Method for Nonlinear BVP
• Finite Difference Method
• Partial Differential Equations
          Shooting Method for Nonlinear
                   ODE-BVPs
• Nonlinear        y   f ( x , y , y  ) , a  x  b
                   y( a )  y , y( b )  y
  ODE
                               a                b


• Consider u  f ( x , u , u ) with guessed slope t
            u(a)  y , u( a )  t
                    a
• Use the difference between u(b) and yb to adjust u’(a)
• m(t) = u(b, t) - yb is a function of the guessed value t
• Use secant method or Newton method to find the correct t
  value with m(t) = 0
     Nonlinear Shooting Based on Secant
                  Method

• Nonlinear  y   f ( x , y , y  ) , a  x  b
             y( a )  y , h( y( b ), y ( b ))  0
   ODE                   a


   Step1 : use u(a)  y a , u (a)  t(1)  Error  m(1)
   Step 2 : use u(a)  y a , u (a)  t(2)  Error  m(2)
   Step 3 : use secantmethodto obtain a new estimate
                             t(i  1)  t(i  2)
          t(i)  t(i  1)                       m(i  1)
                            m(i  1)  m(i  2)
   Step 4 : iterate until t(i) - t(i - 1)  tol
MATLAB Example in Nonlinear Shooting
            Method
• Nonlinear shooting with secant method
    y   2 yy  , 0  x  1
   
    y( 0 )  1 , h( y( 1 ), y ( 1 ))  y( 1 )  y ( 1 )  0 .25  0
               exact solution y  1 /( x  1 )
• Convert to two first-order ODE-IVPs
                     let z 1  y , z 2  y 
                   
               z1  z 2 ,              z1 ( 0 )  1
               
                z   2 z 1 z 2 ,
                   2                    z2 ( 0 )  t
• Update t using the secant method
Nonlinear Shooting - Secant Method



                              y(x)



                           y’(x)
                Nonlinear Shooting
             Based on Newton’s Method

• Nonlinear  y   f ( x , y , y  ) , a  x  b
            
   ODE       y( a )  y a , y( b )  y b
u  f(x,u, u)                           u(a)  y a , u(a)  t k
                                          
v   vf u (x,u, u)  v f u (x,u, u)  v(a)  0, v (a)  1
Check for convergence of m(t)
                               if m  tol , stop
                               
 m( t k )  u(b, t k )  y b                                  m
                               Otherwise , t k  1  t k  v ( b , t )
                                                                    k
               Nonlinear Shooting
            Based on Newton’s Method
• Nonlinear u  f ( x , u , u ) , u( a )  y a , u( a )  t
            
  ODE-IVP m( t )  u( b , t )  y b
                                             0
• Chain u f ( x , u , u ) f x f u f u
                                                     
  Rule t          t              x t       u t u t
                               x and t are independent
                       u              u          u
           let v         , then v       , v  
                       t              t            t
          u f u f u
                                      v   f u v  f u v 
           t       u t u t
          u( a , t )  t  v' (a)  u/t  1
Nonlinear Shooting with Newton’s Method
• Solve ODE-IVP
                 ( y  )2                ( u ) 2
        y             ,      u            ,
                    y                      u
        y( 0 )  1 , y( 1 )  2  u( 0 )  1 , u( 0 )  1
                                 
• Construct the auxiliary equations
                      f ( u  )2
 f u ( x , u , u )                     ( u ) 2     2 u
                      u   u  2
                                    v            v       v
                                           u 2
                                                          u
 f ( x , u , u )  f   2 u    v ( 0 )  0 , v ( 0 )  1
 u                                
                      u     u
Nonlinear Shooting with Newton’s Method

• Calculate m(t) -- deviation from the exact BC
                    m( t )  u(b, t )  y b
                    
                     m t ( t )  u t ( b, t )
• Update t by Newton’s method
                    m(t i 1 )                    m
   t i  t i 1                  ; v  m t (t)      u t ( b, t )
                    m t (t i 1 )                 t
         Finite-Difference Methods

• Divide the interval of interest into subintervals
                               ba
        x0  a , x n  b , h       x i 1  x i
                                h
• Replace the derivatives by appropriate finite-difference
  approximations in Chapter 11
• Solve the system of algebraic equations by methods in
  Chapters 3 and 4
• For nonlinear ODEs, methods in Chapter 5 may be used
              Finite-Difference Method
• General Two-Point BVPs
      y   p( x ) y   q( x ) y  r ( x ), a  x  b
     
      y(a)   , y(b)  
• Replace the derivatives by appropriate finite-difference
  approximations
                y i  1  2 y i  y i 1                y i  1  y i 1
  y ( x i )                           , y ( x i ) 
                          h2                                   2h
                  h            h         h          h
                        xi-1       xi        xi+1
             Finite-Difference Method
• Central difference approximations
                                  y i 1  y i 1
                    y ( x i ) 
                                        2h
                                  y i 1  2 y i  y i 1
                    y ( x i ) 
                                            h2
• Tridiagonal system
        y i  1  2 y i  y i 1      y i 1  y i 1
                    2
                                  pi                  q i y i  ri
                  h                         2h
   
   
        h 
                        
     1  pi  y i 1  2  h q i y i   1  pi  y i  1  h 2 ri
                                 2
                                                  h 
       2                                         2 
                 Finite-Difference Method
• Central Difference ==> Tridiagonal system
   ( 2  h 2 q1 )         0                    0                  0            y1 
                                                                               
   1  (h/2)p 2  ( 2  h 2 q 2 ) 1  (h/2)p 2                     0            y2  
         0           1  (h/2)p 3  ( 2  h q 3 ) 
                                                   2
                                                                     0
                                                                                  
                                                                                 y3   
                                                                                      
                                                                            
                                                         ( 2  h 2 q n 1 )  y n 1 
         0                 0                    0                                     
               h 2 r1  ( 1  hp1 /2 ) y 0   h 2 r1  ( 1  hp1 /2 )  
                                                                            
                          h 2 r2                           h 2 r2           
                                                                            
                           2
                           h r3                              2
                                                              h r3             
                                                                          
                                                                            
              h rn 1  ( 1  hpn -1 /2 ) y n  h rn 1  ( 1  hpn -1 /2 ) 
              
                 2
                                                
                                                     2
                                                                               
          Finite-Difference Method for
                Nonlinear BVPs
• Nonlinear ODE-BVPs
            y   f ( x , y , y  ) , a  x  b
           
            y(a)   , y(b)  
• Evaluate fi by appropriate finite-difference
  approximations
                 y i 1  2 y i  y i 1
                            2
                                          fi  0
                          h
                 h           h        h          h
                      xi-1       xi       xi+1
        Finite-Difference Method for
              Nonlinear BVPs
• SOR method
                1
       yi             ( y i  1  2 y i  y i  1  h 2 f i )
            2( 1   )
• Iterative solution
• Convergence criterion

    h Q / 2
          2        0  Q  f y ( x , y , y' )  Q 
                    
              with 
  h  2 / P
             
                     f y ( x , y , y  )  P
                    
              Example 14.12 - MATLAB

• Note error in Text
                                 ( y )2
  y   f ( x , y , y  )             , y(0)  1 , y(1)  2
                                    y
              1
  yi                  ( y i  1  2 y i  y i  1  h 2 f i )
         2( 1   )

  f       
             ( y i 1  y i 1 ) / 2 h 2   
                                                ( y i 1  y i 1 )2
      i
                           yi                          4h2 yi

              fi : negative sign
       Chapter 15

Partial Differential Equations
       Classification of PDEs
• General form of linear second-order PDEs with
  two independent variables

au xx  bu xy  cu yy  du x  eu y  fu  g  0

• linear PDEs: a, b, c,….,g = f(x,y) only

b 2  4 ac  0 , Hyperbolic (2 real roots)
 2
b  4 ac  0 , Parabolic (1 double root)
b 2  4 ac  0 , Elliptic  (2 complex roots)

  Heat Equation: Parabolic PDE
• Heat transfer in a one-dimensional rod
      x=0                                     x=a
  g1(t)                                         g2(t)

       u       2u
           c       ; 0  x  a, 0  t  T
       t     x  2

          u(x,0)  f(x),        0 xa
          u(0,t)  g 1 ( t )
                                 , 0 t T
          u( a , t )  g 2 ( t )
   Discretize the solution domain in space
   and time with h = x and k = t
     t
         10

            9

            8

            7

            6

            5

            4

            3


  Time      2


(j index)   1

            0
                0   1   2   3   4   5   6   7   8   9   10


                                 space                  x
                                (i index)
Initial and Boundary Conditions
                            Explicit Euler method
          10

          9

          8

          7



u(0, t)
          6

          5
                                                              u(a, t)
= g1(t)   4                                                   = g2(t)
          3

          2

          1

          0
               0    1   2    3   4   5   6   7   8   9   10



                   Initial conditions : u(x,0) = f(x)
                 Heat Equation
 • Finite-difference
                                tj+1
                                       t
 u(x,t)                                              (i,j+1)
             t                   tj
                                             (i-1,j) (i,j)     (i+1,j)




                                x
                                           xi-1     xi xi+1              x
                           1
Forward-difference   ut  ( ui , j  1  ui , j )
                           k
Central-difference            c
  at time level j    cu xx  2 ( ui  1 , j  2 ui , j  ui  1 , j )
                             h
                         Explicit Method
  • Explicit Euler method for heat equation
              h   x  a / n , x i  ih
          let 
              k   t  T / m , t j  jk
                 1                            c
ut  cu xx       ( u i , j  1  u i , j )  2 ( ui  1 , j  2 u i , j  ui  1 , j )
                 k                           h

  • Rearrange
                              ck
   ui , j  1    ui , j  2 ( ui  1 , j  2 ui , j  ui  1 , j )   ck   c t
                             h                                      r 2 
                                                                      h     x2
                 ru i 1 , j  ( 1  2 r )ui , j  ru i  1 , j

                      Stability: 0  r  0.5
             Explicit Euler Method
           ui , j  1  ru i  1 , j  ( 1  2 r )ui , j  ru i  1 , j
• Stable
        r  0.01  ui , j  1  0.01 ui  1 , j  0.98 ui , j  0.01 ui  1 , j
        r  0.1  ui , j  1  0.1 ui  1 , j  0.8 ui , j  0.1 ui  1 , j
        r  0.4  ui , j  1  0.4 ui  1 , j  0.2 ui , j  0.4 ui  1 , j
        r  0.5  ui , j  1  0.5 ui  1 , j  0.5 ui  1 , j
• Unstable (negative coefficients)
       r  1  ui , j  1  ui  1 , j  ui , j  ui  1 , j
       
       r  10  ui , j  1  10 ui 1 , j  19 ui , j  10 ui  1 , j
       r  100  ui , j  1  100ui 1 , j  199 ui , j  100ui  1 , j
       
Heat Equation: Explicit Euler Method

      r = 0.5
Example: Explicit Euler Method
• Heat Equation (Parabolic PDE)
                 ut  cu xx ; 0  x  1
           u(x,0)  20  40 x
                       t                  2t
           u(0,t)  20e , u( 1 , t )  60e
• c = 0.5, h = 0.25, k = 0.05

      2
20e   -t                                              60e -2t
      1

      0
           0        1        2         3          4
                         20 + 40 x
                                   Example
• Explicit Euler method
                      ck ( 0.5 )( 0.05 )
                    r 2            2
                                          0.4
                      h     ( 0.25 )
             ui , j  1  ru i  1 , j  ( 1  2 r )ui , j  ru i  1 , j
                        0.4 ui 1 , j  0.2 ui , j  0.4 ui  1 , j
• First step: t = 0.05
 u0 ,1     20e 0.05  19.02458849

 u 1 ,1    0.4 u0 ,0  0.2 u1 ,0  0.4 u 2 ,0  0.4( 20 )  0.2( 30 )  0.4 ( 40 )  30

 u 2 ,1    0.4 u1 ,0  0.2 u 2 ,0  0.4 u 3 ,0  0.4 ( 30 )  0.2( 40 )  0.4( 50 )  40
u          0.4 u 2 ,0  0.2 u 3 ,0  0.4 u4 ,0  0.4 ( 40 )  0.2( 50 )  0.4 ( 60 )  50
 3 ,1
u          60e 0.10  54.29024508
 4 ,1
• Second step: t = 0.10
 u0 , 2    20e 0.10  18.09674836

 u1 , 2    0.4 u0 ,1  0.2 u1 ,1  0.4 u2 ,1
           0.4( 19.02458849 )  0.2( 30 )  0.4( 40 )  29.6098354


u 2 ,2     0.4 u1 ,1  0.2 u2 ,1  0.4 u3 ,1  0.4( 30 )  0.2( 40 )  0.4( 50 )  40
u          0.4 u2 ,1  0.2 u3 ,1  0.4 u4 ,1
 3 ,2
           0.4( 40 )  0.2( 50 )  0.4( 54.2924508 )  47.71609803

u4 ,2
           60e 0.20  49.12384518

                             29.61          40         47.72
20e -t                                                               60e -2t
                             30             40         50


              0            1              2           3             4
                                    20 + 40 x
Heat Equation: Time-dependent BCs

      r = 0.4
       Numerical Stability
• Stability for Explicit Euler Method
• It can be shown by Von Neumann
  analysis that
         1             1x   2
      r      or  t 
         2             2 c

• Switch to Implicit method to avoid
  instability
Explicit Euler Method: Stability

    r=1




                       Unstable !!
                  Implicit Euler method
                       Unconditionally Stable
      10

          9

          8

          7




u(0, t)
          6

          5
                                                              u(a, t)
= g1(t)   4                                                   = g2(t)
          3

          2

          1

          0
              0    1     2   3   4   5   6   7   8   9   10



          Initial conditions : u(x,0) = f(x)
                    Implicit Method
 • Finite-difference
                                          tj+1
                                                 t (i-1,j+1) (i,j+1) (i+1,j+1)
 T(x,t)
                   t                       tj
                                                               (i,j)




                                          x
                                                    xi-1      xi xi+1            x
Forward-difference u  1 ( u
                    t        i , j  1  ui , j )
                                  k
Central-difference         c
 at time level j+1 cu xx  2 ( ui  1 , j  1  2 ui , j  1  ui  1 , j  1 )
                                      h
             Implicit Euler Method
• Implicit Euler method for heat equation
   1                            c
     ( ui , j  1  ui , j )  2 ( ui  1 , j  1  2 u i , j  1  ui  1 , j  1 )
   k                           h
        ru i 1 , j  1  ( 1  2 r )ui , j  1  ru i  1 , j  1  ui , j

• Tridiagonal matrix (Thomas algorithm)
1  2 r   r                                  u1 , j  1   u1 , j  ru0 , j  1 
 r     1  2r   r                          u                      u2 , j          
                                              2 , j 1  
                                                                                      
                                                                                         
          r   1  2r           r            u3 , j  1           u3 , j          
                                            
                                     r                                          
                                                                                      

                                 r 1  2 r  un 1 , j  1  un1 , j  ru n , j  1 
                                                                                     

• Unconditionally stable
Implicit Euler Method

r=2




             Unconditionally stable
Example: Implicit Euler Method
  • Heat Equation (Parabolic PDE)
                   ut  cu xx ; 0  x  1
             u(x,0)  20  40 x
                         t                  2t
             u(0,t)  20e , u( 1 , t )  60e
  • c = 0.5, h = 0.25, k = 0.1

        1
  20e   -t                                              60e -2t

        0
             0        1        2         3          4
                           20 + 40 x
                             Example
• Implicit Euler method
                 ck ( 0.5 )( 0.10 )
               r 2            2
                                     0.8
                 h     ( 0.25 )
 (  r )ui 1 , j  1  ( 1  2 r )ui , j  1  (  r )ui  1 , j  1  ui , j
 ( 0.8 )ui 1 , j  1  ( 2.6 )ui , j  1  ( 0.8 )ui  1 , j  1  ui , j


 1  2 r         r              0   u1, 1  u1,0  ru0 , 1 
  r           1  2r                  u   u
                                  r   2 , 1   2 ,0
                                                                       
                                                                      
  0
                  r          1  2 r  u 3 , 1  u 3 ,0  ru 4 , 1 
                                                                   
 • Solve the tridiagonal matrix
   2.6  0.8    0   u1 ,1   30  0.8( 20e 0.1 )
    0.8 2.6  0.8  u   40                     
                     2 ,1                       
   0
          0.8 2.6  u3 ,1  50  0.8( 60e 0.2 )
                                                 
     u1 ,1   28.95515793
                         
   u2 ,1    38.50751457
    u  46.19426454
     3 ,1                

    1           28.96    38.51     46.19

20e -t                                         60e -2t

     0
         0     1        2         3        4
                    20 + 40 x
                   Crank-Nicolson method
   Implicit Euler method : first-order in time
    Crank-Nicolson : second-order in time
          10

          9

          8



u(0, t)
          7
                                                              u(a, t)
          6


= g1(t)   5
                                                              = g2(t)
          4

          3

          2

          1

          0
               0     1   2   3   4   5   6   7   8   9   10


                   Initial conditions : u(x,0) = f(x)
                   Crank-Nicolson Method
      • Crank-Nicolson method for heat equation
      • Average between two time levels
         1                            c
           ( ui , j  1  ui , j )     2
                                          ( ui  1 , j  2 ui , j  ui  1 , j )
         k                           2h
                                      c
                                    2 ( ui  1 , j  1  2 ui , j  1  ui  1 , j  1 )
                                     2h

      • Tridiagonal matrix
    r                                      r                r                             r
     ui 1 , j  1  ( 1  r )ui , j  1  ui  1 , j  1  ui 1 , j  ( 1  r )ui , j  ui  1 , j
    2                                      2                2                             2

      • Unconditionally stable (neutrally stable)
      • Oscillation may occur
     General Two-Level Method
• General two-stage method for heat equation
1                           c
  ( ui , j  1  ui , j )  2 ( ui  1 , j  2 ui , j  ui  1 , j )
k                           h
                            c( 1   )
                                 2
                                         ( ui  1 , j  1  2 ui , j  1  ui  1 , j  1 )
                                h

• Weighted-average of spatial derivatives
  between two time levels n and n+1

    λ  0 : im plicit Euler schem e
   
    λ  1 : explicit Euler schem e
    λ  1/2 : Crank - Nicolson schem e
   
Example: Crank-Nicolson Method
   • Heat Equation (Parabolic PDE)
                   ut  cu xx ; 0  x  1
             u(x,0)  20  40 x
                         t                  2t
             u(0,t)  20e , u( 1 , t )  60e
   • c = 0.5, h = 0.25, k = 0.1

        1
  20e   -t                                              60e -2t

        0
             0        1        2         3          4
                           20 + 40 x
                                  Example
• Crank-Nicolson method
                                       ck ( 0.5 )( 0.10 )
                                 r      2
                                                             2
                                                                      0.8
                                       h             ( 0.25 )
  r                                   r                  r                          r
 ui  1, j  1  (1  r)ui, j  1  ui  1, j  1  ui  1, j  (1  r)ui, j  ui  1, j
  2                                   2                  2                          2
  0.4ui  1, j  1  1.8ui, j  1  0.4ui  1, j  1  0.4ui  1, j  0.2ui, j  0.4ui  1, j


• Tridiagonal matrix (r = 0.8)
         r                   r                          r       r      
 1 r        0                  u0 ,0  ( 1  r )u1 ,0  u2 ,0  u0 ,1 
 r
          2
                   u 1 ,1   2
                               r
                                                           2       2
                                                                          
                r                                        r              
     1  r   u2 ,1    u1 ,0  ( 1  r )u2 ,0  u3 ,0              
   2           2              2                         2
 0                  u 3 ,1   r
                              
                                                                          
             1  r
          r                                                r       r
                                u2 ,0  ( 1  r )u3 ,0  u4 ,0  u4 ,1 
         2                   2
                                                          2       2      
                                                                          
 • Solve the tridiagonal matrix
 1.8  0.4    0   u1 ,1  0.4( 20 )  0.2( 30 )  0.4( 40 )  0.4( 20e 0.1 )
  0.4 1.8  0.4  u   0.4( 30 )  0.2( 40 )  0.4( 50 )                     
                   2 ,1                                                      
 0
        0.4 1.8  u3 ,1  0.4( 40 )  0.2( 50 )  0.4( 60 )  0.4( 60e 0.2 ) 
                                                                              
   37.23869934                 u1 ,1   29.42144598
                                                   
 40                         u2 ,1    39.29975855
  69.64953807                 u  47.42746748
                               3 ,1                




       1                29.42        39.30         47.43

20e -t                                                             60e -2t

       0
           0          1            2             3             4
                              20 + 40 x
Implicit Euler method

r=2




              Unconditionally stable
Heat Equation with Insulated
        Boundary
 • No heat flux at x = 0 and x = a
      x=0                                 x=a


     ux(0,t) = 0                        ux(a,t) = 0

       u      2u
          c       ; 0  x  a, 0  t  T
       t     x 2

        u(x,0)  f(x),    0 xa
       u x (0,t)  0
                         ,   0 t T
       u x ( a , t )  0
           Insulated Boundary
• No heat flux at x = a
                   un  1 , j  un  1 , j
  ux ( a , t )                               0  un  1  un  1
                     x n 1  x n 1
                                               ux(a,t)=0

                                             xn-1    xn     xn+1

                                                    x=a
 un , j  1  ru n1 , j  ( 1  2 r )un , j  ru n 1 , j
           2 ru n1 , j  ( 1  2 r )un , j