Title: Ordinary differential equations
1Ordinary differential equations
- KE-42.4520
- Process modeling methods and tools
- Kari I. Keskinen
- 3.3.2008
2Introduction
- 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.
3Plug flow reactor
4Plug flow reactor
- Let us inspect a volume ?V, where the reaction
rate can be assumed constant.
y?y
y
Fj(y?y)
?V
5Plug flow reactor
- Reaction rate is constant in ?V
Steady state assumption, i.e. no accumulation
6Plug flow reactor
- For the volume ?V we can write the material
balance for component j.
A cross-sectional area
7Plug flow reactor
8Plug flow reactor
Integral form of the dimensioning equation.
9Classification 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).
10ODEs 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.
11ODEs 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
12ODEs 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.
13ODEs 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
14ODEs 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.
15ODEs 6
- Linear ODE
- The coefficient a do not depend on the variable
y. f(x) is a forcing function. - Nonlinear ODE
16ODEs 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 )
17ODEs 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
18ODEs 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
19ODEs 10
- can be transformed into a set of n coupled first
order differential equations
This is the RHS referred later.
20Solution 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
21Solution 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.
22Series 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.
23Series 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.
24Series 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.
25Laplace 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.
26Laplace 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.
27Laplace 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).
28Huokosdiffuusio
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
29Huokosdiffuusio
?r ? 0
30Huokosdiffuusio
- reaktio A ? B
- ? ekvimolaarinen vastavirtaus (lisäksi vakio T
ja p) - aiemmin todettu
- Yhtälö soveltuu myös laimeille seoksille.
- Sijoittamalla
31Huokosdiffuusio
32Huokosdiffuusio
Differentiaaliyhtälö diffuusio reaktio
33Huokosdiffuusio
- Reunaehdot
- CA äärellinen, kun r 0
- CA CAs, kun r R
34Huokosdiffuusio
- ratkaistaan dimensiottomien muuttujien avulla
- uudet reunaehdot
- CA CAs r R
CA äärellinen r 0 ?
äärellinen ? 0
35Huokosdiffuusio
- Differentiaaliyhtälö dimensiottomien muuttujien
avulla
36Huokosdiffuusio
37Huokosdiffuusio
38Huokosdiffuusio
- 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
39Huokosdiffuusio
? 1 ? 1
Reunaehdot
? äärellinen ? 0
Ratkaistaan uuden muuttujan avulla
y ? ? ? (1/?) y
40Huokosdiffuusio reaktio
Sijoittamalla
41Huokosdiffuusio reaktio
42Hyperboliset funktiot
43Huokosdiffuusio reaktio
- Reunaehdot
- ? 0 ? äärellinen
- cosh?1? ?1
sinh?1? ? 0
1/? ? ? ? äärellinen
A1 0
44Huokosdiffuusio reaktio
? sijoittamalla A1 ja B1
45Huokosdiffuusio reaktio
- Pienet ?1n arvot ? pintareaktio kontrolloi
eli An konsentraatio on suuri huokosissa (ennen
kuin A reagoi)
46Non-isothermal porous catalyst particle
- For a non-isothermal porous catalyst particle
the following equation has been derived.
47Non-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
48Non-isothermal porous catalyst particle
- The differential equations can be combined ?
equation between temperature and concentration
kt effective heat transfer coefficient
49Non-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
50Non-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!
51Non-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.
52Non-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 ? ?
53Internal effectiveness as a function of Thiele
modulus
54Example 1
- Classify
- Linear Nonlinear
- Homogeneous Inhomogeneous
- Order
- Degree
- Forcing function
55Example 1
- Classify
- Linear Nonlinear (if nltgt1)
- Homogeneous Inhomogeneous
- Order 2
- Degree 1
- Forcing function? 0
56Example 2
- Classify
- Linear Nonlinear
- Homogeneous Inhomogeneous
- Order
- Degree
- Forcing function?
57Example 2
- Classify
- Linear Nonlinear
- Homogeneous Inhomogeneous
- Order 1
- Degree 2
- Forcing function? x - 1
58Example 3
- Classify
- Linear Nonlinear
- Homogeneous Inhomogeneous
- Order
- Degree
- Forcing function?
59Example 3
- Classify
- This is actually a differential-algebraic
equation as there is also an algebraic equation
in addition to a differential equation.
60Why 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.
61Stability 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.
62Stability 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.
63Stability 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.
64Solutions by Euler method
Model
Initial condition
Closure first order reaction
65Separable form, analytical solution possible
66Numerical solution Eulers method
Explicit Euler if
is calculated at tn
Implicit Euler if
is calculated at tn1
67Numerical solution Eulers method
For the test problem
Explicit Euler
Implicit Euler
68Numerical 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
69Example
70Example
71Example
72Example
73Example
74Example
75Shorter time steps
76Other options
at average of tn and tn1
77Crank-Nicolson method
78Crank-Nicolson method
Stable, but oscillates!
79Solving 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.
80Solving 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.
81Solving 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.
82Runge-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.
83Runge-Kutta 2
- This is the second-order Runge-Kutta method. Can
be extended to fourth-order Runge-Kutta method as
follows
84Runge-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.
85Richardson 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.
86Bulirsch-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).
87Bulirsch-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.
88Stiff 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.
89Stiff 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.
90Stiff sets of equations 3
- Consider
- With boundary conditions
- A solution is found
91Stiff 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 .
92Stiff 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.
93Stiff 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.
94Multistep, 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.
95Multistep, 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.
96Multistep, Multivalue, and Predictor-Corrector
Methods 3
- Adams-Bashforth-Moulton method
- Adams-Bashforth predictor
- Adams-Moulton corrector
97Boundary 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
98Two-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.
99The 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.
100The 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.
101Relaxation 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
102Relaxation 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.
103Selection 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.
104Selection 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.
105Selection 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.
106Selection 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.
107Selection 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.
108Selection 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.
109Selection 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.