Numeric al anal ysis of proces ses
1 / 26


  • Uploaded on

NUMERIC AL ANAL YSIS OF PROCES SES. NAP 4. N ONLINEAR MODELS of systems described by ordinary differential equations Initial problem (Runge Kutta, Adams) Time delay Strange attractors. Rudolf Žitný, Ústav procesní a zpracovatelské techniky ČVUT FS 2010. Model s ODE initial probl e m.

I am the owner, or an agent authorized to act on behalf of the owner, of the copyrighted work described.
Download Presentation

PowerPoint Slideshow about 'NUMERIC AL ANAL YSIS OF PROCES SES' - senona

An Image/Link below is provided (as is) to download presentation

Download Policy: Content on the Website is provided to you AS IS for your information and personal use and may not be sold / licensed / shared on other websites without getting consent from its author.While downloading, if for some reason you are not able to download a presentation, the publisher may have deleted the file from their server.

- - - - - - - - - - - - - - - - - - - - - - - - - - E N D - - - - - - - - - - - - - - - - - - - - - - - - - -
Presentation Transcript
Numeric al anal ysis of proces ses


NONLINEAR MODELS of systems described by ordinary differential equations

Initial problem (Runge Kutta, Adams)

Time delay

Strange attractors

Rudolf Žitný, Ústav procesní a zpracovatelské techniky ČVUT FS 2010

Model s ode initial probl e m
Models ODE initial problem


Systems which are described by a system of ordinary differential equations and their solutions are fully described by the initial state, are for example integral models of mass and enthalpy balances elementary units ("lumped parameter" or "compartment" models). The aim is to determine the evolution of state variables (concentration, temperature) at the time, so it is an evolutionary problems.

These elementary units may be specific devices (reactors, absorbers, columns, tanks, filters, ...) in "spread sheets" or zones in apparatuses in which the state can be characterized by mean temperature, and the concentration of components (ideal mixers, plug flow zone). An example is the description of a stirred vessel system of ideal mixers below and above a turbine stirrer or a replacement shell and tube heat exchanger by a systém of elementary exchangers between the baffles.

Each elementary unit correspond to ordinary differential equations describing the mass balance (one for each component of the differential equation describing the rate of change of concentration)

and a enthalpy balance (or total energy balance) describing temperature evolution in the department

volumetric flowrate from neigbouring compartments

volumetric hear source

wnthalpic flows from neighbouring compartments

Model s ode initial probl e m1
Models ODE initial problem


A system need not be described by ordinary differential equations with only first derivatives. Eg. the equations of vibration include also the second time derivative (acceleration)

But every such equation can be replaced by a system of equations with only first derivatives, e.g.

The initial problem is characterized by the fact that each computed variable corresponds to one differential equation of the first order and to one initial condition (value) corresponding to the initial state.

…or how many equations so many initial conditions corresponding to one value of independent variable.

Model s ode initial problem
Models ODE initial problem


backmixing model (dispersion, diffusiion)

Example of continuous system:










Ci, Ti





model of stagnant zone

recirkulation with transport delay

Differential equations for concentrations

Parametersdetermining distribution of flowrates (s-stagnant zone, r-recirculation, b-backmixing) must be specified, r-from flowmeter measurement, b-from correlations. Also the volume of individual compartments must be known. These parameters are most frequently identiofied from fit of predicted characteristics (for example RTD - residence time distribution) with experiment.As a result information about the size of dead zones, recirculation zones, shortcuts, intensity of dispersion etc. are obtained.

Model s ode fourier transform
Models ODE Fourier transform


The models can be solved by the Fourier transformation method (see previous lecture NAZP3).

It is not enough to apply only a convolution, because the model includes compartments with back mixing, stagnant zone ... However, you can come out of these differential equations, transformed by FFT. For the Fourier transform of derivation applies (assuming that the concentration at time zero is infinity)

Differential equations are transformed to algebraic equations, where instead of time t is frequency f.

This system of algebraic equations (linear!) is already fairly easy to solve, and the result will be the Fourier transform of concentration. Concentration courses in time are then get back by FFT. The solution is accurate and fast (using FFT), but is limited to linear models. The procedure is preferred when the coefficients of the system (r, b, s) are constant. Otherwise, it was necessary to compute the Fourier transform of the product of two functions (and it is not equal to the product of the transformation). In the general case with nonlinear terms is therefore better to stay in the time domain and seek solutions to the numerical integration of the ODE.

Model s ode laplace transform
Models ODE Laplace transform


Similar relationships hold also for Laplace transform: the same convolution, correlation and transformation of derivatives

Inverse Laplace transform is usually searched in dictionaries, for example

Ode initial problem numerical solutions
ODEinitial problem numerical solutions












One-step methods approximate solution at time tk+1=tk+ by estimating mean value of derivativein interval<tk,tk+1>

Euler metodsare easily programmable even in Excel

first order of accuracy means that the error is directly proportional to the step

when you write a Taylor expansion for one time step, the error is O(2), a second-order accuracy. Errors, however, accumulate, so 1/ steps will result in an error O().

Euler method explicit

order of accuracy (first)

stableonly for small integration steps 

Euler method implicit

order of accuracy (first)

stable, but iterations are necessary (because ck+1is unknown)

Ode initial probl e m numeric al solutions
ODEinitial problem numericalsolutions







One-step methods approximate solution at time tk+1=tk+ by estimating mean value of derivativein interval<tk,tk+1>

Runge Kuttamethods calculate mean value of derivative as a weighted average of derivatives f(t,c). The most often used version calculates four points (endpoints tk, tk+1 and two midpoints at tk+1/2)

This variant is implemented in MATLABu as a functionode45

Weight coefficients (1/6,1/3,1/3,1/6) andincrements c aredesigned so that the order of accuracy is the highest (4th order of accuracy in this case). Generally, the order of accuracy of RK methods is the same as the number of evaluated points.

Ode initial probl e m numeric al solutions1
ODEinitial problem numerical solutions










Multistep methods approximate solution at time tk+1on the basis of several previous values for example at times tk, tk-1 (two), or tk-2 (three-step method).

This is usually a combination of predictor (explicit) and corrector (default formula with higher precision and stable)

it is sufficient to evaluatejust one value

Both formulas are second order accuracy but corrector (implicit formula) has the absolute value of the error about 9-fold lower (implicit formulas tend to have an error less than explicit, even if the order errors is the same).

The disadvantage is that for the first two steps a single-step method, for example RK, must be used. On the other hand, number of computed derivatives is smaller.

Ode initial probl e m numeric al solutions2
ODEinitial problem numerical solutions


Dynamic adjustment of integration step

The simplest solution is that each step is performed twice, once with step  and then with increments /2 (twice). If the difference is greater than the specified tolerance, the step  shrinks. Using the two calculated values (for full and half step) the accuracy of the final value can be improved by so called Aitken extrapolation. MATLAB uses a dynamic choice of the time step automatically.

Stiff problem

When you solve a system of equations may be that the integration step is dictated by the equation with an extremely short time constant, and the stability of the solution then requires extremely short integration step (the calculation time is prohibitively increased). We say that the system is stiff. Some methods of integrating ODE this unpleasant impact suppresses (an example is Gearova method).

Transport delay (time shift)

Time delay was included in the previous example. In this case the solver should remember the whole history of calculated results and not only the last step. Even the values for negative times (before the initial time) have to be defined.

Ode initial probl e m numeric al solutions3
ODE initial problem numerical solutions


Stability analysis

Consider a differential equation describing e.g. the rate of a chemical reaction of the first order

stabilitycondition restriction of time step (coef. of gain less than 1)

Euler explicit integration formula of the first order

Euler implicit integration formulaof the first order

gain less than 1 for arbitrary 

Ode examples in matlab
ODEexamples in MATLAB


ode45, ode23Runge-Kutta (one step method with different order of accuracy)

ode113Adams (multistep method)

ode15sStiff (Gear multistep methodwith variable accuracy order-suitable for stiff problems)

dde23transport delay

Rtd s eries of mixers 1 4
RTD series of mixers 1/4




c1 c2 c3 c4


An example has been solved in the previous lectures using FFT. Series 4 ideally mixed vessels (container volume, flow rate Q). The aim is to find the time courses of concentrations c1,…,c4for arbitrary initial conditions and for arbitrary course of concentration at inlet x(t).

Mass balances

tm=2 odpovídá poměru V/Q=0.5 (to je střední doba prodlení jedné nádoby)

As an initial concentration x(t) is applied impulse response of a set of two mixers (the same volume of tanks) and zero initial concentration in all tanks. The reason for such an assignment is that there is an analytical solution that can be compared with the results of numerical solution of the ODE. Stimulus function x(t) is the impulse response of E (t) for N=2 and tm=1; outletconcentration y(t)=c4(t) shou;ld be impulse response for N=6 and tm=3

RTD series of mixers2/4


MATLAB provides for the solution of ODE variety of methods (ode45, ode113, ode15s, ...) which can automatically adjust the integration step so as to achieve the required accuracy. Procedure is always the same, first define a function whose output is the vector of first derivatives for a given value of integration variable (time) and a solution vector.

In our case it is necessary to define a function whose result is a column vector (4x1) of right-hand sides of 4 differential equations. It will be assumed that the stimulus function x(t) is the impulse response of nx series mixers, with a mean timetmx. The parameter tm is the mean time of the modeled system of 4 mixers.

function dy = serie4(t,y,tm,tmx,nx)

dy = zeros(4,1);







It would probably be smarter if the stimulus function x(t) was directly the next formal parameter of the function serie4. I have a little problem with that in Matlab, I need advice.

The Efun function can be defined as a sub-function in the m-file serie4.m, but maybe it's better to define it as a separate m-file efun.m

function y=efun(t,tm,n)


RTD series of mixers3/4


Integration ODE by Runge Kutta method

Sol = ode45(function describing right side, time range, initial conditions)

in our case for example

sol=ode45(@(t,y)serie4(t,y,2,1,1),[0 10],[0 0 0 0])

Vector of initial conditions (zero concentrations)

solver: 'ode45'

extdata: [1x1 struct]

x: [1x29 double]

y: [4x29 double]

stats: [1x1 struct]

idata: [1x1 struct].

Integration from t=0, to t=10. It is not possible to specify a constant integration step in Matlab

Using a function as an actual parameter looks a bit strange in Matlab. It is not possibleto write simply the name of this function (serie4) because ode45 expects only two parameters (t, y) and the function serie4 has 5 parameters. Therefore, the so-called anonymous [email protected](t,y)expression.must be used.

Result is the structuresolcontaining vector of time steps sol.x (in this specific case 29 steps was necessary for prescribed accuracy) and a corresponding matrix sol.y of calculated concentrations (four rows for four concentrations)

RTD series of mixers4/4


Numerical solution by Runge Kutta (ode45) can be compared with the analytical solution E(t) for N=6 and tm=3

This is demonstration how to define a one-line function (without using a m-file)

[email protected](t,tm,n) n^n.*t.^(n-1)./(factorial(n-1)*tm^n).*exp(-n.*t/tm)



Standard function linspace generates vector of 1024 equidistant time steps


these points are results ofode45with automatically adjusted integration time steps

functiondeval interpolatesthe nonequidistant results sol.y according to the prescribed time scale t

RTDtime delay (transport delay)


In the case that there are time delays (long connecting pipelines, recycles) it is not possible to use the standard integration functions (ode23,ode45,…), but a different „family“ of MATLAB functions which is capable to save previous results (with prescribed vector of time delays1 2 …)

RTDdelay in recycle 1/2








Example: continuous system formed by two mixed tanks with time delay in recycle. It is quite common: connecting pipelines are usually modeled by a time delay.

The aim is to find time courses of concentrations c1(t), c2(t), for given flowrate (Q), recirculation ratio (r) and volume of tanks (V1,V2), for prescribed time course of inlet concentration x(t). The system is described by the mass balances:

time delay

RTDdelay in recycle 2/2


Special ODE solvers must be used in MATLAB (dde23fora constant time delay, andddesdfor variabl;e time delay)

sol=dde23(@derivatives, [1, 2…], history for t<tmin, [tmin,tmax])

function dcdt=dclag(t,c,z)





if t<.2









q-flowrate,v1,v2 volumes, r-recycle ratio, impulse stimulus function x(t) with duration 0.2 seconds

z is the matrix of solutions in times t-k. Column index (1) is the index in vector of delays

zero history c1=c2=0 for t<0



vector of delays1=1 (there is only one delay)



The initial problem isalso description of motion of particles under action of forces. Examples are trajectory of droplets in the spray drying chamber, motion of particles in a fluidized bed, burning coal particles in the combustion chamber, particles in the mixer, a settling tank, etc. The equations of motion are ODE of the type:


buoyant, centrifugal, electrical forces

drag for in fluid moving with velocity v

Trajectories of particles


In the last lecture of the course we will pay attention to the type of modern methods of DEM (Discrete Element Method), which simulates the dynamic behavior of particulate systems such as hoppers, silos, mills, mixers, separators, solid and fluid bed. Calculations are based on the solution of ordinary differential equations of motion and the consideration of interactions of individual particles.

Verlet integration of dynamic equations

Theory of chaos strange attractor
Theory of chaos strange attractor


Trajectories of particles calculated alsothe meteorologist E.Lorenz (1963) who tried to model the natural convection in atmosphere, the layer of air that is heated from below and cooled from above.

After drastic simplificationshe obtained the three ordinary differential equations for the coordinates x, y, z of a particle in the atmosphere

Rayleigh-Bénárdova instability

The model has only 3 parameters, Rayleigh number Ra, Prandtl number Pr and the parameter b is a slenderness ratio of rotating cells.

These 3 equations can be solve easily in MATLAB see next slide…

Theory of chaos strange attractor


trajektory pf particle fopr Ra=10, Pr=28, b=8/3

function dy = chaos(t,y)

dy = zeros(3,1);

dy(1) = 10*(y(2)-y(1));

dy(2) = -y(1)*y(3)+28*y(1)-y(2);

dy(3) = y(1)*y(2)-8/3*y(3);

sol=ode45(@chaos,[0 30],[.10001 .1 .1])





initial coordinates x,y,z

integration up to 30 s

plot graph y-z fro automatically calculated integration steps

Theory of chaos strange attractor


Numerically calculated trajectories behave erratically for Pr> 1 (they do not converge to a trivial steady steate solution x = y = z = 0), but they create something that looks like a counter-rotating swirls of final dimensions, and these seemingly random trajectories "jump" from the left to the right wing, and never intersect. The trajectories are extremely sensitive to initial conditions and to any disturbances (including approximation errors of numerical integration).

This limiting state is called „strange attractor“

parametric trajectoryof particle


Such behavior (deterministic chaos) appears at nonlinear systemswhen a stability limit is exceeded (eg Prantl, Reynolds or Rayleigh numbers). Deterministic chaos is typical for turbulent flow.




ODE Initial problem

numerical solutions, stability

Exam - remember


Fourier transform of derivatives

What is it initial problem

What is it order of accuracy O(n)

Explicit and implicit Euler method

Principles of RK and multistep methods (predictor, corrector)

Principle of stability analysis

Principle of integration step optimisation