Ordinary differential equations

1
Ordinary differential equations

• KE-42.4520
• Process modeling methods and tools
• Kari I. Keskinen
• 3.3.2008

2
Introduction

• Differential equations are used to describe
  phenomena that are in movement or change with
  time.
• Pre-examination material handled quite much the
  methods to solve differential equations
  analytically.
• The classification of differential equations is
  important in order to find a suitable solution
  strategy.
• In this lecture we handle mainly ordinary
  differential equations.

3
Plug flow reactor

• Material balance

4
Plug flow reactor

• Let us inspect a volume ?V, where the reaction
  rate can be assumed constant.

y?y
y
Fj(y?y)
?V
5
Plug flow reactor

• Reaction rate is constant in ?V

Steady state assumption, i.e. no accumulation
6
Plug flow reactor

• For the volume ?V we can write the material
  balance for component j.

A cross-sectional area
7
Plug flow reactor
8
Plug flow reactor

• Separation of variables

Integral form of the dimensioning equation.
9
Classification of differential equations

• A differential equation is an ordinary
  differential equation, ODE, when there are
  derivatives only with respect to one variable.
  Often this is a change with time.
• When the same equation contains derivates with
  respect to more than one variable, it is called
  partial differential equation, PDE. Often this is
  with respect to several spatial coordinates
  (x,y,z) and time (t).

10
ODEs 1

• Ordinary differential equations can be classified
  in two groups
• initial value problems the variable value at
  t0, must be known and also for higher order
  differential equations the values of the
  derivates up to degree of highest derivative
  minus one must be known at t0
• boundary value problems the variable value at
  two boundaries must be known.

11
ODEs 2

• Initial value ODE.
• Let x(t) and y(t) denote the number of rabbits
  and foxes at time moment t, respectively. The
  evolution of the rabbit and fox populations can
  then be described as a function of time

12
ODEs 3

• The example shows that a set of ODEs can be
  included in the model. These are still ODEs, as
  they both have a derivate only with respect to
  one variable.
• Generally, of set of ODEs can be expressed in a
  form of a vector equation
• where y is a vector and f is vector valued
  function.

13
ODEs 4

• Boundary value ODE.
• The solution for a set of n ODEs is sought after
• where vectors y and f contain n components. In
  addition the solution is bound by two distinct
  points (a and b) of the interval of interest
• Here the boundary condition is given by

14
ODEs 5

• The solution of boundary value ODE is much more
  difficult than solving initial value ODE.
• The degree of the highest derivative determines
  the order of the ODE.
• One classification of ODEs is
• linear ODEs
• non-linear ODEs
• Normally, the solution of linear systems is much
  easier than that of non-linear ODEs.

15
ODEs 6

• Linear ODE
• The coefficient a do not depend on the variable
  y. f(x) is a forcing function.
• Nonlinear ODE

16
ODEs 7

• The latter is a function that is homogenous of
  degree n satisfying for every value of ?
• and can be transformed into form (by
  substitution )

17
ODEs 8

• Higher order ODEs can be treated as follows
• This is a second order ODE. Let us denote y y1
  for the original variable and take into use a new
  variable y2 y1. Then we obtain a set of first
  order ODEs

18
ODEs 9

• By solving the original ODE as a set of first
  order ODEs one obtains in addition to the
  approximation y1 for the solution curve y(t) the
  approximation y2 for the derivative y(t) .
• Generally a differential equation of order n

19
ODEs 10

• can be transformed into a set of n coupled first
  order differential equations

This is the RHS referred later.
20
Solution methods

• Some solution methods were reviewed in the
  pre-examination material
• I-factor method for first order ODE
• manipulate non-linear equations into linear form
• no general methods are available for non-linear
  systems

21
Solution methods 2

• One can try to use the mathematical software in
  finding a solution to ODE problems. (Mathematica,
  Maxima, Maple etc.) Often, these dont help as
  there is no analytical solution available.
• Other methods are series solutions and the use of
  Laplace transforms. A few slides of these
• The main way to tackle ODEs and systems of them
  is to use numerical methods.

22
Series solutions
Pre-examination material showed that all
linear differential equations with constant
coefficients, like
sustained complementary solutions of the form
where the characteristic roots in this case
satisfy
.
The exponent function exp() can be expressed as
a series.
23
Series solutions 2
Based on this the solution of above
differential equation could have been tried as
Now we know that , but in general
case when analytical solution is not possible,
e.g. non-constant coefficients arise, like
,
this method proves useful.
24
Series solutions 3
The series usage has been developed into
method of Frobenius (see e.g. Rice Do, chapter
3.) Some special forms of differential equations
often arise in engineering and these has been
given names after persons Bessel, Laguerre,
Legendre, Chebysev. The solution of these
equations can be expressed by using special
functions by the names Bessel functions,
Laguerre polynomials, Legendre polynomial and
Chebysev polynomials. These are found in
handbooks in tabulated form.
25
Laplace transforms

• Using the complex variable theory the following
  integral transform has been defined
• s is a complex variable defined by
  , where s and ? are real variables and
  .
• The usage of this means that we map our problem
  for example from time variable t to an abstract
  Laplace space where we have the Laplace variable
  (often used are s and p).
• Then the equations are solved.

26
Laplace transforms 2

• Once the solution is available as a function of
  the Laplace variable, an inverse Laplace
  transform is carried out.
• The solution is normally such that a ready made
  tabulated inverse Laplace transform is not
  available. Then the solution must be manipulated
  for example by making it into a partial fraction
  expression. In order to do that quite often a
  high degree polynomial must be solved for roots.

27
Laplace transforms 3

• The elegant point of the Laplace transform method
  is that the initial values or boundary values are
  taken into account already at the beginning in
  the solution method (and not as the final point
  for finding a particular solution for the general
  solution).

28
Huokosdiffuusio
R
Matemaattinen malli

• An moolitase ?rssä (?V)
• sis. ulos muod. akkum.
• sis. WAr 4?r2?r
• ulos WAr 4?r2?r?r

CAs
r
r ?r
rm on keskisäde välillä r ja r?r
29
Huokosdiffuusio

• ?V approksimoidaan

?r ? 0
30
Huokosdiffuusio

• reaktio A ? B
• ? ekvimolaarinen vastavirtaus (lisäksi vakio T
  ja p)
• aiemmin todettu

• Yhtälö soveltuu myös laimeille seoksille.
• Sijoittamalla

31
Huokosdiffuusio

• Laadut

32
Huokosdiffuusio

• oletus

Differentiaaliyhtälö diffuusio reaktio
33
Huokosdiffuusio

• Reunaehdot
• CA äärellinen, kun r 0
• CA CAs, kun r R

34
Huokosdiffuusio

• ratkaistaan dimensiottomien muuttujien avulla

• uudet reunaehdot
• CA CAs r R

CA äärellinen r 0 ?
äärellinen ? 0
35
Huokosdiffuusio

• Differentiaaliyhtälö dimensiottomien muuttujien
  avulla

36
Huokosdiffuusio
37
Huokosdiffuusio
38
Huokosdiffuusio

• Kun Thiele-moduuli on suuri, diffuusio
  tavallisesti rajoittaa reaktion kokonaisnopeutta
  kun ?n on pieni, pintareaktio määrää nopeuden.
• Esimerkki ratkaisusta
• A ? B

ensimmäinen kertaluku (pintareaktio rajoittaa,
korkea T)
k1 m3/m2 s m/s
39
Huokosdiffuusio

• Differentiaaliyhtälö

? 1 ? 1
Reunaehdot
? äärellinen ? 0
Ratkaistaan uuden muuttujan avulla
y ? ? ? (1/?) y
40
Huokosdiffuusio reaktio
Sijoittamalla
41
Huokosdiffuusio reaktio

42
Hyperboliset funktiot
43
Huokosdiffuusio reaktio

• Reunaehdot
• ? 0 ? äärellinen
• cosh?1? ?1

sinh?1? ? 0
1/? ? ? ? äärellinen
A1 0
44
Huokosdiffuusio reaktio

• Toinen ehto
• ? 1 ? 1

? sijoittamalla A1 ja B1
45
Huokosdiffuusio reaktio

• Pienet ?1n arvot ? pintareaktio kontrolloi
  eli An konsentraatio on suuri huokosissa (ennen
  kuin A reagoi)

46
Non-isothermal porous catalyst particle

• For a non-isothermal porous catalyst particle
  the following equation has been derived.

47
Non-isothermal porous catalyst particle

• From the analogy of mass and heat transfer ?
  correspondingly from the energy balance

Boundary conditions
dT/dr 0 r 0 T Ts r R

• ? To be solved together numerically, because
  kn k0 e-E/RT

48
Non-isothermal porous catalyst particle

• The differential equations can be combined ?
  equation between temperature and concentration

kt effective heat transfer coefficient
49
Non-isothermal porous catalyst particle

• if kt ? f(T) and ?H ? f(T)
• First integration

Boundary conditions
? r 0 dCA/dr 0 and dT/dr 0
? A1 0
50
Non-isothermal porous catalyst particle
? r R CA CAs T Ts
De CAs (kt/?H) Ts A2
A2 De CAs - (kt /?H) Ts
Independent from the kinetic expression!
51
Non-isothermal porous catalyst particle,
temperature gradient

• Maximum ?T is obtained, when all A has been
• converted, i.e. CA 0

We can evaluate, whether ?T inside particle is
important.
If an exothermal reaction and non-isothermal case
? can be gt 1, because the external temperature
is lower than the temperature inside particle.
? Reaction rate inside particle higher that on
the particle surface.
52
Non-isothermal porous catalyst particle,
temperature gradient

• Define two parameters ? and ?

kt thermal conductivity
Possibility for several stationary states.
No multiple stationary states when
4(1 ?) gt ? ?
53
Internal effectiveness as a function of Thiele
modulus
54
Example 1

• Classify
• Linear Nonlinear
• Homogeneous Inhomogeneous
• Order
• Degree
• Forcing function

55
Example 1

• Classify
• Linear Nonlinear (if nltgt1)
• Homogeneous Inhomogeneous
• Order 2
• Degree 1
• Forcing function? 0

56
Example 2

• Classify
• Linear Nonlinear
• Homogeneous Inhomogeneous
• Order
• Degree
• Forcing function?

57
Example 2

• Classify
• Linear Nonlinear
• Homogeneous Inhomogeneous
• Order 1
• Degree 2
• Forcing function? x - 1

58
Example 3

• Classify
• Linear Nonlinear
• Homogeneous Inhomogeneous
• Order
• Degree
• Forcing function?

59
Example 3

• Classify
• This is actually a differential-algebraic
  equation as there is also an algebraic equation
  in addition to a differential equation.

60
Why numerical solutions?

• Sometimes an analytical solution can be found by
  using the mathematical solution methods or by
• Series solution methods, like Frobenius
• Special equation types of Bessel, Laguerre,
  Chebysev, Legendre
• Laplace transforms
• Symbolic mathematic programs
• Often the evaluation of a numerical value from an
  analytic solution might be painful numerically.

61
Stability of ODEs

• The second order differential equation
• has the solution y(t)e-t. If we change the
  starting value of the first initial value by e
  then we have new initial values
• By solving the problem with the new initial
  values we obtain the solution
• The RHS goes to infinity when time goes to
  infinity, no matter how small the positive e is.

62
Stability of ODEs 2

• The previous example is clearly unstable as an
  extremely small deviation in the starting values
  causes the divergence of the solution. The
  solution explodes.
• This is important in numerical computation as
  there are always limited accuracy and the
  rounding errors tend to cumulate and cause the
  same phenomena as the erroneous starting values.

63
Stability of ODEs 3

• The stability of a set of ODEs can be evaluated
  using the Jacobian matrix, J, of the solutions.
• The Jacobian is evaluated for the solution and if
  all the real parts of eigenvalues are negative,
  the set of ODEs is stable.

64
Solutions by Euler method
Model
Initial condition
Closure first order reaction
65
Separable form, analytical solution possible
66
Numerical solution Eulers method
Explicit Euler if
is calculated at tn
Implicit Euler if
is calculated at tn1
67
Numerical solution Eulers method
For the test problem
Explicit Euler
Implicit Euler
68
Numerical solution Eulers method
Implicit Euler
For this linear test problem we have an explicit
formulation
Usually implicit methods require solution of a
nonlinear set of equations at each time step
For example, pth order reaction
69
Example

• c01
• k0.1
• ?t1

70
Example

• c01
• k0.5
• ?t1

71
Example

• c01
• k1
• ?t1

72
Example

• c01
• k1.5
• ?t1

73
Example

• c01
• k2
• ?t1

74
Example

• c01
• k2.5
• ?t1

75
Shorter time steps

• c01
• k2.5
• ?t0.2

76
Other options

• Calculate

at average of tn and tn1
77
Crank-Nicolson method

• c01
• k0.5
• ?t1

78
Crank-Nicolson method

• c01
• k2.5
• ?t1

Stable, but oscillates!
79
Solving initial value ODEs

• The methods are
• Runge-Kutta
• Richardson extrapolation, especially
  Bulirsch-Stoer method
• Predictor-corrector methods, also known as
  multistep methods
• Runge-Kutta methods propagate a solution over an
  interval by combining the information from
  several Euler-style steps, and then using the
  information obtained to match a Taylor series
  expansion up to some higher order.

80
Solving initial value ODEs 2

• Richardson extrapolation uses the idea of
  extrapolating a computed result to the value that
  would have been obtained if the step size had
  been very much smaller than it actually was. In
  particular, extrapolation to zero step size is
  the desired goal. The first practical ODE
  integrator that implemented this idea was
  developed by Bulirsch and Stoer. Thus the
  commonly used name Bulirsch-Stoer methods.

81
Solving initial value ODEs 3

• Predictor-corrector methods or multistep methods
  store the solution along the way, and use those
  results to extrapolate the solution one step
  advanced. They then correct the extrapolation
  using derivative information at the new point.
  These are best for smooth functions.

82
Runge-Kutta

• As Euler method (
  advancing the solution from to
  ), is not very accurate, an improvement
  of it is that a trial step is taken to the
  midpoint of the interval.
• Then values of both x and y at that midpoint are
  used to compute the real step across the whole
  interval.

83
Runge-Kutta 2

• This is the second-order Runge-Kutta method. Can
  be extended to fourth-order Runge-Kutta method as
  follows

84
Runge-Kutta 3

• The 4th-order Runge-Kutta requires four
  evaluations of RHS while second-order Runge-Kutta
  requires less RHS evaluations.
• Provided that 4th-order Runge-Kutta has much more
  a longer step size than second-order Runge-Kutta,
  it is superior in performance.
• Still, higher order method does not always
  guarantee high accuracy.
• Using an adaptive step size in 4th-order
  Runge-Kutta improves the method. There are
  several variants of this class.

85
Richardson Extrapolation

• These methods are good for smooth functions.
  Bulirsch-Stoer method is believed to be the
  best-known way to obtain high accuracy solutions
  for ordinary differential equations with minimal
  computational effort. Three caveats
• For nonsmooth problem, use Runge-Kutta.
• If singular points inside the integration
  interval, then adaptive step size Runge-Kutta or
  special methods.
• If the function is smooth and the RHS expensive
  to evaluate, use predictor-corrector methods.

86
Bulirsch-Stoer method

• A modified midpoint method advances a vector of
  dependent variables y(x) from point x to a point
  x H by a sequence on n substeps of size h.
• hH/n
• The number of RHS evaluations required is n 1.
  The formulas for the method are
• Here the zs are intermediate approximations that
  march along in steps of h, while yn is the final
  approximation to y(xH).

87
Bulirsch-Stoer method 2

• Bulirsch-Stoer attempts to cross the interval H
  using the modified midpoint method with
  increasing value of n, the number of substeps.
  They used the sequence
• n 2, 4, 6, 8, 12, 16, 24, 32, 48, 64, 96,
• nj2nj-2
• After each successive n is tried, a polynomial
  extrapolation is tried. If errors are not
  satisfactory a higher n is attempted.
• If we go too far in increasing n, then the
  interval H is reduced. There also a limit for
  reduction of H.

88
Stiff sets of equations

• If we have more than one first-order ODE, the
  possibility of a stiff set of equations arises.
• Normally stiffness occurs when there are very
  different time scales involved in the system One
  phenomena is extremely fast and requires a very
  small time step, while some other phenomena are
  very slow and determine the total time required
  for the solution interval.

89
Stiff sets of equations 2

• An initial value problem is stiff on interval
  a,b if it is valid for all that
• Re(?i) lt 0, i 1, 2, , s,
• where ?i are the eigenvalues of the Jacobian
  matrix J determined at point (t,y(t)). Term S(t)
  is called the stiffness coefficient at point t.
• The definition means that a stabile set of ODEs
  is stiff when the difference between the largest
  and smallest real parts of the eigenvalues of the
  Jacobian is big.

90
Stiff sets of equations 3

• Consider
• With boundary conditions
• A solution is found

91
Stiff sets of equations 4

• The part e-1000x requires a very small step size
  in the solution algorithm.
• The part e-1000x is negligible when compared to
  the first part of the solution, as soon as one is
  away the origin.
• The accuracy requirement of the solution for
  reasonably large interval would allow much longer
  step based on the first term than the term
  e-1000x .

92
Stiff sets of equations 5

• In stiff problems the step size is largely
  determined by terms that are negligible and fast
  disappearing on the point of view of the whole
  problem.
• This kind of systems are often met when chemical
  reaction kinetics are modeled and some of the
  components (intermediate) react rapidly but the
  overall reactions are slow.

93
Stiff sets of equations 6

• The explicit (or forward) Euler type methods
  would diverge.
• Implicit methods has to be used, like backward
  Euler
• Methods used for stiff systems Rosenbrock and
  semi-implicit Euler. The best methods for stiff
  problems are multistep backward methods.

94
Multistep, Multivalue, and Predictor-Corrector
Methods

• To advance the solution of yf(x,y) from xn to x
  we have
• In a single-step method (like Runge-Kutta or
  Bulirsch-Stoer) the value yn1 at xn1 depends
  only on yn.
• In multistep methods f(x,y) is approximated by a
  polynomial passing through several points xn,
  xn-1, and possibly also through xn1.

95
Multistep, Multivalue, and Predictor-Corrector
Methods 2

• The result of evaluating previous integral at x
  xn1 is then of the form
• Here yn denotes f(x,y) and so on.
• If ß0 0 the method is explicit, otherwise it is
  implicit.
• The methods has been developed to include a
  predictor step and a corrector step ?
  predictor-corrector method.

96
Multistep, Multivalue, and Predictor-Corrector
Methods 3

• Adams-Bashforth-Moulton method
• Adams-Bashforth predictor
• Adams-Moulton corrector

97
Boundary value problems

• Boundary value ODE.
• The solution for a set of n ODEs is sought after
• where vectors y and f contain n components. In
  addition the solution is bound by two distinct
  points (a and b) of the interval of interest
• Here the boundary condition is given by

98
Two-point Boundary Value Problems

• The boundary value ODE problem is basically
  solved by two techniques
• Shooting method
• Relaxation method (also called difference method)
• As the solution must satisfy two points, a point
  in the interval beginning and also at the end,
  this is not so straightforward as the initial
  value problem.
• The boundary value problem solution requires
  iteration.

99
The Shooting Method

• In the shooting method, using the following
  example problem,
• the equation is transformed in to form
• where the initial value y(a)a is fixed.

100
The Shooting Method 2

• The solution of the latter equation at point b is
  a function of the shooting angle ?, i.e. the
  y(a), and we mark it with g(?).
• The equation to be solved is now
• In practice the shooting angle is first guessed
  and problem II is solved. This yields at point
  tb the approximation g(?0).
• If it is close to y(b)ß we have solution,
  otherwise choose a new shooting angle and
  continue. After two approximations secant method
  can be used.

101
Relaxation Methods

• In the relaxation method the interval a,b is
  divided into m equal parts, so that h(b-a)/m,
  tiaih, i0,,m.
• Denote the approximation of the accurate solution
  of y(ti) by yi.
• Then replace the derivatives with approximations

102
Relaxation Methods 2

• This leads into a new difference equation set
• These are easily expressed and manipulated in
  matrix form.
• For different types of boundary value problems
  special methods have been developed.
• Difference methods are used extensively in
  solving sets of partial differential equations,
  where boundary problems are typical.

103
Selection of methods and some recommended
computer programs

• The number of routines for solving ODEs is big.
  There are many routines in IMSL- and
  NAG-packages.
• The best programs can analyze whether the problem
  is stiff and based on the analysis they are able
  to change the solution method, step size or order
  of the method.
• Some programs are also able to recognize slight
  discontinuities and are able to manage with these.

104
Selection of methods and some recommended
computer programs

• Stiff problems
• Only so called backward type programs are
  suitable for stiff problems.
• There are no available program packages of
  implicit Runge-Kutta methods.
• The recommended programs are LSODA, LSODAR. They
  are also applicable to non-stiff problems and can
  find discontinuities.
• DSTIFF program also work with this class of
  problems.

105
Selection of methods and some recommended
computer programs

• Non-stiff problems, the determination of the
  derivative f(x,y) is cheap
• Runge-Kutta works well. There are different
  orders of Runge-Kutta solvers. Some programs
  RKSUITE, DEGRK, XRK.
• In Matlab there are also Runge-Kutta solvers
  ode25 and ode45.

106
Selection of methods and some recommended
computer programs

• Non-stiff problems, the determination of the
  derivative f(x,y) is expensive
• All predictor-corrector methods are useful.
• E.g. DVDQ, STEP and VOAS
• Adams-Moulton technique based programs are
  suitable for this class of problems, examples of
  programs DIFSUB, GEAR, EPISODE and LSODA.

107
Selection of methods and some recommended
computer programs

• High accuracy is required
• The extrapolation methods are excellent provided
  that the solution curves are smooth enough (no
  discontinuities in derivatives).
• Extrapolation methods can be used also for stiff
  problems but then it is expensive as the number
  of derivative determinations is high.
• Extrapolation methods are not suitable if the
  solution is to be found on many points very close
  to each other.
• A recommended program is METAN1.

108
Selection of methods and some recommended
computer programs

• Higher order problems
• The Taylor series methods are suitable, because
  the there is need then to transform the system in
  to a set of first order ODEs.
• Taylor series methods generate the higher
  derivatives utilizing recursive algorithms.
• One on the best known programs in this class in
  ATSMCC.

109
Selection of methods and some recommended
computer programs

• Most of the program codes discussed earlier are
  available in netlib
• http//www.netlib.org/
• Many mathematical software package has its own
  selection of ODE solvers, but rarely these are
  suitable for all cases.