Transcript Slide 1

ORDINARY DIFFERENTIAL
EQUATIONS
(ODE)
Differential Equations
Heat transfer
d (  C pT )
qz
 
A
dz
(4.1)
ODE
Mass transfer
J AZ   DAB
dC A
dz
(4.2)
Conservation of momentum, thermal energy or
mass

 2
 2  R
t
z
(4.3)
PDE
ODE
independent
Definition


f t ,  (t ),  ' (t ),...,  ( n) (t )  0
Example
dependent
  2e     t
'''
(4.4)
t
''
'
4
A 3rd order differential equation for  = (t)
Solution dQ
 kQ
dt
Q(t )  cekt ,  t  
(4.5)
(4.6)
(4.7)
Important Issues
1. Existence of a solution
2. Uniqueness of the solution
3. How to determine a solution
Linear Equation (1)
an (t )  ( n )  an 1 (t )  ( n 1)  ...  a2 (t )  '  a1 (t )  ( 0)  a0 (t )  g (t )
a2 (t )  'a1 (t )   a0 (t )  g (t )
1.
2.
(4.8)
(4.9)
Rewrite 4.9
g (t )  a0 (t )
a1 (t )
 '

; a2 (t )  0
a2 (t )
a2 (t )
for all t
(4.10)
Determine
 a1 (t ) 
m (t )  exp  
dt 
 a2 (t ) 
where m(t) is called an integrating factor
(4.11)
Linear Equation (2)
3.
Multiply both sides of equation 4.10 by m(t)

g (t )  a0 (t )
a1 (t ) 

'


m
(
t
)

m (t )


a2 (t ) 
a2 (t )

(4.12)
Observe that the left-hand side of eqn 4.12 can be written
as
d
m (t )
dt
or
 a1 (t ) 
 a1 (t )  d
a1 (t )
 ' exp  
dt   
exp  
dt   m (t )
a2 (t )
 a2 (t ) 
 a2 (t )  dt
(4.13)
Linear Equation (3)
Equation (4.12) can be rephrase as:
 a1 (t )  g (t )  a0 (t )
 a1 (t ) 
d 


dt  
exp  
dt 
  exp  
dt 
a2 (t )
 a2 (t ) 
 a2 (t ) 
4.
(4.14)
Integrate both sides of Equation (4.14) with respect to
the independent variable:
 a1 (t ) 
 a (t ) 
g (t )  a0 (t )
dt   
exp   1 dt dt  c
a2 (t )
 a2 (t ) 
 a2 (t ) 
 exp  

 a1 (t ) 

a1 (t )  g (t )  a0 (t )
a1 (t ) 
 (t )  exp   
dt 
exp  
dt dt  c exp   
dt 
a2 (t )
 a2 (t ) 
 a2 (t ) 
 a2 (t ) 
where c is the constant of integration
(4.15)
Example 1
Water containing 0.5 kg of salt per liter is poured into a tank at a
rate of 2 l/min, and the well-stirred mixture leaves at the same rate.
After 10 minutes, the process is stopped and fresh water is poured
into the tank at a rate of 2 l/min, with the new mixture leaving at 2
l/min. Determine the amount (kg) of salt in the tank at the end of 20
minutes if there were 100 liters of pure water initially in the tank.
2 l/min
½ kg salt/l
CA
2 l/min, CA (l/min)
Solution
Example 2
Consider a tank with a 500 l capacity that initially
contains 200 l of water with 100 kg of salt in solution.
Water containing 1 kg of salt/l is entering at a rate of 3
l/min, and the mixture is allowed to flow out of the tank
at a rate of 2 l/min. Determine the amount (kg) of salt in
the tank at any time prior to the instant when the solution
begins to overflow. Determine the concentration (kg/l) of
salt in the tank when it is at the point of overflowing.
Compare this concentration with the theoretical limiting
concentration if the tank had infinite capacity.
Solution
THEOREM
If the functions p and g are continuous on an open interval
 < x < b containing the point x = x0, then there exists a
unique function y = f(x) that satisfies the differential
equation
y’ + p(x)y = g(x)
for  < x < b , and that also satisfies the initial condition
y(x0) = y0
where y0 is an arbitrary prescribed initial value.
Higher ODE Reduces to 1st Order
d2y
dy

q
(
x
)
 r ( x)
2
dx
dx
dy
Define z  , we have
dx
 dy
 dx  z

 dz  r ( x)  q( x) z
 dx
y '''( x)  y ( x) y ''( x )  2( y '( x )) 2  y ( x )  0
Define y1  y,









y2  y ',
dy1
 y2
dx
dy2
 y3
dx
dy3
  y1 y3  2( y2 ) 2  y1
dx
In general, it is sufficient to solve first-order ordinary differential
equations of the form
dyi
 fi ( x, y1 ,
dx
y3  y '', we have
, y N ), i  1, 2,
,N
Nonlinear equations can be reduced to linear
ones by a substitution. Example:
y’ + p(x)y = q(x)yn
and if n 0,1 then
n(x) = y1-n(x)
(4.16)
(4.17)
reduces the above equation to a linear equation.
Example 3
Suppose that in a certain autocatalytic chemical
reaction a compound A reacts to form a
compound B. Further, suppose that the initial
concentration of A is CA0 and that CB(t) is the
concentration of B at time t. Then CA0 – CB (t) is
the concentration of A at time t. Determine CB(t)
if CB(0) = CB0.
Solution
NONLINEAR ORDINARY
DIFFERENTIAL EQUATIONS
NONLINEAR EQUATIONS
Rewrite as
d
 f (t ,  )
dt
d
M (t ,  )  N (t ,  )
0
dt
If M is a function of t only, and N is a function of  only, then
d
M (t )  N (  )
0
dt
N (  )d   M (t )dt
Separable
NONLINEAR EQUATIONS
Consider
dC B
 kCB (C A0  C B )
dt
subject to
CB (0)  CB 0
Then, it is separable and results in:
dC B
 kdt
C B (C A0  C B )
(4.16)
Simplifying left-hand side; 1st consider the fraction
1
C B (C A0  C B )


CB


C A0  C B
where  and  are constants to be determined. Then:
 (C A0  CB )  CB  1
If we put
CB  0 : C A0  1
then

1
C A0
(4.17)
If we put
CB  C A0 : C A0  1
then

1
C A0
Rewrite equation 4.17
1
1
C A0
C A0
1


CB (C A0  CB ) CB C A0  CB
And equation 4.16 becomes
1
C A0
 1

1
 
 dCB  kdt
 CB C A0  CB 
which integrates to:
 CB 


 C A0  C B 
1
C A0
 m1 exp( kt)
where m1 is an arbitrary constant to be determined with the
given initial condition. @ t = 0, CB = CB0, then
 CB 


C

C
B
 A0
1
C A0
 m1
C B 0C A0
 C B (t ) 
(C A0  C B 0 ) exp( kCA0t )  C B 0
Example of Problem Setup
Consider the continuous extraction of benzoic acid from
a mixture of benzoic acid and toluene, using water as the
extracting solvent. Both streams are fed into a tank
where they are stirred efficiently and the mixture is then
pumped into a second tank where it is allowed to settle
into two layers. The upper organic phase and the lower
aqueous phase are removed separately, and the
problem is to determine what proportion of the acid has
passed into the solvent phase.
Example (cont…)
List of assumptions
1.
2.
3.
4.
5.
6.
7.
8.
Combine the two tanks into a single stage
Express stream-flow rates on solute-free basis
Steady flowrate for each phase
Toluene and water are immiscible
Feed concentration is constant
Mixing is efficient, the two streams leaving the stage are in
equilibrium with each other given by y = mx
Composition stream leaving is the same with the composition
in the stage
The stage initially contains V1 liter toluene, V2 liter water and
no benzoic acid
Problem 1
Consider an engine that generates heat at a rate of 8,530 Btu/min.
Suppose this engine is cooled with air, and the air in the engine
housing is circulated rapidly enough so that the air temperature can be
assumed uniform and is the same as that of the outlet air. The air is fed
to the housing at 6lb-mole/min and 65oF. Also, an average of 0.20 lbmole of air is contained within the engine housing and its temperature
variation can be neglected. If heat is lost from the housing to its
surroundings at a rate of Ql(Btu/min) = 33.0(T-65oF) and the engine is
started with the inside air temperature equal to 65oF.
1. Derive a differential equation for the variation of the outlet
temperature with time.
2. Calculate the steady state air temperature if the engine runs
continuously for indefinite period of time, using Cv = 5.00 Btu/lbmole oF.
Problem 2
A liquid-phase chemical reaction with stoichiometry A  B takes
place in a semi-batch reactor. The rate of consumption of A per unit
volume of the reactor is given by the first order rate expression
rA (mol/liter.s) = kCA
where CA (mol/liter) is the reactant concentration. The tank is initially
empty. At time t=0, a solution containing A at a concentration
CA0(mol/liter) is fed to the tank at a steady rate f(liters/s). Develop
differential balances on the total volume of the tank contents, V, and
on the moles of A in the tank, nA .
Solving ODEs using Numerical Methods
Initial Value Problem (IVP)
y’’ = -yx
y(0) = 2,
y’(0) = 1
Boundary Value Problem (BVP)
y” = -yx
y(0) = 2,
y’(1) = 1
General Procedure
Re-write the dy and dx terms as Δy and Δx and multiply
by Δx
dyi ( x)
 f i ( x, y1, y2 ,... y N )
for equations i  1,..., N
dx
Literally doing this is Euler’s method
dy( x)
 f ( x, y )
dx
y ( x)
 f ( x, y )
x
yi 1  yi  f ( xi , yi )x
Tank mixing problem
dc
V

(cin  c)
dt Vtank
 V

ci 1  ci  
(cin  ci ) t
Vtank

Mixing tank
t
Error Et
at t=600
300
1.4
150
0.61
100
0.39
50
0.19
30
0.11
15
0.055
10
0.036
5
0.018
3
0.011
Improvements to Euler’s Method
Euler
dy( x)
 f ( x, y )
dx
yi 1  yi  f ( xi , yi )x
Heun’s method (predictor-corrector)
Procedure
calc yi+1 with Euler (predictor)
calc slope at yi+1
calc average slope
use this slope to calc new yi+1 (corrector)
Heun example
dy
y
dx
Analytical
y (0)  0.5
1.5
at x  0.1
y ( x)  y0e x so y (1)  0.5e1  1.36
1
y
Solve
0.5
0
Numerical with Δx  1
0
x
y
 y0  (0.5) 
x
predicted location yi 1  yi  f ( x, y )x 
predictor (Euler)
slope
corrector : slope at predicted location
y
 y1 
x
average slope
updated location
yi 1  yi 
x 
1
Midpoint Method
1
y
Use Euler to calculate a midpoint location
1.5
evaluate slope y’ at the midpoint
0.5
0
use that to calculate full step location
y   f ( x, y )
x
yi 1 / 2  yi  f ( xi , yi )
2
yi 1  yi  f ( xi 1 / 2 , yi 1 / 2 )x
0
x
1
Runge-Kutta
dy
 f ( x, y )
dx
xh

y ( x  h)  y ( x ) 
 f ( x, y) dx
x
Could use higher order (polynomia l) fit for integral
based on values of f evaluated at locations between x  x  h
Or use Gaussian quadrature to evaluate the integral
General form
yi 1  yi  f h
f is increment function (aka a slope)
R-K – General form
General form
yi 1  yi  f h
Write f as :
f  a1k1  a2 k 2    an k n
where
a1 , a2 ,  , an
constants
k1  f ( xi , yi )
k 2  f ( xi  p1h ,
yi  q11k1h)
k3  f ( xi  p2 h ,
yi  q21k1h  q22 k 2 h)

k n  f ( xi  pn 1h ,
yi  qn 1,1k1h  qn 1, 2 k 2 h  ...  qn 1, n 1k n 1h)
R-K – 1st Order Form
General form
yi 1  yi  f h
f  a1k1
where
a1  1 constant
k1  f ( xi , yi )
yi 1  yi  1 f ( xi , yi )h
R-K – 2nd Order Form
yi 1  yi  a1k1  a 2 k 2  h
y(x)
k1  f ( xi , yi )
k 2  f ( xi  p1h ,
a1,a2 ,p1,q11
yi  q11k1h)
constants
a1  a2  1
1
a 2 p1   a2 q11
2
xi
xi+1
x
RK2 – Options
yi 1  yi  a1k1  a2 k 2  h
k1  f ( xi , yi ) ; k 2  f ( xi  p1h, yi  q11k1h)
1
a1  a2  1 ; a2 p1   a2 q11
2
a2  0

a1  0

y(x)
xi
x
xi+1
RK2 – Options
yi 1  yi  a1k1  a2 k 2  h
k1  f ( xi , yi ) ; k 2  f ( xi  p1h, yi  q11k1h)
a1  a2  1 ; a2 p1 
1
a2 
2

2
3

a2 
y(x)
1
 a2 q11
2
xi
xi+1
x
xi
xi+1
x
y(x)
R-K – 2nd Order Form
yi 1  yi  a1k1  a2 k 2  h
k1  f ( xi , yi )
k 2  f ( xi  ph, yi  qk1h)  f  ph
a1,a2 ,p1,q1
constants
f
f
 qfh  O(h 2 )
x
y
RK – 3rd Order Form
1
yi 1  yi  k1  4k 2  k3  h
6
y(x)
k1  f ( xi , yi )
1
1
k 2  f ( xi  h , yi  k1h)
2
2
k3  f ( xi  h , yi  k1h  2k 2 h)
xi
xi+1
x
RK – 4th Order
1
yi 1  yi  k1  2k 2  2k3  k 4  h
6
y(x)
k1  f ( xi , yi )
1
1
k 2  f ( xi  h , yi  k1h)
2
2
1
1
k3  f ( xi  h , yi  k 2 h)
2
2
k3  f ( xi  h , yi  k3h)
xi
xi+1
x
Example y΄=x+y, y(0)=0
x
yo
k1=fi
k2=f(x+h/2,y+h
/2k1)
0
0
0
0.1
0.11
0.222
0.02140
0.2
0.021
4
0.221
0.344
0.356
0.493
0.0918
0.4
0.092
0.492
0.641
0.656
0.823
0.2221
0.6
0.222
0.822
1.004
1.023
1.227
0.4255
0.8
0.426
1.226
1.448
1.470
1.720
0.718
1
0.718
1.718
1.990
2.017
2.322
1.120
1.2
1.120
2.320
2.5
2.652
2.685
3.057
1.655
1.4
1.655
3.055
2
3.461
3.501
3.955
Euler
2.353
1.6
2.353
3.953
1.5 4.448
4.498
5.052
analytical
3.250
1.8
3.250
5.050
1 5.654
5.715
6.393
RK4
4.389
3.5
y
3
k3=f(x+h/2,y+h
/2k2)
k4=f(x+h,y+
hk3)
yn=yo+1/6(k1+2k2+2k
3+k4)h
0.5
0
0
0.5
1
x
1.5
2