Numerical Solution of ODE with Euler & Runge-Kutta

Computer Methods in Chemical Engineering

The following dynamic equations govern a fermentor. 
    dx(t)
    ----- = u(s)*x(t) - D*x(t)     x(0)=1        ... biomass
     dt

    ds(t)             u(s)*x(t)
    ----- = D*[sf-s]- ---------    s(0)=10       ... substrate
      dt                 Y(s)

    where

                1*s
       u(s) = ------            ... specific growth rate
               1 + s

       Y(s) = 0.5 + 0*s         ... yield coefficient
Operating conditions: substrate feed concentration sf=10, dilution rate D=0.1.

Euler's Method for Numerical Integration of ODEs.

Find the values of x and s at t=0.1 in one step (i.e., step size = 0.1). 
With sf=10 and dilution rate D=0.1, the dynamic equations are:
    dx/dt = fx(x,s) = (u(s)-0.1)*x
    ds/dt = fs(x,s) = 0.1*(10-s) - 2*u(s)*x
Start with x=1 and s=10
    u = 10/(1+10) = 0.909
Slopes at x=1 & s=10
    dx/dt = (0.090-0.1)*1           =  0.809
    ds/dt = 0.1*(10-10) - 2*0.909*1 = -1.818
With h=0.1
    x (at t=0.1) = x(0) + dx/dt*h =  1 +  0.809*0.1 = 1.0809
    s (at t=0.1) = s(0) + ds/dt*h = 10 + -1.818*0.1 = 9.8182

Runge-Kutta's Method for Numerical Integration of ODEs.

In the following, we repeat with Runge-Kutta's 4th-Order Method. 
  1st evaluation is at (from Problem 1):
    x1 =  1
    s1 = 10
    u  = 10/(1+10) = 0.909
  1st slope for x and s (from Problem 1) ...
    kx1= fx(x1,s1) =  0.809
    ks1= fs(x1,s1) = -1.818
  2nd evaluation is at:
    x2 =  1 + kx1*h/2 = 1 + 0.809*0.05 = 1.0405
    s2 = 10 + ks1*h/2 = 10 -1.818*0.05 = 9.9091
    u  = 9.9091/(1+9.9091) = 0.9083
  2nd slope for x and s ...
    kx2= fx(x2,s2) = (0.9083-0.1)*1.0405               =  0.8411
    ks2= fs(x2,s2) = 0.1*(10-9.9091) - 2*0.9083*1.0405 = -1.8811
  3rd evaluation is at:
    x3 =  1 + kx2*h/2 = 1 + 0.8411*0.05 = 1.0421
    s3 = 10 + ks2*h/2 = 10 -1.8811*0.05 = 9.9059
    u  = 9.9059/(1+9.9059) = 0.9083
  3rd slope for x and s ...
    kx3= fx(x3,s3) = (0.9083-0.1)*1.0421               =  0.8423
    ks3= fs(x3,s3) = 0.1*(10-9.9059) - 2*0.9083*1.0421 = -1.8837
  4th evaluation is at:
    x4 =  1 + kx3*h = 1 + 0.8423*0.1 = 1.0842
    s4 = 10 + ks3*h = 10 -1.8837*0.1 = 9.8116
    u  = 9.8116/(1+9.8116) = 0.9075
  4th slope for x and s ...
    kx4= fx(x4,s4) = (0.9075-0.1)*1.0842               =  0.8755
    ks4= fs(x4,s4) = 0.1*(10-9.8116) - 2*0.9075*1.0842 = -1.9490
  Weighted average slope for x and s ...
    kxave = (kx1+2*kx2+2*kx3+kx4)/6 =  (0.809+2*0.8411+2*0.8423+0.8755)/6 =   5.0513/6 =  0.8419
    ksave = (ks1+2*ks2+2*ks3+ks4)/6 = (-1.818-2*1.8811-2*1.8837-1.9490)/6 = -11.2966/6 = -1.8828
  Update the values of x and s at t=h=0.1
    x = x + kxave*h =  1 +  0.8419*0.1 = 1.0840
    s = s + ksave*h = 10 + -1.8828*0.1 = 9.8117

Return to Prof. Nam Sun Wang's Home Page
Return to Computer Methods in Chemical Engineering (ENCH250)

Computer Methods in Chemical Engineering -- Numerical Solution of ODE
Forward comments to:
Nam Sun Wang
Department of Chemical & Biomolecular Engineering
University of Maryland
College Park, MD 20742-2111
301-405-1910 (voice)
301-314-9126 (FAX)
e-mail: nsw@umd.edu ©1996-2006 by Nam Sun Wang
UMCP logo