ChE306: Heat and Mass Transfer

1
Chapter 5
Numerical Solution of Ordinary Differential
Equations
2
Ordinary Differential Equations (ODEs)

• Classifications of ODEs
• Order (by highest derivative)
• Linearity (nonlinear if it contains products of
  the dependent variable or its derivatives, or
  both)
• Boundary conditions (initial vs. boundary value)
• Homogeneous vs. non-homogeneous R(x)0
• Variable vs. constant coefficients

3
Transformation to Canonical Form

• n simultaneous 1st order ODEs in canonical form

given initial conditions
solution has form
Expressed in condensed matrix form
4
Transformation to Canonical Form
define transformations

• nth order ODE of the form

substituting
said to be autonomous if RHS ? f(x)
5
Transform to Matrix Form
autonomous example
transforms
6
Transform to Matrix Form
nonautonomous example
transforms
transformed (autonomous)
7
Linear Ordinary Differential Eqns
single ODE set of ODEs
analogous
analogy
8
Linear Ordinary Differential Eqns
matrix exponential function eAt can be expanded as
and can be shown to be a solution of
express in terms of eigenvalues eigen- vectors
of A to avoid evaluation of ? series
9
Linear Ordinary Differential Eqns
For nonsingular matrix A of order n, there are n
eigenvectors and n nonzero eigenvalues
10
Linear Ordinary Differential Eqns
Manipulating the compact form of the eigenvector
and eigenvalue relationship
or by raising to any power n
11
Linear Ordinary Differential Eqns
Starting with the matrix exponential function,
replace each term with the equivalent form of An
solution is thus
12
Example
Series equilibrium reactions
solution must be
reduces to
13
LinearODE.m code
case 1 Matrix exponential method for k
1nt if t(k) gt 0 y(,k) expm(At(k))y0
else y(,k) y0 end end
case 2 Eigenvector method X,D eig(A)
Eigen vectors/values IX inv(X) e_lambda_t
zeros(nA,nA,nt) Building the matrix
exp(LAMBDA.t) for k 1nA e_lambda_t(k,k,)
exp(D(k,k) t) end Solving the set of
equations for k 1nt if t(k) gt 0 y(,k)
X e_lambda_t(,,k) IX y0 else
y(,k) y0 end end
14
Nonlinear ODE Initial Value Problems

• Consider the set of ODEs in canonical form
• Treat 1 ODE as follows

treat by finite differences
15
Euler Method
16
Euler Method
17
Euler Method
18
Euler Method
19
Euler Method
20
Euler Method
21
Explicit Euler Method

• Forward marching by finite difference

EQN 3.53
next value obtained from previous value plus
tangential direction step of width h
22
Implicit Euler (Backward) Method

• Accuracy can be improved by combining forward and
  backward differences
• Forward marching by backward difference

EQN 3.32
function is evaluated at an unknown value, (xi1
,yi1)
23
Modified Euler (Predictor/Corrector)

• must be solved as set of simultaneous algebraic
  equations
• Linear ? Gauss Elimination
• Nonlinear ? Newton's Method
• Modified Euler ? Predictor/Corrector (also known
  as Crank-Nicolson)
• Apply Explicit Euler to predict
• Apply Implicit Euler using predicted value for f
• Expand as

term is in ? - in ? ? O(h3)
24
Modified Euler (Predictor/Corrector)
25
Crank-Nicolson Method

• Generalize Predictor/Corrector equation as
• Nature of weighting functions and positions
  dictated by number of terms retained.
• Forms the basis for a series of integration
  formulas with increasingly higher accuracy for
  ordinary differential equations.

26
Runge-Kutta Method

• Most widely used ODE integration method
• Single step, self-starting method

27
Runge-Kutta Method

• In compact notation
• The value of m is fixed by retaining m1 terms of
  the infinite series

28
Derivation of Runge-Kutta Method

• Step 1
• Chose m (fixes accuracy)
• m 2 for second-order Runge-Kutta
• truncate series after m1 terms

29
Derivation of Runge-Kutta Method

• Step 2
• Replace each derivative of y by its equivalent in
  f, remembering f is a function of both x and y(x)

30
Derivation of Runge-Kutta Method

• Step 3
• Write weighted trajectories formula with m terms
• Step 4
• Expand the f function in a Taylor series
• combine

31
Derivation of Runge-Kutta Method

• Step 5
• In the 2nd order Runge-Kutta method, there are
  three equations in four unknowns. There is
  always degrees of freedom more than the order of
  the relationship, allowing some flexibility in
  selection.

(2) define
(1) Compare Taylor Series expansion to weighted
trajectories
(3) chose
(4) 2nd order Runge-Kutta is identical to
Crank-Nicholson
32
Adams Method

• Multiple step ? non-self starting
• Must use a self-starting method to initiate
• Pass quadratic polynomial through (xi-2, yi-2),
  (xi-1, yi-1), and (xi, yi), extrapolate to
  f(xi1, yi1)

33
Adams Method

• Chose a uniform step size, apply 2nd degree
  backward Gregory-Newton interpolating polynomial
• where
• upon integration

34
Adams Method

• Replace ? with expansions, rearrange
• To use requires two evaluations by Runge-Kutta
  before method can begin.

35
Adams-Moulton Predictor/Corrector

• Applies 3rd degree interpolating polynomial
• Interpolation performed with cubic Gregory-Newton
  over the range xi-2 to xi1

36
Example 5.3 - Non-isothermal PFR
mole balance
energy balance
rate law
heat of reaction
heat capacity
37
Simultaneous Differential Equations

• Expand to 4th order RK is easy to code

38
Nonlinear ODEs Boundary Value

• Canonical form
• r equations with initial conditions
• n-r equations with final conditions

39
Shooting Method

• Converts boundary value problem to an initial
  value problem
• Initial conditions guessed
• Final conditions calculated
• Calculated conditions compared to boundaries
• Guess conditions adjusted
• Procedure repeated until converged

40
Shooting Method
split boundary conditions
initial guess
desired
rearrange
expand
41
Shooting Method
expand
for convergence
which is of the form of Newton-Raphson
apply limit to expansion
recall definition of ?(?)
defined
substitute
42
Shooting Method
suggests an interative equation of the form
iterate until the condition is achieved
43
Shooting Method

• Generalized to a set of n simultaneous eqns

(1) guess n-r initial conditions
(2) integrate forward simultaneously
(3) evaluate J (Jacobian Matrix)
44
Shooting Method

• Generalized to a set of n simultaneous eqns

(4) correct initial condition guess
(5) iterate (repeat steps 2 - 4) until
45
Finite Difference Method

• Replace derivatives with difference equations

46
Finite Difference Method

• Divide integration into n segments of equal
  length, write difference eqns for i0,1,2,n-1
• 2n simultaneous nonlinear algebraic equations in
  (2n2) variables
• BCs provide equalities for remaining 2 variables
• Solve using Newton's Method
• More difficult to apply than Shooting Method
• Only use when Shooting Method is too unstable
• If DEs are linear, solution by Difference
  Equation Method is simple by matrix inversion or
  Gauss Elimination

47
Collocation Methods

• Again
• Now suppose solutions y1(x) and y2(x) can be
  approximated by polynomials (trial functions)

48
Collocation Methods

• Take derivatives of the trial functions and
  substitute into the differential equations
• Form residuals as
• Determine coefficients cm,i to minimize residuals

method of weighted residuals
49
Collocation Methods

• Collocation choses the weighting function as the
  Dirac Delta function (unit impulse)
• which has the property
• method of weighted residuals becomes

50
Collocation Methods

• This expression contains (2n2) undetermined
  coefficients, cm,i i 0, 1, , n m 1, 2,
  
• cm,i calculated by choosing (2n2) collocation
  points
• 2 collocation points fixed by boundary conditions

51
Collocation Methods

• Chose the remaining 2n internal collocation points

y1,f and y2,0 are unknown
52
Collocation Methods

• complete set of (2n2) simultaneous nonlinear
  equations in (2n2) unknowns
• Solution by Newton's method
• If collocation points are chosen at equidistant
  intervals, collocation method is equivalent to
  polynomial interpolation of equally spaced points
  and to the finite difference methods (not
  preferred)
• It is advantageous to locate the collocation
  points at the roots of the appropriate orthogonal
  polynomials

53
Orthogonal Collocation

• Chose trial functions y1(x) y2(x) as linear
  combinations

in general,
54
Orthogonal Collocation

• Chose ci,k such that
• Choosing Pi(x) as the Legendre Polynomials
  w(x)1interval is -1, 1, therefore transform
  as
• 2-point BV problem has (2n2) collocation points,
  zj j 0, 1, , n1 with 2 known boundary
  values z0 -1 and zn1 1

55
Orthogonal Collocation

• Location of the internal n collocation points z1
  to zn are determined from the roots of Pn(z) 0
• 
  
  must be determined so
  that the BCs are satisfied. Rewrite in notation
  as

d1 d2 matricies of unknown coefficients
56
Orthogonal Collocation

• The derivatives of y can be expressed in terms of
  d and z

57
Orthogonal Collocation

• The restatement of the boundary value problem
  becomes

with boundary conditions
This represents 2n4 simultaneous nonlinear
equations that can be solved by Newton's method.
58
Orthogonal Collocation

• Presentation in matrix form (still for 2 ODE
  system)

First and last equations are not included in
matrix set because they are established by a
boundary condition rather than a collocation
point.
59
Orthogonal Collocation

• Generalize to m simultaneous 1st order ODEs

dependent variables yij i 1,2,,m j
0,2,,n1 are evaluated from the simultaneous
solution of the following set of nonlinear
equations plus boundary conditions
60
Orthogonal Collocation

• For a second-order 2-point boundary value
  problem,
• Transform independent variable x ? z at n2
  points
• Derivatives are

61
Orthogonal Collocation

• In matrix form

62
Example 5.5 Orthogonal Collocation

• Optimal temp profile for batch penicillin
  fermentation
• cell mass production
• penicillin synthesis
• bi rate constants (gt0) f(?)
• y1 dimensionless conc of cell mass
• y2 dimensionless conc of penicillin
• yi(0) 0.03, 0, 0, 1
• ? temperature
• w 13.1, 0.005, 30, 0.94, 1.71, 20

Hamiltonian
necessary condition
63
Stability

• Inherent Stability
• Determined by the mathematical formulation of the
  problem. Dependent on Eigen values of the
  Jacobian matrix of the ODE system
• Numerical Stability
• Nature of error propagation in the numerical
  integration method. Behavior of error
  propagation dependent on values of the
  characteristic roots of the difference equations
  used to derive the solution.

64
Error

• Truncation error
• Function of the number of terms retained in the
  approximation of the infinite series expansion
• Can be reduced by retaining more terms or by
  reducing step size (h).
• Lower truncation error (higher accuracy) comes at
  the cost of increased computational requirements,
  which also leads to an accumulation of round-off
  error.
• Round-off error
• Related to the retention of a finite number of
  digits.
• Retention of more digits reduces round-off error.
• Error Accumulation and Propagation
• Accumulation of error grows exponentially or in
  an oscillatory pattern.
• Propagation of error causes deviation from
  correct solution.

65
Stability Error Propagation in Euler's

• Consider the initial value differential equation
• Apply explicit Euler, ignore error

analytical solution
?? is real and y0 is finite
Suggests solution is inherently stable
1st order, homogeneous difference equation
66
Stability Error Propagation in Euler's

• Identify characteristic equation
• n increases unbounded, therefore it's behavior is
  determined by value of (1?h). Numerical
  solution is absolutely stable if
• Because (1?h) is the root of the characteristic,
  an alternative definition of stability is

w/ initial cond gives solution
with root
requires
67
Stability Error Propagation in Euler's

• Step size ? limited for stable solution. h must
  be gt0.
• Explicit Euler is conditionally stable
• Any method with an infinite general stability
  boundary is considered unconditionally stable.
• ? may be complex, in which case solution is
  stable, converging with damped oscillations when
  moduli of the roots is ?? 1 (i.e., r ? 1).

A finite general stability boundary