Section 4: Dynamic Systems and Differential Algebraic Equations

Section 4.1: Ordinary Differential Equations

Section 4.1Ordinary Differential Equations

In the dynamic modelling of chemical engineering systems, ordinary differential equations (o.d.e.s) occur almost exclusively as unsteady state material balances. We will therefore deal only with the kind of o.d.e.s that occur in this context.

Single O.D.E.s

The simplest single equation of the appropriate form is:

dy(t)/dt = -y

In this equation t, representing time, is called the independent variable. The single unknown or dependent variable is y. The solution is strictly speaking y(t), i.e. y as a continuous function of t. Where the equation is solved numerically this will be a set of discrete values of y at corresponding values of t.

To enable o.d.e.s to be solved it is necessary to have a known fixed point in the solution space. For the present type of problem this will invariably be a starting or initial value of y at an initial time normally designated as t=0, i.e.:

y = y₀ at t=0

This equation is one of the few o.d.e.s that has an analytical solution. If we solve it with initial condition:

y = 1 at t=0

We obtain:

y(t) = exp(-t)

To obtain a numerical solution we write:
dy/dt ~ dy/dt

Here dy and dt represent small changes in y and t rather than differentials. Hence we can calculate the change in y for a given change in t or timestep dt approximately as in general:

dy = dy/dt dt

Hence we can calculate a `new' y value since:
y(t + dt ) = y (t) + dy = y (t) + dy(t)/dt dt

For this particular equation therefore:
dy = - y dt

And so:
y(t + dt ) = y(t) - dt y(t)

Starting with a known value of y = y₀ = 1 at t = 0:
y₁ = y(dt) = - y₀ dt + y₀
y₂ = y(2dt) = - y₁ dt + y₁
y₃ = y(3dt) = - y₂ dt + y₂
... etc.

This numerical solution is implemented in a simple spreadsheet. Try changing the size of the timestep in cell B7 between say 0.01 and 0.5 and see how the numerical solution compares to the analytical one.

Timeconstants, Timesteps and Stability

The equation used as an example above is a special case, with k = 1:

dy(t)/dt = - k y

Time Constants and Time Steps

Note that k has dimensions of time^-1 and so we can alternatively write:

dy(t)/dt = - y / T

T is called the time constant of the equation. The analytical solution of this equation is now:

y(t) = exp(-t / T )

Here is a spreadsheet for the modified problem. You can change the time constant. If you experiment with different time constants and time steps, you will observe that the accuracy of the solution depends on the ratio of the time step. Be careful at this stage to keep the time step always less than the timeconstant.

In discussion of lags in `black box' models it was shown that the time constant of a linear o.d.e. is a measure of how fast its solution progresses towards a final steady state. E.g. 63% in a time equal to the time constant and 90% in approximate 2.3 times the time constant. If we wanted to observe 99% of the response towards steady state this would require a time given by:
0.01 = exp (-t/T)

This gives gives a time approximately 4.6 times T.

Stability

Now make the time constant and time step exactly the same, e.g. 1. Note that the solution goes immediately to zero after the first timestep and remains there. The solution behaves in general as it should, i.e. it starts at 1.0 and reduces monotonically to zero, the correct final steady state vale, and stays there.

Now make the time step just larger than the time constant, e.g. 1.1. Observe that the solution now starts to behave in a quite different, and completely incorrect oscillatory manner. If the time step is made exactly twice the time constant the solution oscillates between a maximum and minimum value. If the timestep is made any larger, then the oscillations will grow and eventually become infinite.

Sets of O.D.E.s

It is possible to solve several equations which make up a set of o.d.e.s using the above method applied to each of them. For example to solve:

dx/dt = - x
dy/dt = x - y

Starting from initial values y₀ and x₀ we just apply repeatedly the two relationships:
x₁ = - x₀ dt + x₀
y₁ = (x₀ - y₀) dt + y₀

Here is a spreadsheet which solves a slightly more general form of these equations where each has a separately adjustable time constant, denoted by T1 and T2 respectively in the spreadsheet program.

`Stiffness'

With sets of equations the choice of timestep becomes a significant issue.

It is clear from our investigation of stability that we must not use a time step greater than the time constant, and for reasonable accuracy it should normally be less than one tenth of the timestep. Since we would usually wish to follow most of the behaviour of the system towards its final steady state, since 4.6 times the time constant takes us 90% of the way there for a single equation a total simulation time of perhaps 5 times the time constant should be adequate,

For a single equation, or a pair of equations like those above which have the same time constant, a total of around 50 time steps would be adequate to provide sufficient accuracy and to cover the range of `interesting' behaviour.

Unfortunately, for a set of equations, the shortest time constant governs the acceptable choice of time step and longest determines the length of time for which the simulation must be run. Experiment with this using the program supplied. With the second time constant of say 5 the system will still be a long way from its final steady state, which is both variable and equal to zero, after a simulation time of 2.5. However, making the first time constant too small, i.e. less than 0.1 with the time step of 0.1, still causes instability.

This phenomenon is call stiffness and becomes a problem whenever time constants differ by more than an order of magnitude. Since in practice these may differ by several factors of 10, more sophisticated numerical methods than the one described above often have to be used to allow larger timesteps to be taken.

Section 4.2: The Systematic Construction of Dynamic Models

Section 4.2.1: Introduction

A complete discussion of how to develop chemical engineering models is really an complete course in chemical engineering, with practical experience besides! The objective of this section is thus only to show that there is a good, as opposed to a bad, way of constructing models. If this approach, which is systematic, is followed, then, assuming that the model builder's chemistry, physics and chemical engineering understanding of the process being modelled is correct, then there is a better chance than otherwise that the final model will be correct.

The basis of this approach to model building is that the equations which constitute a model are not arbitrary mathematical entities, but have a consistent physical basis. There are certain types of equation which describe different aspects of a model. A knowledge of this helps to ensure that all equations are written down.

Another point is that the set of model equations should be written without regard to how they might be solved. A confusing feature of some of the classical textbooks in chemical engineering is that the author, who knows the answer he wants to reach, slips steps in the solution process in amongst the construction of the model. The way to avoid confusion is to write all the equations down without attempting to perform any algebraic manipulation or simplification.

This section might therefore be subtitled: How to avoid algebra. The final models may be less elegant than those in some text books, but the computer will have no difficulty in solving them.

Section 4.3: Differential Algebraic Equation Systems

Section 4.3: Differential Algebraic Equations Systems

As should be clear from the discussion of systematic model construction in section 4.2, all `real' models involve a mixture of differential and algebraic equation, which is formally described as a system of differential-algebraic equations or d.a.e.s.

Introduction

D.a.e.s have properties which are a mixture of those of o.d.e.s and algebraic equations, and in general may require a combination of solution techniques. We will deal here only with the simplest case of d.a.e.s, for which solution is straightforward provided certain conditions are met. There are possible complications which we will not discuss here.

If the following conditions are met, then solution of a set of d.a.e.s can be handled simply as described below.

The system contains only algebraic equations and o.d.e.s with a single independent variable (normally time).
There are as many equations in total, either differential or algebraic, as there are variables excluding the independent variable.
Initial values are given for as many variables as there are o.d.e.s.
All variables which do not appear as derivatives occur at least once in the algebraic equations.
The algebraic equation subset may be solved using the normal methods available for a.e.s. if values are known for all variables which appear as derivatives.

We will also assume the following additional restrictions which are not fundamental but which make solution very straightforward:

initial values are given for all variables which appear as derivatives, and
the algebraic subset is amenable to non-iterative `direct' solution as described in section 3.1.

Algorithm for D.A.E. Solution

Consider first the following arbitrary set of 4 d.a.e.s:

dx/dt = w - y
dy/dt = y - x
z = y - x
w = z - xº

With initial conditions:
x=1 and
y=0 at t=0

Note that:

There are only a.e.s and o.d.e.s with one independent variable t.
There are 4 equations and 4 variables w, x, y, z
There are 2 o.d.e.s and two initial conditions are given.

This satisfies conditions 1-3 given above. Furthermore:

In the 2 o.d.e.s x and y appear in derivatives. These will be referred to as the differential variables.
The remaining variables w and z, which we will call the algebraic variables both appear in the two a.e.s.
Thus if values were known for the differential variables, the algebraic subset would consist of two equations in two unknowns and would, in principle, be solvable.

We have thus satisfied all the main conditions. Additionally in this example:

Initial conditions are given for both the differential variables.
The a.e.s are already both written as formulas which enable the two algebraic variables to be calculated directly.
The solution algorithm involves three steps, one of which is performed only once, at the start of the solution procedure. The other two must each be carried out at each timestep.
Initialisation (once only):

Initial conditions being given for the differential variable, the a.e.s can be solved for the algebraic variables.
The values of all variables are now known at time t=0.
The next two steps must be performed at each timestep. The procedure starts with values for all variables known at a general time t and finishes with them calculated at the `new' time (t+dt).
O.D.E solution:

All variables are known at time t (time t = 0 on the first step)
The right hand sides of the o.d.e.s, from which values of derivatives can be calculated, are evaluated and the derivatives determined.
The ``new' values of the differential variables are calculated by a suitable o.d.e. solving procedure, e.g. using:
y(t + dt ) = y (t) + dy = y (t) + dy(t)/dt dt
After this has been performed for all o.d.e.s the values of the differential variables are known at time (t+dt).
Algebraic equation solution:

The differential variable values being calculated above, the a.e.s. now involve only the algebraic variables as unknowns.
Solve the a.e.s
All variables are now known at time (t+dt).
Notice that the final step is actually the same as the initialisation, namely solving a set of algebraic equations. If we are given initial values for the differential variables then it in fact involves solving the identical set of equations for the identical variables. In many practical problems the steps differ in detail because the initial conditions are not given as values of the differential variables. This is quite acceptable so long as the correct number of conditions is given. The Repeated a.e. solution step however always involves solving for the algebraic variables. In many practical problems an iterative solution will be required.
Example
A simple system of two o.d.e.s and one a.e. is given under the description of our modelling language. (section L2.3).
This illustrates both the nature of a d.a.e. system and the approach to its solution.
Section 4.4: Additional Theory and Notes
Section 4.4: Additional Theory and Notes

This section has a more mathematical discussion of ordinary differential equations and is provided to complement the other notes in the section.

4.4.1 Steady State Model

4.4.2 The Richness of Dynamics

4.4.2.1 Analytic Solution of Gravity Outflow Problem with Constant Feed

4.4.3 Numerical Solution of an ODE

4.4.4 Time constant, Accuracy and Step Size

4.4.5 Numerical Instability

4.4.6 Systems of ODE's

4.4.6.1 Choosing a Step Size for More than One ODE

4.4.7 Other Methods: Euler is not the Only Method

4.4.7.1 Demonstration that Improved and Modified Euler are Second Order Methods

4.4.8 Applications of ODE's

4.4.8.1 Implicit Differentiation of Redlich-Kwong equation

4.4.9 DAE systems

4.4.9.1 DAE Formulation of Flow Controller problem

4.4.10 Towards Process Simulation

4.4.10.1 Material and energy balance for a tank

Section 4.4.1: Steady State Model
Section 4.4.1: Steady State Model

At any place in the process,nothing changes with time.
Described by Algebraic Equations (AEs).

Material balance:
$\sum$ flows in = $\sum$ flows out

Dynamic Model

Quantities may change with time.

Described by Ordinary Differential Equations (ODEs) and AEs.

Material balance:
rate of change of mass
= $\sum$ flows in - $\sum$ flows out
In a molar or volumetric balance, mass is replaced by moles or volume.

Flow through a tank

Steady state material balance

F₁ - F₂ = 0

Dynamic material balance

dM/dt = F₁ - F₂

Dynamic balance:

Introduces holdup M
ODE not AE

More equations needed

Dynamic material balance
dM/dt = F₁ - F₂

Need two constraint equations e.g. for flows

F₁ - a₁ = 0

F₂ - a₂ = 0

Now we have a DAE system. Get back to ODE by substituting the AEs into the material balance.
Some DAE systems are hard! Care must be used in formulating DAE systems. (See later:sections 4.4.9 and section 4.4.10).
Solve ODE for holdup
Substitute flow constraints $\rightarrow$ dynamic material balance

dM/dt = a₁ - a₂

Solution:

M(t) = (a₁ - a₂)t + c

or

M(t) = bt + c

Initial Conditions
ODEs need initial conditions

one for each differential variable
valid at initial time (t₀)

Here : one differential variable, M
Give initial holdup M₀ at time t₀ = 0
Substitute into ODE:

$\begin{displaymath}M_{0} = b \times 0 + c \end{displaymath}$

gives

c = M₀

Hence particular solution

M(t) = M₀ + bt

Section 4.4.2: The Richness of Dynamics
Section 4.4.2: The Richness of Dynamics

A steady state model is described by AEs, a dynamic model by ODE(s) or DAEs.
A dynamic model requires initial condition(s) (ICs) for the differential variable(s)
More constraints are usually required in a dynamic model than in a steady state model of equivalent scope.In the tank model
steady state only one of the flows (in or out) may be constrained.
dynamic model both flows are constrained.

A dynamic model usually contains more detail (more variables) than a steady state model of equivalent scope. In the tank model

introduce the holdup in the dynamic model.
may have to introduce an extra stream to provide a correct mass balance if the tank overflows.

Avoid infeasible specifications

If we specify a constant outflow greater than the inflow (also constant) the tank may empty.
A careless simulation could then lead to the holdup becoming negative.
Care must be taken to arrange process specifications which affect dynamics so that they do not lead to infeasible situations.

A more difficult example

Tank Outlet flow proportional to $\sqrt M$
.. from Bernoulli, if

liquid has constant density
tank has constant cross sectional area

Dynamic material balance

$\begin{displaymath}dM/dt = a_{1} - k \sqrt M \end{displaymath}$

Can we solve it analytically?

$\begin{displaymath}dM/dt = a_{1} - k \sqrt M \end{displaymath}$

$\begin{displaymath}dM / ( a_{1} - k \sqrt M ) = dt \end{displaymath}$

$\begin{displaymath}\int_{M_{0}} dM / ( a_{1} - k \sqrt M ) = \int_{0} dt = t \end{displaymath}$

Needs a change of variable.
Details in section 4.4.2.1 if interested.
Result:

$\begin{displaymath}k^{2}t/2 = k (\sqrt M_{0} - \sqrt M) - a_{1} ln (\frac{a_{1} - k \sqrt M}{a_{1} - k \sqrt M_{0}} ) \end{displaymath}$

We get t as function of M, but we want M as function of t: Awkward!

Summary

Dynamic models of process systems $\rightarrow$ ODEs or DAEs
An initial condition is needed for each differential variable
Some simple process models can be solved exactly
Most need numerical solution

Section 4.4.2.1: Analytic Solution of Gravity Outflow Problem with Constant Feed
Section 4.4.2.1: Analytic solution of gravity outflow problem with constant feed

The dynamic material balance is

$\begin{displaymath}dM/dt = a_{1} - k \sqrt M \end{displaymath}$

Separation of variables gives

$\begin{displaymath}dM / ( a_{1} - k \sqrt M ) = dt \end{displaymath}$

Then integrate from M = M₀ and t = 0 to get

$\begin{displaymath}\int_{M_{0}} dM / ( a_{1} - k \sqrt M ) = \int_{0} dt = t \end{displaymath}$

Try a change of variable of integration to u, where

$\begin{displaymath}u = a_{1} - k \sqrt M \end{displaymath}$

and

$\begin{displaymath}du = - (k/2 \sqrt M) dM \end{displaymath}$

or

$\begin{displaymath}dM = - (2 \sqrt M / k) du = - (2 (a_{1} - u) / k^{2}) du = 2 ((u - a_{1}) / k^{2}) du \end{displaymath}$

With this substitution the integral over M is transformed to

$\begin{displaymath}\int_{u_{0}} du (2 (u - a_{1})) / (k^{2} u) = (2 / k^{2}) \int_{u_{0}} du ( 1 - a_{1} / u) \end{displaymath}$

where $u_{0} = a_{1} - k \sqrt M_{0}$ .
The integral is

(2/k²) ( u - u₀ - a₁ ( ln u - ln u₀ )

or

(2/k²) ( u - u₀ - a₁ ln (u / u₀ )

This integral is equal to t - 0, or just t. Thus with a little rearrangement we get

k²t/2 = u - u₀ - a₁ ln (u/u₀)

or expressed in terms of M again,

$\begin{displaymath}k^{2}t/2 = k (\sqrt M_{0} - \sqrt M) - a_{1} ln (\frac{a_{1} - k \sqrt M}{a_{1} - k \sqrt M_{0}} ) \end{displaymath}$

Section 4.4.3: Numerical Solution of an ODE
Section 4.4.3: Numerical Solution of an ODE

Finite differences
The key to any numerical method for ODEs is to make a finite difference approximation to the derivative.

m₀ = slope of tangent at t₀
m₁ = slope of tangent at t₁
m = slope of secant from t₀ to t₁

$\begin{eqnarray*}m_{0} = \frac{dx}{dt} \vert _{t=t_{0}} \\ m_{1} = \frac{dx}{dt} \vert _{t=t_{1}} \\ m = \frac {x_{1} - x_{0}}{t_{1} - t_{0}} \end{eqnarray*}$

m is well approximated by m₀ or m₁ if

curvature of x(t) slight over the time step t₀ to t₁;
no discontinuities or sharp kinks in x(t).

Forward difference

$\begin{eqnarray*}\frac{dx}{dt} \vert _{t=t_{0}} \approx \frac {x_{1} - x_{0}}{t_{1} - t_{0}} \\ x_{1} \approx x_{0} + (t_{1} - t_{0}) m_{0} \end{eqnarray*}$

Backward difference

$\begin{eqnarray*}\frac{dx}{dt} \vert _{t=t_{1}} \approx \frac {x_{1} - x_{0}}{t_{1} - t_{0}} \\ x_{1} \approx x_{0} + (t_{1} - t_{0}) m_{1} \end{eqnarray*}$

Explicit Euler method
Forward difference $\rightarrow$ Explicit Euler method
(Backward difference $\rightarrow$ Implicit Euler - see later)
Solve

$\begin{displaymath}\frac {dx}{dt} = f(x) \end{displaymath}$

Specify x₀ at initial time t₀.
Calculate f(x) at t₀ $\rightarrow$ f₀
Finite difference formula:

x₁ = x₀ + h f₀

Algorithm: Euler's method

! The ODE is in general $\frac{dx}{dt} = f(x,t)$
! The right hand side function may
! depend on time t explicitly
! as well as implicitly through
! the behaviour of x.
Specify initial conditions:

time t₀
value of differential variable x₀
Specify final time t_f and step size h.
Initialise counter n = 0 to count steps.
do while t_n < t_f
! we still have some way to go

set t_n+1 = t_n + h
if t_n+1 > t_f then
! truncate the final step set t_n+1 = t_f
! for the last step only redefine h = t_f - t_n
end if
evaluate f_n = f(x_n, t_n)
calculate x_n+1 = x_n + h f_n
increment n by 1
end do
return

Numerical solution of a tank flow model
ODE: $\frac{dM}{dt} = b$ with b = 1 kg/s CORRECTION!
IC: M₀ = 500 kg at t₀ = 0 s
Final time: t_f = 100 s
Step size: h = 20 s

Solution

step number time t (s) holdup M (kg)

0 0 500

1 20 520

2 40 540

3 60 560

4 80 580

5 100 600

Numerical = Exact ?

The numerical solution equals the exact solution because the exact solution is linear.
This statement must be qualified because computers do not work in exact arithmetic. Real numbers are stored in a finite form.
Perhaps 500.00000 + 20.000000 is calculated as 519.999999.
Rounding error accumulates in computer arithmetic because of the finite representation of numbers.
Termination logic

A further consequence of rounding error.
The algorithm contains tests to decide whether

to truncate the last step,
to do another step.

The latter may have unexpected effects if, say,
80 + 20 = 99.99999

The algorithm will execute the loop again after reaching `100' (or as the computer thinks, 99.99999).
It will chop h to 0.00001.
Probably not serious, but output may look odd.

Another example: linear ODE

$\begin{displaymath}\frac{dx}{dt} = - k x \end{displaymath}$

This does not mean that the solution is linear. Linear means that x appears in the equation to first degree only. (No x², x³.) Many applications in chemical engineering.
Exact solution is known:

x(t) = x₀ e^-kt

k is a rate constant for decay of x.
Units: (1/seconds).

Numerical solution

Parameter values:

k = 1 s^-1
IC: x₀ = 100
at: t₀ = 0 s
final time: t_f = 2 s
step size: h = 0.25 s

Result

step time (s) x(Euler) x(Exact) abs error

0 0.00 100.0000 100.0000 0.0000

1 0.25 75.0000 77.8801 2.8801

2 0.50 56.2500 60.6531 4.4031

3 0.75 42.1875 47.2367 5.0492

4 1.00 31.6406 36.7879 5.1473

5 1.25 23.7305 28.6505 4.9200

6 1.50 17.7979 22.3130 4.5151

7 1.75 13.3484 17.3774 4.0290

8 2.00 10.0113 13.5335 3.5222

Errors

Absolute error

e_abs,n = | x_n,Euler - x_n,exact |

Relative error

e_rel,n = | (x_n,Euler - x_n,exact) / x_n,exact |

step time (s) abs error rel error

0 0.00 0.0000 0.0000

1 0.25 2.8801 0.0370

2 0.50 4.4031 0.0726

3 0.75 5.0492 0.1069

4 1.00 5.1473 0.1399

5 1.25 4.9200 0.1717

6 1.50 4.5151 0.2024

7 1.75 4.0290 0.2319

8 2.00 3.5222 0.2603

Summary

Finite difference approximation $\rightarrow$ Euler's method for ODEs
Easy to apply explicit Euler: If $\frac{dx}{dt} = f(x,t)$ , calculate f(x,t) at start of each step where x and t already known
If exact solution is linear, no errors
Linear ODE, exponential solution, $\rightarrow$ numerical solution has errors.

Section 4.4.4: Time constant, Accuracy and Step Size
Section 4.4.4: Time constant, Accuracy and Step Size

Error in numerical solution related to curvature in true solution:
second derivative $\frac{d^2x}{dt^2}$

Solve $\frac{dx}{dt} = - k x$
One step, size h.
Euler: x₁ = x₀ (1 - hk)
Exact: $x_1 = x_0 \exp (-hk)$

Choice of time step

k : rate constant (s^-1)
1/k: time constant
or characteristic time (s)
h : time step (s)
$\epsilon$ : desired local truncation error (units of x)

$\begin{displaymath}x_0 (hk)^2 /2 = \epsilon \end{displaymath}$

or

$\begin{displaymath}h = \sqrt{(2 \epsilon / x_{0})} / k \end{displaymath}$

General non linear equation: solution x(t)

$\begin{displaymath}(h^{2} / 2) d^{2}x/dt^{2} = \epsilon \end{displaymath}$

or

$\begin{displaymath}h = \sqrt{2 \epsilon/(d^{2}x/dt^{2})} \end{displaymath}$

Accumulation of error
k=1, x₀ = 100
One step: h = 0.04s
Exact: $x(0.04) = 100 \exp(-0.04) \\ = 96.079$
Euler: $x(0.04) = 100 - 0.04 \times 100 \\ = 96.000$
Error = 0.079

Two steps: h = 0.02s
First step: $x(0.02) = 100 - 0.02 \times 100 \\ = 98.000$ (exact: 98.020)
Second step: $x(0.04) \\ = 98.000 - 0.02 \times 98.000 = 96.040$
Error = 0.039

Global error

Local truncation error proportional to h².Gives error over one step. Subsequent steps start with an error. Global error proportional to h.
Decaying exponential (linear ODE):
Local error declines as steady state is approached.
Global error also declines.
Maximum global error at $t \approx 1/k$ .

Material balance for dilute species flowing through mixing tank:

V dc/dt = q c₁ - q c

or

T dc/dt = c₁ - c

where

T = V/q

Example: V = 0.5 m³, q = 0.002 m³/s
so T = 250 s
c₁ = c_f = 0.1 kmol/m³
c₀ = 0 at t₀ = 0
Exact solution: $c = 0.1 ( 1 - \exp (-0.004 t))$

Choose step:
for $\epsilon = 0.001$ kmol/m³ (1% of st.st)
$0.1 \times (0.004h)^{2} / 2 = 0.001$
h = 35.355

Results

step time c Euler c exact abs error

0 0.000 0.000000 0.000000 0.000000

1 35.355 0.014142 0.013188 0.000954

2 70.710 0.026284 0.024636 0.001648

3 106.065 0.036709 0.034575 0.002134

4 141.420 0.045660 0.043203 0.002457

5 176.775 0.053344 0.050693 0.002652

6 212.130 0.059942 0.057195 0.002747

7 247.485 0.065607 0.062840 0.002767

8 282.840 0.070471 0.067741 0.002731

9 318.195 0.074647 0.071995 0.002652

Summary

For desired local accuracy we can choose a step size which depends on time constant 1/k for linear system, or alternatively some estimate of second derivative.
Errors accumulate
Global error > local error after one step.
Euler's method:
local error $\propto h^2$
global error $\propto h$
linear ODE - decaying exponential
max global error at $t \approx 1/k$

Variable step size
Approximate d²c/dt² from previous three points.

$\begin{displaymath}\frac{d^2c}{dt^2} \approx \frac{2}{h_n + h_{n-1}} [\frac{c_n - c_{n-1}}{h_n} - \frac{c_{n-1} - c_{n-2}}{h_{n-1}}] \end{displaymath}$

where

h_n = t_n - t_n-1

h_n-1 = t_n-1 - t_n-2

We then estimate a new step size using

$\begin{displaymath}h = \sqrt{2 \epsilon/(d^{2}c/dt^{2})} \end{displaymath}$

step time c step size abs error

0 0.000 0.000000

1 35.355 0.014142 35.355000 0.000954

2 70.710 0.026284 35.355000 0.001648

3 106.065 0.036709 35.355339 0.002134

4 144.222 0.046369 38.156248 0.002533

5 186.208 0.055376 41.986497 0.002857

6 231.751 0.063505 45.542669 0.003079

7 281.040 0.070700 49.289334 0.003193

8 335.044 0.077030 54.003740 0.003210

9 394.952 0.082534 59.907804 0.003135

10 462.030 0.087220 67.078139 0.002974

11 537.973 0.091102 75.943619 0.002729

12 625.322 0.094211 87.348124 0.002409

13 727.867 0.096586 102.545520 0.002025

14 851.443 0.098273 123.575870 0.001592

15 1000.000 0.099299 148.557098 0.001131

Section 4.4.5: Numerical Instability
Section 4.4.5: Numerical Instability

High accuracy needs small steps: $hk \ll 1$ But if accuracy is not important, can we take large steps?
Example:

$\begin{displaymath}\frac{dc}{dt} = 0.004 (0.1 - c) \end{displaymath}$

with c₀ = 0.099 kmol/m³ at t = 0 (nearly at steady state).
Try h = 900s = 15 min.

step time (s) c

0 0.0 0.099000

1 900.0 0.102600

2 1800.0 0.093240

3 2700.0 0.117576

4 3600.0 0.054302

5 4500.0 0.218814

6 5400.0 -0.208916

7 6300.0 0.903181

8 7200.0 -1.988271

Numerical instability if time step is large.

Reason for the instability:
1 - hk negative if $hk > 1 \Longrightarrow$ oscillations;
1 - hk < -1 if $hk > 2 \Longrightarrow$ oscillations grow.
Implicit Euler method
Avoid extrapolation of derivative from start of step.
Approximate dx/dt by backward difference.

$\begin{displaymath}\frac{dx}{dt} \vert _{t=t_{1}} \approx \frac {x_{1} - x_{0}}{t_{1} - t_{0}} \end{displaymath}$

So

$\begin{displaymath}\frac{x_1 - x_0}{h} \approx f(x_1) \end{displaymath}$

Implicit Euler (for first step):

x₁ = x₀ + h f(x₁)

Must be treated as an equation for x₁. If f(x) is linear, solution is easy. If f(x) is non-linear, we have a NLAE; solution needs iteration.

$\begin{displaymath}\frac{dx}{dt} = - k x \end{displaymath}$

Implicit Euler:

x₁ = x₀ - h k x₁

So

$\begin{displaymath}x_1 = \frac{x_0}{1 + hk} \end{displaymath}$

Does this overcome the problem of instability?

Implicit Euler:

$\begin{displaymath}\frac{x_1}{x_0} = \frac{1}{1 + hk} \end{displaymath}$

Always between one and zero for positive k.
Tends to zero at large k.
`Strongly A-stable' method.
Example: large steps
Rerun

$\begin{displaymath}\frac{dc}{dt} = 0.004 (0.1 - c) \end{displaymath}$

with c₀ = 0.099 kmol/m³ at t = 0 and h = 900s.

step time (s) Implicit Euler exact error

0 0.0 0.099000 0.099000 0.000000

1 900.0 0.099783 0.099973 0.000190

2 1800.0 0.099953 0.099999 0.000047

3 2700.0 0.099990 0.100000 0.000010

4 3600.0 0.099998 0.100000 0.000002

5 4500.0 0.100000 0.100000 0.000000

6 5400.0 0.100000 0.100000 0.000000

7 6300.0 0.100000 0.100000 0.000000

8 7200.0 0.100000 0.100000 0.000000

Non-linear equation: solve a NLAE at each time step.
For

$\begin{displaymath}\frac{dx}{dt} = f(x,t) \end{displaymath}$

implicit Euler at time step n+1 is

x_n+1 - x_n = h f(x_n+1,t_n+1)

We need to solve

g(x_n+1) = x_n+1 - x_n - h f(x_n+1,t_n+1) = 0

for x_n+1.
Choice of NLAE method:

Newton
bisection
secant
direct substitution

For short time steps direct substitution may work well enough:

x^(k) = x_n + h f(x^(k-1),t_n+1)

where x^(k) is the kth iterate for solving x_n+1.
NB: x_n and t_n+1 are known and stay fixed while iterating.
For a linear ODE the limit at which direct substitution fails to converge is hk = 1.
Newton's method is recommended. Use the value at x_n as a starting guess for x_n+1, i.e.

x⁽⁰⁾ = x_n

Algorithm is similar to explicit Euler except for the point within the time step loop where we calculate x_n+1.
Algorithm: Implicit Euler method

Specify initial conditions: t₀ and x₀.
Specify final time t_f and
step size h.
Initialise step counter n = 0.
Do while t_n < t_f
set t_n+1 = t_n + h
if t_n+1 > t_f then
set t_n+1 = t_f
redefine h = t_f - t_n
end if
! THIS IS DIFFERENT:
solve
$g(x_{n+1}) = \\ x_{n+1} - x_n - h f(x_{n+1},t_{n+1}) = 0$ for x_n+1.
increment n by 1
end do
return

Newton algorithm for NLAE

Initialise x⁽⁰⁾ = x_n.
Do k = 1,max_iterations
! max_iterations is maximum number of Newton iterations at each time step
call deriv_and_fn
(x^(k-1),t_n+1, f,dfdx)
! you must supply a subroutine
! deriv_and_fn (x,t,f,dfdx) which
! when given x and t will calculate
! and return
! f = f(x,t) and dfdx = $\frac{\partial f}{\partial x}(x,t)$
calculate g = x^(k-1) - x_n - h f
if |g| <nlae_tolerance exit
calculate $x^{(k)} = x^{(k-1)} - \frac{g}{1 - h * \frac{\partial f}{\partial x}}$
end do
set x_n+1 = x^(k-1)
Remarks

nlae_tolerance is the equation-solving tolerance for converging Newton's method. This number should be very small compared with the desired integration accuracy for the dynamic simulation (2 or 3 orders of magnitude smaller).
Usually Newton's method converges quickly in 2 or 3 iterations, because the starting guess is good.
Set max_iterations to be about 10 or 20.
For each specific f(x,t) a different version of the subroutine deriv_and_fn (x,t,f,dfdx) will be needed. Given x and t, return f and dfdx.

Example

$\begin{displaymath}\frac{dM}{dt} = 4 - 0.2 \sqrt{M} \end{displaymath}$

with M₀ = 100 kg at t = 0.
Integrate to t = 500 s with h = 50 s.
NLAE tolerance = 0.01 kg.

step time M iterations

0 0.0 100.000 0

1 50.0 169.719 3

2 100.0 221.043 3

3 150.0 259.846 3

4 200.0 289.654 3

5 250.0 312.794 3

6 300.0 330.890 3

7 350.0 345.117 3

8 400.0 356.346 3

9 450.0 365.235 3

10 500.0 372.280 2

Total number of iterations 29

Summary
With large steps (h k > 2) explicit Euler is unstable.
Implicit Euler uses backward difference and is always stable for problems which tend to a steady state, no matter how large the time step.
If the ODE is non-linear, then with implicit Euler we must solve a NLAE at each time step.
Newton's method will converge quickly in most cases, using x at the start of the step as a guess for the end of the step.

Section 4.4.6: Systems of ODEs
Section 4.4.6: Systems of ODEs

Flow through two tanks $\Longrightarrow$ two ODEs

$\begin{displaymath}\frac{dM_1}{dt} = F_0 - F_1 \end{displaymath}$

$\begin{displaymath}\frac{dM_2}{dt} = F_1 - F_2 \end{displaymath}$

Assume F₀ = a, constant;
$F_1 = k_1 \sqrt{M_1}$ , $F_2 = k_2 \sqrt{M_2}$
Two ODEs in two unknowns M₁, M₂:

$\begin{displaymath}\frac{dM_1}{dt} = a - k_1 \sqrt{M_1} \end{displaymath}$

$\begin{displaymath}\frac{dM_2}{dt} = k_1 \sqrt{M_1} - k_2 \sqrt{M_2} \end{displaymath}$

Explicit Euler for M ODEs

General problem:

$\begin{displaymath}\frac{dx_1}{dt} = f_1(x_1,x_2,..,x_M,t) \end{displaymath}$

$\begin{displaymath}\frac{dx_2}{dt} = f_2(x_1,x_2,..,x_M,t) \end{displaymath}$

..

$\begin{displaymath}\frac{dx_M}{dt} = f_M(x_1,x_2,..,x_M,t) \end{displaymath}$

Approximate the derivative at the start of the step (t_n):

$\begin{displaymath}\frac{dx_i}{dt}\vert _{t_n} = \frac{x_i(n+1) - x_i(n)}{h} \end{displaymath}$

NEW NOTATION!
x_i(n) is the ith differential variable at time step n.
Explicit Euler:

x₁(n+1) = x₁(n) + h f₁(x₁(n),x₂(n),...,x_M(n),t_n)

x₂(n+1) = x₂(n) + h f₂(x₁(n),x₂(n),...,x_M(n),t_n)

..

x_M(n+1) = x_M(n) + h f_M(x₁(n),x₂(n),...,x_M(n),t_n)

Algorithm: explicit Euler

Specify initial conditions: time t₀, value of each differential variable x_i(0).
Specify final time t_f and step size h.
Initialise step counter n = 0.
do while t_n < t_f
set t_n+1 = t_n + h
if t_n+1 > t_f then
set t_n+1 = t_f
redefine h = t_f - t_n
end if
! DIFFERENT
do i=1,M
calc $f_i(n) = f_i({\bf x}(n), t_{n})$
calc x_i(n+1) = x_i(n) + h f_i(n)
end do
increment n by 1
end do
return
Vector notation: ${\bf x} = (x_1,x_2,...,x_M)$

Example: two tanks
a = 4 kg/s,
k₁ = 0.2 kg^0.5/s,
k₂ = 0.4 kg^0.5 /s.
Initial conditions:
M₁(0) = 100 kg,
M₂(0) = 225 kg.
Choose step size for initial local error of 1 kg on M₁ and M₂.
Calculation for M₁ $\Longrightarrow$ h = 10s.
Calculation for M₂ $\Longrightarrow$ h = 7.746s.

Details in section 4.4.6.1
So choose h < 7.746 s, say h = 7 s.
Results
Note minimum in M₁.
Dynamic simulation reveals interesting, possibly unexpected behaviour.

Stiff ODE systems

Large tanks and vessels
$\rightarrow$ large time constant
$\rightarrow$ want large steps.

Small volumes or fast rate processes
$\rightarrow$ small time constant
$\rightarrow$ need small steps.

Step size is limited by the shortest time constant in the process:
(a) for accuracy of the fast dynamics,
(b) for numerical stability,even if fast dynamics no longer needs to be accurate.

If a system contains widely different time constants it is said to be stiff. Usually this means a range of at least 3 orders of magnitude in the time constants.

Two tanks with recycle
Tank 1 : a small mixing device
Tank 2 : a large storage tank.

Assume total holdups in steady state. Recycle flow (F₃) is fraction r of feed F₀.
Model mole fraction of a dilute solute. Assume tanks are well stirred with mole fractions x₁ and x₂. Let T₁ and T₂ be tank residence times.
The solute mole fractions are modelled by

$\begin{displaymath}T_1 \frac{dx_1}{dt} = x_0 + r x_2 - (1 + r) x_1 \end{displaymath}$

$\begin{displaymath}T_2 \frac{dx_2}{dt} = (1 + r) (x_1 - x_2) \end{displaymath}$

Model derivation in section 4.4.6.1
Simulate with T₁ = 10s, T₂ = 2000s, r = 0.6, x₀ = 0.1 = constant.
ODEs:

$\begin{displaymath}\frac{dx_1}{dt} = 0.01 + 0.06 x_2 - 0.16 x_1 \end{displaymath}$

$\begin{displaymath}\frac{dx_2}{dt} = 0.0008 (x_1 - x_2) \end{displaymath}$

Results with explicit Euler
1. h = 1 s ( = 10% of smaller time constant)

step time x(1) x(2)

0 0.0 0.00000 0.00000

1 1.0 0.01000 0.00000

2 2.0 0.01840 0.00001

3 3.0 0.02546 0.00002

4 4.0 0.03138 0.00004

5 5.0 0.03637 0.00007

10 10.0 0.05160 0.00024

20 20.0 0.06075 0.00069

40 40.0 0.06295 0.00167

100 100.0 0.06410 0.00457

200 200.0 0.06585 0.00922

300 300.0 0.06751 0.01364

800 800.0 0.07469 0.03271

900 900.0 0.07592 0.03599

2. h = 200 s ( = 10% of larger time constant)

step time x(1) x(2)

0 0.0 0.000 0.000

1 200.0 2.000 0.000

2 400.0 -60.000 0.320

3 600.0 1865.840 -9.331

4 800.0 -57951.014 290.696

Numerical instability destroys x(1) and then x(2).

Implicit Euler: M simultaneous ODEs
Approximate the derivative at the end of the step (t_n+1):

$\begin{displaymath}\frac{dx_i}{dt}\vert _{t_{n+1}} = \frac{x_i(n+1) - x_i(n)}{h} \end{displaymath}$

Implicit Euler: solve a system of M NLAEs in M unknowns, x₁(n+1),..,x_M(n+1):

$\begin{displaymath}{\bf g}({\bf x}(n+1)) = {\bf0} \end{displaymath}$

where $g_i({\bf x}(n+1))$ $\; = x_i(n+1) - x_i(n) - h f_i({\bf x}(n+1), t_{n+1})$
Algorithm: implicit Euler

Specify initial conditions: time t₀, and each x_i(0).
Specify t_f and h.
Initialise step counter n = 0.
do while t_n < t_f
set t_n+1 = t_n + h
if t_n+1 > t_f then
set t_n+1 = t_f
redefine h = t_f - t_n
end if
! DIFFERENT FROM
! EXPLICIT EULER
solve simultaneous system of M NLAEs
$g_1({\bf x}(n+1),t_{n+1}) = 0$
...
$g_M({\bf x}(n+1),t_{n+1}) = 0$
increment n by 1
end do
return

Algorithm for NLAE solution insection 4.4.6.1
Example: two tanks with recycle.Linear ODEs and initial conditions as for explicit Euler.
h = 200 s

step time x(1) x(2)

0 0.0 0.00000 0.00000

1 200.0 0.06381 0.00880

2 400.0 0.06875 0.01707

3 600.0 0.07163 0.02460

4 800.0 0.07421 0.03144

5 1000.0 0.07655 0.03766

10 2000.0 0.08543 0.06126

15 3000.0 0.09094 0.07592

20 4000.0 0.09437 0.08504

Implicit Euler solves the problem very effectively with a large time step.
Values are slightly behind (lower than) explicit Euler with h = 1 s.
Summary

Systems of ODEs

explicit Euler algorithm is like single ODE version.
implicit Euler needs solution of simultaneous NLAE system

Stiff ODEs

characterised by wide range of time constants
explicit methods must take steps $\approx$ shortest time constant; impractical
implicit methods highly effective even with step $\gg$ shortest time constant

Section 4.4.6.1: Choosing a Step Size for More than One ODE
Section 4.4.6.1: Choosing a Step Size for More than One ODE

For each ODE we have a desired local error, say at t₀ = 0.
Calculate step size h to meet the desired error for each equation.
Choose the smallest of the step sizes calculated.
Example
ODEs:

$\begin{displaymath}\frac{dM_1}{dt} = 4 - 0.2 \sqrt{M_1} \end{displaymath}$

$\begin{displaymath}\frac{dM_2}{dt} = 0.2 \sqrt{M_1} - 0.4 \sqrt{M_2} \end{displaymath}$

with M₁(0) = 100 kg, M₂(0) = 225 kg.
Let desired local error be 1 kg for both M₁ and M₂.
Local errors are:

$\begin{displaymath}e_1 \approx (h^2/2) \frac{d^2M_1}{dt^2} \end{displaymath}$

$\begin{displaymath}e_2 \approx (h^2/2) \frac{d^2M_2}{dt^2} \end{displaymath}$

M₁:
$d^2M_1/dt^2 = d( 4 - 0.2 \sqrt{M_1} )/dt = - 0.2 d(\sqrt{M_1})/dt = -(0.1/\sqrt{M_1}) dM_1/dt\\ = -(0.1/\sqrt{M_1}) ( 4 - 0.2 \sqrt{M_1} )$ .
At t = 0 this gives
$\vert d^2M_1/dt^2\vert = \vert(0.1/\sqrt 100) ( 4 - 0.2 \sqrt 100 ) \vert = 0.02.$
We require
$e_1 \approx 1$ , or $(h^2/2) \times 0.02 = 1$ .
So h = 10 s.
M₂:
$d^2M_2/dt^2 = d( 0.2 \sqrt{M_1} - 0.4 \sqrt{M_2} )/dt = (0.1/\sqrt{M_1}) dM_1... ...( 4 - 0.2 \sqrt{M_1} ) - (0.2/\sqrt{M_2}) ( 0.2 \sqrt{M_1} - 0.4 \sqrt{M_2} )$ .
At t = 0,
$\vert d^2M_2/dt^2\vert = \vert (0.1/\sqrt 100) ( 4 - 0.2 \sqrt 100 ) - (0.2/... ... \sqrt 100 - 0.4 \sqrt 225 ) \vert = \vert 0.02 - 0.053333 \vert\\ = 0.033333$ .
We require
$e_2 \approx 1$ , or $(h^2/2) \times 0.033333 = 1$ .
So $h = \sqrt 60 = 7.746$ s.
Overall, step is restricted to be less than 7.746 s.

Two tanks with recycle: derivation of model

Steady state flow balances, no overall accumulation in tanks:

F₀ + F₃ = F₁

F₁ = F₂ + F₃

Adding equations gives F₂ = F₀ and if we say F₃ = r F₀ to define r, then

F₁ = (1+r) F₀

Component balance on solute, M_i being constant:

$\begin{displaymath}M_1 \frac{dx_1}{dt} = F_0 x_0 + F_3 x_2 - F_1 x_1 \end{displaymath}$

$\begin{displaymath}M_2 \frac{dx_2}{dt} = F_1 x_1 - F_2 x_2 - F_3 x_2 \end{displaymath}$

Express all flows in terms of F₀ and r:

$\begin{displaymath}M_1 \frac{dx_1}{dt} = F_0 (x_0 + r x_2 - (1 + r) x_1) \end{displaymath}$

$\begin{displaymath}M_2 \frac{dx_2}{dt} = F_0 ( (1 + r) x_1 - x_2 - r x_2 ) \end{displaymath}$

Introduce time constants: T₁ = M₁ / F₀ and T₂ = M₂ / F₀.
Rewrite equations (after division by F₀ and a little rearrangement) as

$\begin{displaymath}T_1 \frac{dx_1}{dt} = x_0 + r x_2 - (1 + r) x_1 \end{displaymath}$

$\begin{displaymath}T_2 \frac{dx_2}{dt} = (1 + r) (x_1 - x_2 ) \end{displaymath}$

Algorithm for NLAE solution

do i = 1,M
initialise x_i⁽⁰⁾ = x_i(n).
end do
do k = 1,max_iterations
call multi_der_fn $(M,{\bf x}^{(k-1)},t_{n+1},$ f,DFDX)
do i = 1,M
calculate g_i = x_i^(k-1) - x_i(n) - h f_i
do j = 1,M
calculate DGDX_ij = $\delta_{ij} - h$ DFDX_ij
! $\delta_{ij} = 1$ if i = j, or 0 otherwise
end do
end do
if |g_i| <nlae_tolerance for all i = 1 to M, exit
solve the linear equation system

$\begin{displaymath}DGDX \; {\bf\Delta x} = -{\bf g} \end{displaymath}$

! use solve_linear
do i=1,M
calculate $x_i^{(k)} = x_i^{(k-1)} + \Delta x_i$
end do
end do
do i = 1,M
set x_i(n+1) = x_i^(k-1)
end do

Remark:
You must supply a subroutine
multi_der_fn (x,t,f,DFDX)
which when given x and t will calculate and return f and DFDX where DFDX_ij = $\frac{\partial f_i}{\partial x_j}$ .

Section 4.4.7: Euler is not the Only Method
Section 4.4.7: Euler is not the Only Method

Euler gives exact solution if and only if exact solution is linear function of time.
More sophisticated methods match exact solution to higher orders:

quadratic
cubic
quartic
....

Can use larger step sizes for specified accuracy, but still shorter than smallest time constant.

Graph suggests two possible approaches:

1.
Secant slope is well approximated by average of slopes m₀ and m₁ at start and end of step:

$\begin{displaymath}\frac{1}{2} (\frac{dx}{dt}\vert _{t_n} + \frac{dx}{dt}\vert _{t_{n+1}}) \approx \frac{x(n+1) - x(n)}{h} \end{displaymath}$

2.
Secant slope is good approximation to slope half way across step:

$\begin{displaymath}\frac{dx}{dt}\vert _{t_{n+\frac{1}{2}}} \approx \frac{x(n+1) - x(n)}{h} \end{displaymath}$

Explicit second order methods - 1
Solve

$\begin{displaymath}\frac{dx}{dt} = f(x,t) \end{displaymath}$

Calculate f(n) = f(x(n),t(n)), slope at start of step
Explicit Euler to give first order estimate, `Stage I':
x(I) = x(n) + h f(n)
Stage I derivative estimate at end of step:
f(I) = f(x(I),t_n+1)
Stage II derivative estimate :average of f(n), f(I):
$f(II) = \frac{1}{2} (f(I) + f(n))$

Improved Euler step:
x(n+1) = x(n) + h f(II)

Explicit second order methods - 2
Another possibility - take the first order step (explicit Euler) half way across the step and estimate a derivative for the full step.

$x(I) = x(n) + \frac{1}{2} h f(n)$
Stage I derivative estimate at the half step:
$f(I) = f(x(I),t(n)+\frac{1}{2}h)$
Use this derivative estimate for the Modified Euler step:
x(n+1) = x(n) + h f(I)

Accuracy
General second order explicit method - section 4.4.7.1
Local error
$e_n = \frac{1}{6} (hk)^3 x(n)$ for linear ODE

Case 1: for initial local error = 0.01 x(0)
Euler step: $\frac{1}{2} (hk)^2 = 0.01$ , or
$hk = \sqrt{0.02} = 0.1414$
2nd order step: $\frac{1}{6} (hk)^3 = 0.01$ , or
$hk = \sqrt[3]{0.06} = 0.3915$
Two function evaluations per step, but more than twice the step size ( $\times$ 2.77).
Case 2: for initial local error = 0.0001 x(0)
Euler step: $\frac{1}{2} (hk)^2 = 0.0001$ , or
$hk = \sqrt{0.0002} = 0.01414$
2nd order step: $\frac{1}{6} (hk)^3 = 0.0001$ , or
$hk = \sqrt[3]{0.0006} = 0.08434$
2nd order has nearly six times the step, so nearly three times as efficient overall.

Stability

$1 - h k + \frac{1}{2} (h k)^2$ has a minimum at h k = 1, value $= \frac{1}{2}$
For h k > 2, $1 - h k + \frac{1}{2} (h k)^2$ exceeds 1.
Stability limit at h k = 2, same as explicit Euler.
Similar limits for all higher order explicit methods (usually limiting h k in a range around 2 - 4).

Simultaneous ODEs

Apply each stage in the time step to all equations at the same time.
Example

$\begin{eqnarray*}\frac{dx_1}{dt} = f_1(x_1,x_2) \\ \frac{dx_2}{dt} = f_2(x_1,x_2) \end{eqnarray*}$

Modified Euler:

Calculate

$\begin{eqnarray*}x_1(I) = x_1(n) + \frac{1}{2} h f_1(n) \\ x_2(I) = x_2(n) + \frac{1}{2} h f_2(n) \\ \end{eqnarray*}$

where $f_i(n) = f_i({\bf x}(n))$ .
Stage I derivative:

$\begin{eqnarray*}f_1(I) = f_1({\bf x}(I)) \\ f_2(I) = f_2({\bf x}(I)) \end{eqnarray*}$

Full step:

$\begin{eqnarray*}x_1(n+1) = x_1(n) + h f_1(I) \\ x_2(n+1) = x_2(n) + h f_2(I) \end{eqnarray*}$

Runge-Kutta methods

Extend the idea of the 2nd order explicit methods. Exactly match higher degree polynomials by making more estimates of the derivatives at suitably chosen fractions of the time step and averaging them with appropriate weights.
Up to order 4 or 5 are commonly used.
Number of function evaluations per step is usually equal to or greater than the order. But overall efficiency is better than 2nd order.
Some methods contain built-in error estimates, for only a little extra work (1 or 2 more function evaluations per step.)
Explicit RK methods all have limited stability ( $hk \leq 3 - 4$ ).

Multistep methods

Another way of getting higher order polynomial approximations is to use information from previous steps.
Use only one function evaluation or (for implicit methods) NLAE solution at each time step, like Euler.

Adams-Bashforth methods: explicit
x(n+1) = x(n) + h(a₀ f(n) + a₁ f(n-1) + .. + a_N f(n-N) )
Adams-Moulton methods: implicit
x(n+1) = x(n) + h(a₀ f(n+1) + a₁ f(n) + .. + a_N f(n+1-N) )
Backward Difference Formula (BDF) methods: implicit
$x(n+1) = \alpha_1 x(n) + \alpha_2 x(n-1) + ... + \alpha_N x(n-N) + h \beta f(n+1)$

BDF methods are popular. Good stability. More efficient than RK, BUT they need help getting started. Euler must be used initially, then 2nd order and so on. Also, changing step size is less flexible than with RK.

Summary

Higher order methods than Euler allow
longer steps,
greater overall efficiency

especially at high accuracy.
Second order methods use central derivative approximations
Runge-Kutta methods make several function evaluations per time step
Multistep methods use past information and keep to one function evaluation per time step.
Multistep methods are slightly less flexible and need to start cautiously.

Section 4.4.7.1: Demonstration that Improved and Modified Euler are Second Order Methods

Section 4.4.7.1: Demonstration that Improved and Modified Euler are Second Order methods

Second order means:
the method matches the exact solution if the exact solution is a quadratic:
x(t) = x(n) + a₁ (t - t_n) + a₂ (t - t_n)².
For such a function,
$\frac{dx}{dt} = f(t) = a_1 + 2 a_2 (t - t_n)$ .
See how the proposed methods work on $\frac{dx}{dt} = f(t)$ .

1.
f(n) = f(t_n) = a₁
2.
Take a first order step for a fraction p of the full step (p = 1 or $\frac{1}{2}$ for Improved or Modified Euler)
x(I) = x(n) + p h f(n) = x(n) + a₁ p h;
t(I) = t_n + p h.
3.
Stage I derivative estimate at fraction p across step
f(I) = f(t(I)) = a₁ + 2 a₂ (t(I) - t_n) = a₁ + 2 a₂ p h.
4.
For Stage II derivative take a weighted average of f(I) (assuming weight q) and f(n) (weight (1-q)):
$f(II) = q f(I) + (1-q) f(n) = q(a_1 + 2 a_2 p h) + (1-q) a_1 , \\ f(II) = a_1 + 2 a_2 p q h$
5.
Full step using Stage II derivative f(II)
x(n+1) = x(n) + h f(II) = x(n) + a₁ h + 2 a₂ p q h².

Compare this with the exact solution

x(t_n+1) = x(n) + a₁ (t_n+1 - t_n) + a₂ (t_n+1 - t_n)²

x(t_n+1) = x(n) + a₁ h + a₂ h²

For an exact match a second order method must have

2 p q = 1

Thus if p=1 (Improved Euler with full first order step),
we need $q=\frac{1}{2}$ (average of start and end of step derivative estimates).
If $p=\frac{1}{2}$ (Modified Euler with half first order step),
we need q=1 (use the mid point derivative for the full step).
These values of q are what we assumed before, so the methods are indeed second order.

Performance of Second Order Explicit Methods on a Linear ODE

Apply the general second order explicit method (with parameters p and $q = \frac{1}{2p}$ ) to the linear ODE

$\begin{displaymath}\frac{dx}{dt} = - k x\end{displaymath}$

1.
f(n) = f(x(n)) = - k x(n)
2.
First order step to fraction p of the full step
x(I) = x(n) + p h f(n) = x(n) - p h k x(n)
3.
Stage I derivative estimate
f(I) = f(x(I)) = - k x(I) = - k x(n) + p h k² x(n)
4.
Stage II derivative (weighted average):
$f(II) = q f(I) + (1-q) f(n) = \frac{1}{2p} f(I) + (1- \frac{1}{2p}) f(n) , \\ ... ...II) = - k x(n) + \frac{1}{2p} p h k^2 x(n) = - k x(n) + \frac{1}{2} h k^2 x(n)$
5.
Full step:
$x(n+1) = x(n) + h f(II) = x(n) - h k x(n) + \frac{1}{2} (h k)^2 x(n)$ .

Thus

$\begin{displaymath}\frac{x(n+1)}{x(n)} = 1 - h k + \frac{1}{2} (h k)^2 \end{displaymath}$

Local error

$\begin{displaymath}e_n = x(n+1) - x(n) exp(- h k) = [ 1 - h k + \frac{1}{2} (h k... ... h k + \frac{1}{2!} (h k)^2 - \frac{1}{3!} (h k)^3 ...) ] x(n) \end{displaymath}$

$\begin{displaymath}\vert e_n \vert \approx \frac{1}{6} (h k)^3 x(n) \end{displaymath}$

In general

$\begin{displaymath}\vert e_n \vert \approx \frac{1}{6} h^3 \vert\frac{d^3x}{dt^3}\vert \end{displaymath}$

Section 4.4.8: Applications of ODE's
Section 4.4.8: Applications of ODE's

Control

In the gravity flow tank we might control the outflow with a flow controller.
Without control:

$\begin{displaymath}F_2 = k \sqrt{M} = c \sqrt{p_M} \end{displaymath}$

where $p_M = \frac{Mg}{A}$ ,
A is the area of the tank.
With control:

$\begin{displaymath}F_2 = \sqrt{\frac{p_M}{\frac{1}{c^2} + \frac{1}{c_v^2}}} = \sqrt{\frac{M}{\frac{1}{k^2} + \frac{1}{k_v^2}}} \end{displaymath}$

k_v is a parameter which varies with valve position:

$\begin{displaymath}0 \leq k_v \leq k_{v,max} \end{displaymath}$

Flow Controller

measure flow
compare with desired or setpoint flow, F_s.
generate output signal
transmit to valve actuator
move valve to better position: if flow
too high, partly close the valve
too low, partly open the valve

Simplest type: proportional controller.

Equations:
error = set point - measured value
output = proportional gain * error
or

$\begin{eqnarray*}e = 100 (F_s - F_2) / F_{range} \\ y = K e \end{eqnarray*}$

Usually scale signals 0 - 100 % of range.
Assume valve 50% open at zero output signal, and linear response.
Might think

$\begin{displaymath}k_v = \frac{(y+50)}{100} k_{v,max} \end{displaymath}$

But valve cannot close beyond 0 or open beyond 100%:

$\begin{displaymath}k_v = k_{v,1} \equiv \frac{(y+50)}{100} k_{v,max} \end{displaymath}$

if 0 < k_v,1 < k_v,max
k_v = 0 if $k_{v,1} \leq 0$
k_v = k_v,max if $k_{v,1} \geq k_{v,max}$
Model contains discontinuities!
k_v depends on the flow F₂ through the error e and output y. If a single ODE is used:

$\begin{displaymath}\frac{dM}{dt} = F_1 - F_2 (= f(M)) \end{displaymath}$

we need to solve an equation iteratively for F₂ given M at each function evaluation:

$\begin{displaymath}F_2 = \sqrt{\frac{M}{\frac{1}{k^2} + \frac{1}{(k_v(F_2))^2}}} \end{displaymath}$

It might be better to leave this as an algebraic equation and treat the system as a DAE (see section 4.4.9).
So F₂, k_v, e, y etc. would be kept explicitly in the model and an equation would be added to define each extra variable.

Level Controller
Instead of trying to control the flow, try to control the level.(For this example, assume level $\propto$ holdup.)
Similar equations, but

e = 100 (M_s - M) / M_range

where M_s is the `holdup setpoint', and M_range is the range of holdup for 0 - 100 % measured value signal.
Also, y = K e,
$k_{v,1} = \frac{(y+50)}{100} k_{v,max}$ and
k_v = f_k(k_v,1), as before.
Direct evaluation of F₂ given M in this case - no iteration.

Example

k = 0.2 kg^0.5/s, k_v,max = 1.0 kg^0.5/s.
Assume M_range = 500 kg, M_s = 300 kg, and controller gain K = -2.
Minus because low level $\rightarrow$ positive
error $\rightarrow$ negative output, which is needed to close the outflow valve.
Initial condition: M₀ = 100 kg
Feed flow: F₁ = 4 kg/s until t = 180 s, then F₁ = 2 kg/s

Results
Results are tabulated at the end of this section.
Controller fails to reach set point (300 kg) but does achieve steady state. This is normal for `proportional only' controllers. If the current steady state is not exactly the one designed for, a continuing offset is needed give a non-zero error and maintain the control action.
In this case, zero output gives 50% valve opening: k_v = 0.5 kg^0.5/s
For 300 kg holdup, this needs an outlet flow of $\sqrt{\frac{300}{\frac{1}{0.2^2} + \frac{1}{0.5^2}}} = \sqrt{\frac{300}{29}} = 3.216$ kg/s
Since the feed flow is finally 2 kg/s, we need a smaller final flow, more closed valve (18.2%) and hence positive 'error': set point > measured holdup.

Thermodynamics

The independent variable need not be time!
ODE solving methods are useful in thermodynamics, where the independent variable may be a thermodynamic variable.
Example: Pressure dependence of enthalpy at constant temperature

dh = T ds + v dp

h - specific enthalpy (J kmol^-1),
T - temperature (K),
s - specific entropy (J kmol^-1 K^-1),
v - specific volume (m³ kmol^-1),
p - pressure (Nm^-2)
Find change in enthalpy from a standard pressure p₀ (e.g. 1 atm = 1.013 $\times$ 10⁵ Nm^-2), to required pressure p₁.
Treat T and p as independent variables:

$\begin{displaymath}dh = T \left(\frac{\partial s}{\partial T}\right)_p dT + T \left(\frac{\partial s}{\partial p}\right)_T dp + v dp \end{displaymath}$

Thus at constant T = T₁,

$\begin{displaymath}\left(\frac{\partial h}{\partial p}\right)_T = T_1 \left(\frac{\partial s}{\partial p}\right)_T + v \end{displaymath}$

Maxwell relations give

$\begin{displaymath}\left(\frac{\partial s}{\partial p}\right)_T = - \left(\frac{\partial v}{\partial T}\right)_p \end{displaymath}$

So

$\begin{displaymath}\left(\frac{\partial h}{\partial p}\right)_T = - T_1 \left(\frac{\partial v}{\partial T}\right)_p + v \end{displaymath}$

Equation of state

We would like to get

$\begin{displaymath}v = \Phi(T,p)\end{displaymath}$

from the equation of state (EOS), with T = T₁. If EOS is in the above form, then no problem - just evaluate whenever volume is needed. But usually EOS is pressure explicit:

$\begin{displaymath}p = \Pi(T,v)\end{displaymath}$

e.g. Redlich-Kwong EOS

$\begin{displaymath}p = \frac{RT}{v-b} - \frac{a}{T^{\frac{1}{2}} v(v+b)} \end{displaymath}$

(a,b constants for given chemical).
We could solve iteratively for v at each value of p in the integration. But this means many iterative calculations between p₀ and p₁.

A smarter way
Use implicit differentiation of the pressure-explicit EOS to get expressions for $(\partial v/\partial p)_T$ and $(\partial v/\partial T)_p$ as functions of specific volume at temperature T₁:

$\begin{displaymath}(\partial v/\partial p)_T = \pi_1(T_1,v) \end{displaymath}$

$\begin{displaymath}(\partial v/\partial T)_p = \pi_2(T_1,v) \end{displaymath}$

Do just one iterative calculation to get v₀ at pressure p₀. (h₀ is the enthalpy at T₁ and p₀.) Solve the ODE system

$\begin{displaymath}\frac{dh}{dp} = - T_1 \pi_2( T_1,v) + v \end{displaymath}$

$\begin{displaymath}\frac{dv}{dp} = \pi_1(T_1,v)\end{displaymath}$

for h and v as functions of p at constant T₁.
Section 4.4.8.1 shows how to get functions $\pi_1$ and $\pi_2$ for Redlich-Kwong.

Summary

Control
Dynamic simulation can show effect of control strategy and tuning constants (gain K).Latter not covered in this section. Proportional-only control leads to offset:set point is not usually achieved.

Thermodynamics
Coupled ODEs in a thermodynamic variable (such as pressure) can reduce the computational effort in solving for thermodynamic functions (e.g. enthalpy).

Table of output: Holdup controller example

step time Holdup Outlet flow % control signal

0 0.00 100.0000 0.0000 0.0000

1 10.00 140.0000 0.0000 0.0000

2 20.00 178.6650 0.1954 1.4660

3 30.00 208.2578 1.5985 13.3031

4 40.00 228.7844 2.2156 21.5138

5 50.00 244.8290 2.5444 27.9316

6 60.00 258.2513 2.7553 33.3005

7 70.00 269.8923 2.9068 37.9569

8 80.00 280.2064 3.0238 42.0825

9 90.00 289.4720 3.1183 45.7888

10 100.00 297.8756 3.1972 49.1502

11 110.00 305.5505 3.2647 52.2202

12 120.00 312.5971 3.3235 55.0388

13 130.00 319.0932 3.3752 57.6373

14 140.00 325.1017 3.4213 60.0407

15 150.00 330.6740 3.4627 62.2696

16 160.00 335.8534 3.5001 64.3414

17 170.00 340.6767 3.5341 66.2707

18 180.00 345.1755 3.5651 68.0702

19 190.00 339.3776 3.5250 65.7510

20 200.00 324.6830 3.4181 59.8732

21 210.00 311.0609 3.3109 54.4243

22 220.00 298.5167 3.2030 49.4067

23 230.00 287.0576 3.0945 44.8230

24 240.00 276.6902 2.9854 40.6761

.. ...... ........ ...... .......

27 270.00 252.1062 2.6644 30.8425

30 300.00 236.6112 2.3888 24.6445

33 330.00 228.1177 2.1995 21.2471

36 360.00 223.9512 2.0938 19.5805

39 390.00 222.0400 2.0421 18.8160

42 420.00 221.1926 2.0185 18.4771

45 450.00 220.8229 2.0080 18.3292

48 480.00 220.6627 2.0035 18.2651

54 540.00 220.5636 2.0006 18.2255

60 600.00 220.5452 2.0001 18.2181

66 660.00 220.5418 2.0000 18.2167

72 720.00 220.5411 2.0000 18.2165

Section 4.4.8.1: Implicit Differentiation of Redlich-Kwong equation
Section 4.4.8.1: Implicit Differentiation of Redlich-Kwong equation

In general a pressure explicit EOS is

$\begin{displaymath}p = \Pi(T,v)\end{displaymath}$

Differentiate the EOS with respect to p at constant T.

$\begin{displaymath}1 = \frac{\partial \Pi}{\partial v} \left(\frac{\partial v}{\partial p}\right)_T \end{displaymath}$

so

$\begin{displaymath}\left(\frac{\partial v}{\partial p}\right)_T \equiv \pi_1(T,v) = \frac{1}{\frac{\partial \Pi}{\partial v}} \end{displaymath}$

For $(\partial v/\partial T)_p$ , differentiate with respect to T at constant p.

$\begin{displaymath}0 = \frac{\partial \Pi}{\partial T} + \frac{\partial \Pi}{\partial v} \left(\frac{\partial v}{\partial T}\right)_p \end{displaymath}$

so

$\begin{displaymath}\left(\frac{\partial v}{\partial T}\right)_p \equiv \pi_2(T,v... ...i}{\partial v}} = - \pi_1(T,v) \frac{\partial \Pi}{\partial T} \end{displaymath}$

For the Redlich-Kwong equation

$\begin{displaymath}\Pi(T,v) = \frac{RT}{v-b} - \frac{a}{T^{\frac{1}{2}} v(v+b)} \end{displaymath}$

$\begin{displaymath}\frac{\partial \Pi}{\partial v} = \frac{-RT}{(v-b)^2} + \frac{a}{T^{\frac{1}{2}}} \frac{(2 v + b)}{v^2(v+b)^2} \end{displaymath}$

$\begin{displaymath}\pi_1(T,v) = \frac{1}{\frac{-RT}{(v-b)^2} + \frac{a}{T^{\frac{1}{2}}} \frac{(2 v + b)}{v^2(v+b)^2}} \end{displaymath}$

$\begin{displaymath}\frac{\partial \Pi}{\partial T} = \frac{R}{v-b} + \frac{a}{2 T^{\frac{3}{2}} v(v+b)} \end{displaymath}$

$\begin{displaymath}\pi_2(T,v) = - \frac{ \frac{R}{v-b} + \frac{a}{2 T^{\frac{3}{... ...)^2} + \frac{a}{T^{\frac{1}{2}}} \frac{(2 v + b)}{v^2(v+b)^2}} \end{displaymath}$

Though rather complex to look at, these expressions are easy to evaluate given v. Non-integer powers of T = T₁ would be calculated before the integration over p, to avoid repetition of the more expensive computations.

Section 4.4.9: DAE Systems
Section 4.4.9: DAE Systems

The Flow Controller revisited

So far, ODEs only:

$\begin{displaymath}\frac{d{\bf x}}{dt} = {\bf f}({\bf x},t) \end{displaymath}$

E.g.

$\begin{displaymath}\frac{dM}{dt} = f(M) \end{displaymath}$

If f(M) involves non-trivial calculation of other variables we may be better to `open out' the problem and make all these variables explicitly part of the formulation.
E.g. flow controller:
f(M) = F₁ - F₂ where F₂ depends on M and k_v, the valve opening.
But k_v depends on F₂ because of the control algorithm.

DAE formulation of FC

Full set of equations in section 4.4.9.1
Instead of 1 ODE (mass balance) in 1 variable (M),
we have 12 equations in 12 variables.
1 equation is an ODE (the mass balance).
The other 11 are algebraic equations.
The 6 simple specification equations are included to show the full structure that is possible, though we might prefer to regard the corresponding variables as constants or explicit functions of time instead.
Terminology

General DAE system:
$\frac{d{\bf x}}{dt} = {\bf f}({\bf x,y},t)$ : M_D ODEs
${\bf g}({\bf x,y},t) = {\bf0}$ : M_A algebraic equations.

x: variables whose time derivatives appear in the model,
y: variables whose time derivatives do not appear in the model.
x are called differential variables,
y are called algebraic variables.
We have M_D differential and M_A algebraic variables.
For a DAE system to be easy to solve,
it must be possible to solve the AEs

$\begin{displaymath}{\bf g} = {\bf0} \end{displaymath}$

for y, if fixed values of x and t are given.
If not, `high-index' difficulties occur.
Will not be discussed in this course.
Solution of high-index problems remains an active research topic.
It's always best to make sure that a DAE system you formulate satisfies the above condition.

Flow Controller example

M is a differential variable because $\frac{dM}{dt}$ appears in the model.
The other 11 variables are algebraic variables: F₁, F₂, k, k_v, $\phi$ , k_v,max, y, K, e, F_s and F_range.
If M is given the AEs can be used to solve for the algebraic variables.
For example:

1.
Guess F₂.
2.
Calculate e, then y, then $\phi$ , then k_v; substituting the specification equations as required.
3.
Calculate the residual $F_2 - \sqrt{\frac{M}{1/k^2 + 1/k_v^2} }$
4.
Use a secant method (say) to get a new estimate for F₂.
5.
Repeat from 2. until converged.
This loop gives all the algebraic variables except F₁ which is easily calculated from its specification equation.
Explicit solution

Explicit Euler:
Discretise the ODE as

M(n+1) = M(n) + h ( F₁(n) - F₂(n) )

Needs values for all the algebraic variables at t_n.
Integration gives only the differential variable(s) at t_n+1.
Then solve the AEs for the algebraic variables at t_n+1 using the value calculated for M(n+1) by explicit Euler.

Solution is no different from treating the system as an ODE:

M(n+1) = M(n) + h f(n)

where $f(n) \equiv f(M(n))$ is calculated in exactly the same way as F₁(n) - F₂(n), using for example the iterative scheme on the previous slide.

Initial conditions

If the flow controller is treated as an ODE the initial condition is M(0) at time t₀.
This must remain true if the problem is treated as a DAE system.
General Rule for ICs:
One IC for each differential variable.
A safe choice is to give the initial value of each differential variable (as here). In other cases this may be inconvenient - see next section.
What about the algebraic variables?
Consistent initial values must be found by solving the AEs at time t₀ with specified values ${\bf x}(0)$ for the differential variables.
I.e. solve

$\begin{displaymath}{\bf x} - {\bf x}(0) = 0 \end{displaymath}$

$\begin{displaymath}{\bf g}({\bf x,y},t_0) = 0 \end{displaymath}$

This is equivalent to replacing the ODEs with IC specifications.

DO NOT GIVE INDEPENDENT
INITIAL CONDITIONS
FOR THE ALGEBRAIC VARIABLES.
However, initial guesses may be needed for the algebraic variables when solving the initial point equations above, e.g. by Newton's method.

Implicit solution

Implicit Euler:
Discretise the ODE as

M(n+1) = M(n) + h ( F₁(n+1)) - F₂(n+1) )

M(n+1), F₁(n+1) and F₂(n+1) are all unknown. So it makes sense to solve the discretised ODE simultaneously with the AEs at t_n+1.
The unknowns in the AEs are the differential and algebraic variables at time t_n+1:

$\begin{displaymath}{\bf x}(n+1), {\bf y}(n+1) \end{displaymath}$

This is different from the method used if we regard the system as just an ODE.
If the flow controller is treated as an ODE

$\begin{displaymath}\frac{dM}{dt} = f(M) \end{displaymath}$

and solved with implicit Euler, we must solve

M(n+1) = M(n) + h f(n+1)

for M(n+1), given that f is a function of M.
Solving for M(n+1) involves an iterative loop, if f is non-linear.
Each iteration involves calculation of f at an iterated value M^(k) for M(n+1).
But f = F₁ - F₂ and, in the flow control example, iterative solution for F₂ is required at a given M.
Thus we have an inner iterative loop for F₂,
and an outer iterative loop for M.
Which is better?
Nested loops each for one unknown,
or a single loop for 12 variables
(though a partitioning scheme can easily be devised which reduces the number of simultaneous unknowns in the implicit DAE solution to 2: M and F₂.)

An open question.
Probably depends on the problem.
Experimental results

The flow controller was solved for 200 time steps of 10 s, by which time (2000 s) steady state had been reached.
(Feed F₁ and initial holdup M₀ as in section 7 level control example).
Solving as a DAE needed 418 simultaneous linear equation solutions for the 12 variables, i.e. 418 Newton iterations of the full system.
(This could have been reduced to a 2 variable system with special programming.)
3 iterations per step to t=540 s, 2 iterations per step to 1610 s, 1 iteration per step thereafter.
Solving as an ODE took 464 outer loop iterations:
3 outer to 710 s, 2 outer to 1910 s, then 1 outer. All inner loops converged in 3 secant iterations.
Conclusion

A simple flow control problem was seen to be expressible either as a
single ODE or as a DAE system.
For explicit solution of the ODEs there is no practical difference in the solution for the two formulations.
For implicit solution of the ODEs
the ODE formulation tends to give nested iterative loops, but
the DAE formulation allows fully simultaneous solution.
There must be the same number of initial conditions as ODEs, and they may safely be given as initial values for the differential variables.

Section 4.4.9.1: DAE Formulation of Flow Controller problem<

Section 4.4.9.1: DAE Formulation of Flow Controller problem
The flow controller model presented in lecture 7 can be formulated as a DAE system of 12 equations as follows.

$\frac{dM}{dt} = F_1 - F_2$ ODE: material balance

$F_1 - \theta(t) = 0$ specify feed flow F₁ as function of time

$F_2 - \sqrt{\frac{M}{1/k^2 + 1/k_v^2} } = 0$ outflow - pressure drop equation

k - (0.2) = 0 specify pipe flow constant

$k_v - \phi k_{v,max} = 0$ valve characteristic equation

k_v,max - (1.0) = 0 specify fully-open valve flow constant

$g_{\phi}(\phi,y) = 0$ valve opening as function of controller output

y - K e = 0 proportional controller algorithm

$e - 100 \frac{(F_s - F_2)}{F_{range}} = 0$ controller error definition

K - (1.2) = 0 specify proportional gain

F_s - (2) = 0 specify flow set point

F_range - (5) = 0 specify flow range

Constants in brackets could equally well take other values, or indeed some could be specified functions of time (e.g. F_s, perhaps K).
$g_{\phi}(\phi,y)$ is a `multiple choice' function:
$\phi = (y+50)/100$ if -50 < y < +50
$\phi = 1$ if y > +50
$\phi = 0$ if y < -50
This ensures that $\phi$ , the fractional valve opening, is always in the range 0 to 1.
There are 1 differential and 11 algebraic equations. 6 of the algebraic equations are simple specifications (constant or function of time):
F₁, k, k_v,max, K, F_s, F_range.
There are 12 variables (`in order of appearance'):

M tank holdup

F₁ feed flow

F₂ outlet flow

k pipe flow constant

k_v valve characteristic

$\phi$ valve open fraction

k_v,max valve flow constant (fully open)

y controller output

K controller proportional gain

e controller error signal

F_s flow set point

F_range flow range for controller

Consistent initial values
Replace ODE (1) by

x = 1

and solve with

$\begin{eqnarray*}y_1 y_2 - x^2 = 0 \\ y_1 - x y_2^3 - y_2 = 0 \end{eqnarray*}$

This is of course the same as solving

$\begin{eqnarray*}y_1 y_2 - 1= 0 \\ y_1 - y_2^3 - y_2 = 0 \end{eqnarray*}$

If using Newton we need initial guesses for y₁ and y₂, say
y₁⁽⁰⁾ = 1, y₂⁽⁰⁾ = 1
These are NOT initial conditions for the y variables, because they do not satisfy both the AEs. Applying Newton we get

y₁ 1.0000 1.2000 1.2710 1.2720 = y₁(0)

y₂ 1.0000 0.8000 0.7860 0.7862 = y₂(0)

First time step by implicit Euler
Solve simultaneously the discretised ODE and the AEs:

$\begin{eqnarray*}x(1) - x(0) - h (- x(1)/y_1(1)) = 0 \\ y_1(1) y_2(1) - x(1)^2 = 0 \\ y_1(1) - x(1) y_2(1)^3 - y_2(1) = 0 \end{eqnarray*}$

Suppose h = 0.01. Then solve

$\begin{eqnarray*}x(1) - 1 - 0.01 (- x(1)/y_1(1)) = 0 \\ y_1(1) y_2(1) - x(1)^2 = 0 \\ y_1(1) - x(1) y_2(1)^3 - y_2(1) = 0 \end{eqnarray*}$

Note that the only `time step 0' value to appear in these equations is x(0) = 1. However, we can also use the y₁(0), y₂(0) values as initial guesses to solve these equations. Using Newton the step converges after 1 iteration starting from the t=0 values as guesses.

x 1.0000 0.9921

y₁ 1.2720 1.2578

y₂ 0.7862 0.7825

Section 4.4.10: Towards Process Simulation
Section 4.4.10: Towards Process Simulation

Build a material and energy balance model for a tank with water flows in and out.
Tank is assumed to have non-constant cross section: sloping sides.
Gravity outflow model should be modified, and expressed in terms of level instead of holdup. Therefore we need tank geometry equations.
Energy balance introduces energy and enthalpy. Temperature is a more convenient quantity to specify, and to express thermal equilibrium. Hence thermodynamic and physical property equations are needed.
See section 4.4.10.1, `Material and energy balance for a tank'.

Modelling strategy
Write equations which describe the process. Aim for a logical grouping of equations:

unit by unit
within units, e.g.

Fundamental balances (conservation)
Specific model for unit: `rate equations', geometry,assumptions (well-stirred..)
Thermodynamic relations and
definitions
Physical property correlations
Specifications (`process constraints')

Choosing specifications
Also tracking whether we have suitable modelling equations.
For each equation

note any new variables
add to list of `unsolved' variables
assume the equation solves for one of the unsolved variables: tick it off and remove it from the `unsolved' list.

Choose something sensible. If there are no suitable `unsolved' candidates, go back to where there was a choice, make a different choice and repeat from there. If making a different choice further back is not possible, consider whether the new equation repeats earlier information.
When all modelling equations have been written there will generally be some `unsolved' variables left. This is because there are usually more variables than equations in the basic model of a unit. The remaining degrees of freedom in the model can be satisfied by specifying the variables which are left on the `unsolved' list.
Dynamic DAE models

For a dynamic model you get a DAE system.
In counting equations and variables exclude the independent variable (time).
Apply the following `ticking off' rule for the ODEs in a DAE system:
For an ODE $\frac{dx_i}{dt} = f({\bf x,y},t)$
always tick off the differential variable x_i.

This ensures that the specifications obtained from the final unsolved variable list do not cause a difficult (high index) DAE system.

Initial conditions

The ticking-off rule for ODEs also means that if we replace the ODEs with `standard' initial conditions at t₀, the AEs can still be ticked off in exactly the same way, because if we replace

$\begin{displaymath}\frac{d{\bf x}}{dt} = {\bf f}({\bf x,y},t) \end{displaymath}$

with

$\begin{displaymath}{\bf x} - {\bf x}(0) = 0 \end{displaymath}$

then the latter standard ICs must be ticked off for the differential variables x.
It is always safe to use standard ICs for a well-behaved DAE, but not always convenient.
E.g. do we really know the initial holdup and internal energy in the tank? Liquid level and temperature might be preferable for ICs.

Alternative ICs

In a DAE system it may be acceptable to give ICs for some of the algebraic variables instead of differential variables, as long as

the number of ICs is the same as the number of ODEs, and
the variables for which ICs are given are linked to the differential variables through the AE part of the model.

Two checks

Having constructed a model, including proposed specifications and initial conditions, you should perform two checks.

Assume fixed values of the differential variables x.
Ensure that the AEs can be solved for the algebraic variables y.
If your ICs are not the standard ones ${\bf x} - {\bf x}_0 = {\bf0}$ then also ensure that for the proposed alternative ICs, the AEs can be solved for the remaining variables, which now include both some algebraic and some differential variables.

Check 1 verifies that the DAE system with standard ICs is well-behaved (not high index).
Check 2 verifies that the proposed ICs are in some sense equivalent to the standard ICs for x.
In the event of a problem ..
If check 1 fails, the specifications and perhaps the modelling equations themselves need to be re-examined.
If check 1 succeeds but check 2 fails, you need different ICs.
The standard ICs can always be used.

An important distinction

It is important not to confuse specifications and initial conditions. Specifications are valid throughout the simulation: for example F₁, T₁, p. ICs are valid only at the starting point t₀: for example L, T.
Never try to give ICs for the same variables as specifications.

Discussion

The tank with material and energy balances is perhaps the simplest non-trivial dynamic process model.
General modelling must also take account of composition dynamics, where there is more than one chemical component.
Thus we should also have equations for the individual component balances. The `well mixed' assumption requires extra equations equating concentrations in the outlet stream to those in the tank.
More complex models (e.g. for an enclosed vessel) should introduce the possibility of multiple phases (e.g. vapour + liquid), with some way of relating them (e.g. phase equilibrium or a mass transfer model).
Then there is an endless range of possible reactor models!
DAE formulation fits well with the `equation-based approach' to process simulation, which at steady state views the process model as a system of NLAEs:

$\begin{displaymath}{\bf g}({\bf y}) = 0 \end{displaymath}$

It is natural to extend such simulators to dynamics by adding discretisations of ODEs. These are often implicit, to cope with potentially stiff dynamic flowsheet simulation problems. Hence the DAE approach.

Conclusion

`Sequential modular simulators' view the process as a collection of unit operations. Each unit is solved independently in sequence (or nowadays sometimes in parallel). This approach may lend itself more to the ODE formulation for some (simple) units, and a DAE approach in others.
The nested iterations of the ODE method may be slower, but also tend to converge more reliably, than the simultaneous solution of the DAE method.

Summary

Dynamic process simulation $\rightarrow$ DAE systems
Use logical techniques for constructing models:

group equations by unit, and by type within a unit operation,
notionally solve each equation for a designated variable; always `solve' an ODE for its differential variable.
Get a sensible set of specifications from unsolved variables.

For DAEs the ICs need not be the differential variables.
Check that the model you develop is well behaved

with standard ICs (differential variables)
with the proposed ICs, if different.

Section 4.4.10.1: Material and energy balance for a tank
Section 4.4.10.1: Material and energy balance for a tank

Assume water properties.

Equation name

Equation

Variables (units)

Balances

Material balance $\frac{dM}{dt} = F_1 - F_2$ M holdup (kg)

F₁ feed flow (kg/s)

F₂ outlet flow (kg/s)

Energy balance $\frac{dU}{dt} = h_1 F_1 - h_2 F_2$ U internal energy of contents (kJ)

(note 1) h₁ feed specific enthalpy (kJ/kg)

h₂ outlet specific enthalpy (kJ/kg)

Tank Model

gravity flow model $vF_2 - 0.005 \sqrt{L} = 0$ v specific volume ( m³/kg)

L level of contents (m)

`well mixed' model T - T₂ = 0 T tank temperature (⁰C)

T₂ outlet temperature (⁰C)

mass - volume rel'n M v - V = 0 V volume of contents (m³)

volume - level rel'n V - 0.2 ((1 + 0.8 L)³ - 1) = 0

(note 2)

Thermo and definitions

energy - enthalpy rel'n u - h + pv = 0 u specific internal energy (kJ/kg)

h specific enthalpy (tank) (kJ/kg)

p pressure (kPa) (note 3)

specific energy def'n U - M u = 0

Phys props

enthalpy correl'n h₁ - 4.1846 T₁ = 0 T₁ feed temperature (⁰C)

(feed)

enthalpy correl'n h - 4.1846 T = 0

(tank)

enthalpy correl'n h₂ - 4.1846 T₂ = 0

(outlet)

density correl'n (998.23 - 0.255 (T - 20)) v - 1 = 0

Specifications

feed flow spec $F_1 - \left(2 + \frac{0.02 t}{1 + 0.01 t}\right) = 0$

feed temp. spec T₁ - 20 = 0

pressure spec p - 100 = 0

Notes:

1.
rate of change of energy holdup = net enthalpy flow (in - out)
2.
tank assumed to have sloping sides (radius a linear function of level).
3.
u = h - pv: assume p is in kPa (1 bar = 100 kPa) and v in m³/kg; then units of u and h are consistent (kJ/kg).

Section 4.5: Case Study
Section 4.5: Case Study

4.5.1 Modelling a Simple Control System

4.5.2 Modelling a Vaporiser

4.5.3 Exercise: Modelling the vaporiser with a control system

Section 4.5.1: Modelling a Simple Control System
Section 4.5.1: Modelling a Simple Control System
This is essentially the same example used to illustrate the modelling language and the solution of DAEs.
dx/dt = (u-x) / T₁ dy/dt = (x-y) / T₂ u = mu (s - y)
The initial conditions are: x = 1 ; y = 0 at t=0 This represents a process consisting of two lags in series, with its out put y regulated at setpoint s with a proportional controller with gain mu.
You can copy the model here generate a spreadsheet from it.
Notice the form of the simple proportional controller equation. If we wanted to model a P+I controller this would require a further differential equation to model the integral action.
We can write the equations for a P+I controller with input measurement y and output u as follows:

e = s - y ; ! error u = mu * (e + ei/ti) ! controller algebraic equation .ei = e ! o.d.e. for integral action
Here mu is the controller gain and ti is the integral action time.
The modified program is here. Save it and try it out.
Click here to go to the model generator.
Section 4.5.2: Modelling a Vaporiser
Section 4.5.2: Modelling a Vaporiser
The following is a model of a single component vaporiser or boiler, developed using the approach described in section 4.2.

Section 4.5.3: Exercise: Modelling the vaporiser with a control system
Section 4.5.3: Exercise: Modelling the vaporiser with a control syste
The vaporiser model has no control system and will eventually either boil dry or go off the boil.

First add a proportional controller to regulate the holdup in the boiler (i.e. the level) at a setpoint of 0.8 kmol.
Then add a proportional controller to regulate the pressure at 1.5 bara.
Finally make the pressure a controller a PI controller.
Produce graphs of the level and pressure against time. Email the worksheet with graphs to us as an attachment.
Click here to go to the model generator.

step number	time t (s)	holdup M (kg)
0	0	500
1	20	520
2	40	540
3	60	560
4	80	580
5	100	600

step	time (s)	x(Euler)	x(Exact)	abs error
0	0.00	100.0000	100.0000	0.0000
1	0.25	75.0000	77.8801	2.8801
2	0.50	56.2500	60.6531	4.4031
3	0.75	42.1875	47.2367	5.0492
4	1.00	31.6406	36.7879	5.1473
5	1.25	23.7305	28.6505	4.9200
6	1.50	17.7979	22.3130	4.5151
7	1.75	13.3484	17.3774	4.0290
8	2.00	10.0113	13.5335	3.5222

step	time	c Euler	c exact	abs error
0	0.000	0.000000	0.000000	0.000000
1	35.355	0.014142	0.013188	0.000954
2	70.710	0.026284	0.024636	0.001648
3	106.065	0.036709	0.034575	0.002134
4	141.420	0.045660	0.043203	0.002457
5	176.775	0.053344	0.050693	0.002652
6	212.130	0.059942	0.057195	0.002747
7	247.485	0.065607	0.062840	0.002767
8	282.840	0.070471	0.067741	0.002731
9	318.195	0.074647	0.071995	0.002652

step	time	c	step size	abs error
0	0.000	0.000000
1	35.355	0.014142	35.355000	0.000954
2	70.710	0.026284	35.355000	0.001648
3	106.065	0.036709	35.355339	0.002134
4	144.222	0.046369	38.156248	0.002533
5	186.208	0.055376	41.986497	0.002857
6	231.751	0.063505	45.542669	0.003079
7	281.040	0.070700	49.289334	0.003193
8	335.044	0.077030	54.003740	0.003210
9	394.952	0.082534	59.907804	0.003135
10	462.030	0.087220	67.078139	0.002974
11	537.973	0.091102	75.943619	0.002729
12	625.322	0.094211	87.348124	0.002409
13	727.867	0.096586	102.545520	0.002025
14	851.443	0.098273	123.575870	0.001592
15	1000.000	0.099299	148.557098	0.001131

step	time (s)	c
0	0.0	0.099000
1	900.0	0.102600
2	1800.0	0.093240
3	2700.0	0.117576
4	3600.0	0.054302
5	4500.0	0.218814
6	5400.0	-0.208916
7	6300.0	0.903181
8	7200.0	-1.988271

step	time (s)	Implicit Euler	exact	error
0	0.0	0.099000	0.099000	0.000000
1	900.0	0.099783	0.099973	0.000190
2	1800.0	0.099953	0.099999	0.000047
3	2700.0	0.099990	0.100000	0.000010
4	3600.0	0.099998	0.100000	0.000002
5	4500.0	0.100000	0.100000	0.000000
6	5400.0	0.100000	0.100000	0.000000
7	6300.0	0.100000	0.100000	0.000000
8	7200.0	0.100000	0.100000	0.000000

step	time	M	iterations
0	0.0	100.000	0
1	50.0	169.719	3
2	100.0	221.043	3
3	150.0	259.846	3
4	200.0	289.654	3
5	250.0	312.794	3
6	300.0	330.890	3
7	350.0	345.117	3
8	400.0	356.346	3
9	450.0	365.235	3
10	500.0	372.280	2

step	time	x(1)	x(2)
0	0.0	0.00000	0.00000
1	1.0	0.01000	0.00000
2	2.0	0.01840	0.00001
3	3.0	0.02546	0.00002
4	4.0	0.03138	0.00004
5	5.0	0.03637	0.00007

10	10.0	0.05160	0.00024
20	20.0	0.06075	0.00069
40	40.0	0.06295	0.00167

100	100.0	0.06410	0.00457
200	200.0	0.06585	0.00922
300	300.0	0.06751	0.01364

800	800.0	0.07469	0.03271
900	900.0	0.07592	0.03599

step	time	Holdup	Outlet flow	% control signal
0	0.00	100.0000	0.0000	0.0000
1	10.00	140.0000	0.0000	0.0000
2	20.00	178.6650	0.1954	1.4660
3	30.00	208.2578	1.5985	13.3031
4	40.00	228.7844	2.2156	21.5138
5	50.00	244.8290	2.5444	27.9316
6	60.00	258.2513	2.7553	33.3005
7	70.00	269.8923	2.9068	37.9569
8	80.00	280.2064	3.0238	42.0825
9	90.00	289.4720	3.1183	45.7888
10	100.00	297.8756	3.1972	49.1502
11	110.00	305.5505	3.2647	52.2202
12	120.00	312.5971	3.3235	55.0388
13	130.00	319.0932	3.3752	57.6373
14	140.00	325.1017	3.4213	60.0407
15	150.00	330.6740	3.4627	62.2696
16	160.00	335.8534	3.5001	64.3414
17	170.00	340.6767	3.5341	66.2707
18	180.00	345.1755	3.5651	68.0702
19	190.00	339.3776	3.5250	65.7510
20	200.00	324.6830	3.4181	59.8732
21	210.00	311.0609	3.3109	54.4243
22	220.00	298.5167	3.2030	49.4067
23	230.00	287.0576	3.0945	44.8230
24	240.00	276.6902	2.9854	40.6761
..	......	........	......	.......
27	270.00	252.1062	2.6644	30.8425
30	300.00	236.6112	2.3888	24.6445
33	330.00	228.1177	2.1995	21.2471
36	360.00	223.9512	2.0938	19.5805
39	390.00	222.0400	2.0421	18.8160
42	420.00	221.1926	2.0185	18.4771
45	450.00	220.8229	2.0080	18.3292
48	480.00	220.6627	2.0035	18.2651
54	540.00	220.5636	2.0006	18.2255
60	600.00	220.5452	2.0001	18.2181
66	660.00	220.5418	2.0000	18.2167
72	720.00	220.5411	2.0000	18.2165

$\frac{dM}{dt} = F_1 - F_2$	ODE: material balance
$F_1 - \theta(t) = 0$	specify feed flow F₁ as function of time
$F_2 - \sqrt{\frac{M}{1/k^2 + 1/k_v^2} } = 0$	outflow - pressure drop equation
k - (0.2) = 0	specify pipe flow constant
$k_v - \phi k_{v,max} = 0$	valve characteristic equation
k_v,max - (1.0) = 0	specify fully-open valve flow constant
$g_{\phi}(\phi,y) = 0$	valve opening as function of controller output
y - K e = 0	proportional controller algorithm
$e - 100 \frac{(F_s - F_2)}{F_{range}} = 0$	controller error definition
K - (1.2) = 0	specify proportional gain
F_s - (2) = 0	specify flow set point
F_range - (5) = 0	specify flow range

M	tank holdup
F₁	feed flow
F₂	outlet flow
k	pipe flow constant
k_v	valve characteristic
$\phi$	valve open fraction
k_v,max	valve flow constant (fully open)
y	controller output
K	controller proportional gain
e	controller error signal
F_s	flow set point
F_range	flow range for controller

y₁	1.0000	1.2000	1.2710	1.2720	= y₁(0)
y₂	1.0000	0.8000	0.7860	0.7862	= y₂(0)

Equation name	Equation		Variables (units)
Balances
Material balance	$\frac{dM}{dt} = F_1 - F_2$	M	holdup (kg)
		F₁	feed flow (kg/s)
		F₂	outlet flow (kg/s)
Energy balance	$\frac{dU}{dt} = h_1 F_1 - h_2 F_2$	U	internal energy of contents (kJ)
(note 1)		h₁	feed specific enthalpy (kJ/kg)
		h₂	outlet specific enthalpy (kJ/kg)
Tank Model
gravity flow model	$vF_2 - 0.005 \sqrt{L} = 0$	v	specific volume ( m³/kg)
		L	level of contents (m)
`well mixed' model	T - T₂ = 0	T	tank temperature (⁰C)
		T₂	outlet temperature (⁰C)
mass - volume rel'n	M v - V = 0	V	volume of contents (m³)
volume - level rel'n	V - 0.2 ((1 + 0.8 L)³ - 1) = 0
(note 2)
Thermo and	definitions
energy - enthalpy rel'n	u - h + pv = 0	u	specific internal energy (kJ/kg)
		h	specific enthalpy (tank) (kJ/kg)
		p	pressure (kPa) (note 3)
specific energy def'n	U - M u = 0
Phys props
enthalpy correl'n	h₁ - 4.1846 T₁ = 0	T₁	feed temperature (⁰C)
(feed)
enthalpy correl'n	h - 4.1846 T = 0
(tank)
enthalpy correl'n	h₂ - 4.1846 T₂ = 0
(outlet)
density correl'n	(998.23 - 0.255 (T - 20)) v - 1 = 0
Specifications
feed flow spec	$F_1 - \left(2 + \frac{0.02 t}{1 + 0.01 t}\right) = 0$
feed temp. spec	T₁ - 20 = 0
pressure spec	p - 100 = 0