Velten K. Mathematical Modeling and Simulation: Introduction for Scientists and Engineers

Подождите немного. Документ загружается.

3.2 Introductory Examples 125

rectangle symbolizing the alarm clock, T

outside that rectangle. But although we

already know the system S and the question Q of the mathematical model to be

developed, Model A is still not a complete description of a mathematical model (S,

Q, M) since mathematical statements M that could be used to compute the state

variables are missing.

3.2.3.3 Setting Up the Equations

In this case, it is better to improve Model A a little bit before going into a for-

mulation of the mathematical statements, M.BasedonModelA,theonlything

that the sensor (represented by T

) ‘‘sees’’ is the ambient air temperature. But if

this were true, then the initial decrease in the temperature data in Figure 3.2b

would be hard to explain. If the sensor ‘‘sees’’ only the ambient air temperature,

then its temperature should increase as soon as we enter the warm room at

time t = 0. The sensor obviously somehow memorizes the temperatures of the

near past, and this temperature memory must be included into our alarm clock

model if we want to reproduce the data of Figure 3.2b. Now there are several

possibilities how this temperature memory could be physically realized within

the alarm clock. First of all, the temperature memory could be a consequence of

the temperature sensor’s internal construction. As a ﬁrst, simple idea one might

hypothesize that the temperature sensor always ‘‘sees’’ an old ambient temper-

ature T

(t − t

lag

) instead of the actual ambient temperature T

(t). If this were

true, the above phenomenological model for temperature adaption, Equation 3.1,

could be used as follows. First, let us write down the ambient temperature T

for

this case

(t) =



t < t

lag

t ≥ t

lag

(3.5)

Here, T

is the ambient temperature before t = 0, that is, before we enter the

warm room. T

is the ambient temperature in the warm room. Since we assume

that the temperature sensor always sees the temperature at time t − t

lag

instead of

the actual temperature at time t, Equation 3.5 describes the ambient temperature

as seen by the temperature sensor. In Equation 3.1, T

(t) corresponds to the body

temperature, T

. This means that for t < t

lag

we have

(t) = T

− (T

− T

) · e

−r · t

(3.6)

and, for t ≥ t

lag

(t) = T

− (T

− T

lag

)) · e

−r · (t−t

lag

)

(3.7)

The parameters in the last two equations have the same interpretations as above

in Equation 3.1. Note that T

lag

) has been used as the initial temperature in

Equation3.7sinceweareshiftingfromT

to T

at time t

lag

,andT

lag

)isthe

actual temperature at that time. Note also that t − t

lag

appears in the exponent of

126 3 Mechanistic Models I: ODEs

Equation3.7tomakesurethatwehaveT

lag

) = T

lag

). The overall model can

now be written as

T(t) =



(t) t < t

lag

(t) t ≥ t

lag

(3.8)

3.2.3.4 Comparing Model and Data

Some of the parameters of the last equations can be estimated based on Figure 3.2b

and the corresponding data in

Room.csv:

≈ 18.5 (3.9)

lag

≈ 2.5 (3.10)

≈ 21 (3.11)

In principle, the remaining parameters (T

and r) can now either be determined

by heuristic arguments as was done above in the context of Equation 3.1, or by

nonlinear regression methods as described in Section 2.4. However, before this

is done, it is usually efﬁcient to see if reasonable results can be achieved using

hand-tuned parameters. In this case, a hand tuning of the remaining parameters

and r shows that no satisfactory matching between Equation 3.8 and the data can

be achieved, and thus any further effort (e.g. nonlinear regression) would be wasted.

Figure 3.4 shows a comparison of Equation 3.8 with the data of Figure 3.2b based

on the hand-ﬁtted values T

= 16.7andr = 0.09. The ﬁgure has been produced

using the Maxima code

RoomExp.mac and the data room.csv in the book software.

Looking at

RoomExp.mac, you will note that the if...then command is used to

implement Equation 3.8 (see Maxima’s help pages for more information on this

and other ‘‘conditional execution’’ commands).

RoomExp.mac computes a coefﬁcient of determination R

= 92.7%, reﬂecting

the fact that data and model are relatively close together. Nevertheless, the result

20.5

19.5

Temperature (°C)

18.5

t (min)

10 0

18.2

18.4

18.6

18.8

19.2

19.4

19.6

19.8

20.2

2 4 6 8 10 2 1420 30 40 50 60

(a)

t (min)

(b)

Fig. 3.4 (a) Comparison of Equation 3.8 (line) with the data

of Figure 3.2b (triangles) using T

= 16.7andr = 0.09.

(b) Same picture on a different scale, showing a dissimilarity

between model and data.

3.2 Introductory Examples 127

is unsatisfactory since the qualitative pattern of the model curve derived from

Equation 3.8 differs from the data. This is best seen if model and data are plotted

for t < 14 as shown in Figure 3.4b. As can be seen, there is a sharp corner in

the model curve, which is not present in the data. Also, the model curve is bent

upward for t > 3 while the data points are bent downward there. Although the

coefﬁcient of determination is relatively high and although we might hence be able

to compute reasonable temperature predictions based on this model, the qualitative

dissimilarity of the model and the data indicates that the modeling approach based

on Equation 3.8 is wrong. As it was mentioned in Section 1.2.2, the qualitative

coincidence of a model with its data is an important criterion in the validation of

models.

3.2.3.5 Validation Fails – What Now?

We have to reject our ﬁrst idea of how temperature memory could be included

into the model. Admittedly, it was a very simple idea to assume that the sensor

sees ‘‘old’’ temperatures T

(t − t

lag

)shiftedbyaconstanttimet

lag

.Oneofthe

reasons why this idea was worked out here is the simple fact that it led us to a

nice example of a model rejected due to its qualitative dissimilarity with the data.

Equation 3.8 also is a nice example showing that one cannot always distinguish in a

strict sense between phenomenological and mechanistic models. On the one hand,

it is based on the phenomenological model of temperature adaption, Equation 3.1.

On the other hand, Equation 3.1 has been used here in a modiﬁed form based on

our mechanistic considerations regarding the temperature memory of the system.

As was already mentioned in Section 1.5, models of this kind lying somewhere

between phenomenological and mechanistic models are also called semiempirical

or gray-box models.

Before going on, let us spend a few thoughts on what we did so far in terms of the

modeling and simulation scheme (Note 1.2.3). Basically, our systems analysis above

led us to the conclusion that our model needs some kind of temperature memory.

Equation 3.8 corresponds to the modeling step of Note 1.2.3, the simulation and

validation steps correspond to Figure 3.4. After the validation of the model failed,

we are now back in the systems analysis step. Principally, we could go now into a more

detailed study of the internal mechanics of the temperature sensor. We could, for

example, read technical descriptions of the sensor, hoping that this might lead us

on the right path. But this would probably require a considerable effort and might

result into unnecessarily sophisticated models. Before going into a more detailed

modeling of the sensor, it is better to ask if there are other, simple hypotheses that

could be used to explain the temperature memory of the system.

Remember that Note 1.2.2 says that the simplest model explaining the data is

the best model. If we ﬁnd such a simple hypothesis explaining the data, then

it is the best model of our data. This holds true even if we do not know that

this hypothesis is wrong, and even if the data could be correctly explained only

based on the temperature sensor’s internal mechanics. If both the models – the

wrong model based on the simple hypothesis and a more complex model based

on the temperature sensor’s internal mechanics – explain the data equally and

128 3 Mechanistic Models I: ODEs

indistinguishably well, then we should, for all practical purposes, choose the simple

model – at least based on the data, in the absence of any other indications that it is

awrongmodel.

3.2.3.6 A Different Way to Explain the Temperature Memory

Fortunately, it is easy to ﬁnd another, simple hypothesis explaining the temperature

memory of the system. In contrast to our ﬁrst hypothesis above (‘‘temperature

sensor sees old temperatures’’), let us now assume that the qualitative difference

between the data in the body temperature and alarm clock examples (Figures 3.1a

and 3.2b) is not a consequence of differences in the internal mechanics of the

temperature sensors, but let us assume that the temperature sensors used in both

the examples work largely the same way. If this is true, then the difference between

the data in the body temperature and alarm clock examples must be related to

differences in the construction of the clinical thermometer and the alarm clock

as a whole, namely to differences in their construction related to temperature

measurements. There is indeed an obvious difference of this kind: when the

clinical thermometer is used, the temperature sensor is in direct contact with the

body temperature that is to be measured. In the alarm clock, on the other hand,

the temperature sensor sits somewhere inside, not in direct contact with the

ambient temperature that is to be measured. This leads to the following.

Hypothesis:

The temperature memory of the alarm clock is physically realized in terms of the

temperature of its immediate surroundings inside the alarm clock, for example,

as the air temperature inside the alarm clock or as the temperature of internal

parts of the alarm clock immediately adjacent to the temperature sensor.

To formulate this idea in mathematical terms, we need one or more state

variable(s) expressing internal air temperature or the temperatures of internal parts

immediately adjacent to the temperature sensor. Now it was emphasized several

times that we should start with the simplest approaches. The simplest thing that

onecandohereistouseaneffective internal temperature T

, which can be thought of

as some combination of internal air temperature and the temperature of relevant

internal parts of the alarm clock. A more detailed speciﬁcation of T

would require

a detailed investigation of the alarm clock’s internal construction, and of the way in

which internal air and internal parts’ temperatures affect the temperature sensor.

This investigation would be expensive in terms of time and resources, and it

would hardly improve the results, which is achieved below based on T

as a largely

unspeciﬁed ‘‘black box quantity’’.

Note 3.2.4 (Effective quantities) Effective quantities expressing the cumulative

effects of several processes (such as T

) are often used to achieve simple model

formulations.

3.2 Introductory Examples 129

Introducing T

as a new state variable in Model A, we obtain Model B as an

improved conceptual model of the alarm clock (Figure 3.3). As a next step, we now

need mathematical statements (the M of the mathematical model (S, Q, M)tobe

developed) that can be used to compute the state variables. Since ODEs are required

here, we will go on with this example at the appropriate place below (Section 3.4.2).

It will turn out that Model B explains the data in Figure 3.2b very well.

3.2.3.7 Limitations of the Model

Remember that as a mechanistic modeler you are a ‘‘systems archaeologist’’,

uncovering the internal system mechanics from data similar to an archaeol-

ogist who derives subsurface structures from ground-penetrating radar data.

Precisely in this ‘‘archaeological’’ way, Figure 3.3 was derived from the data

in Figure 3.2b by our considerations above. Our starting point was given by the

data in Figure 3.2b with the puzzling initial decrease in the temperatures, an

effect that obviously ‘‘wants to tell us something’’ about internal system me-

chanics. Model B now represents a hypothesis of what the data in Figure 3.2b

might tell us about the internal mechanics of the alarm clock during temperature

measurement.

Note that Model B is a really brutal simpliﬁcation of the alarm clock as a real

system (Figure 3.2a), similar to our brutal simpliﬁcation of the system ‘‘car’’

in Section 1.1. In terms of Model B, nothing remains of the initial complexity

of the alarm clock except for the three temperatures T

, T

,andT

.Itmust

be emphasized that Model B represents an hypothesis about what is going on

inside the alarm clock during temperature measurement. It may be a right or

wrong hypothesis. The only thing we can say with certainty is that Model B

probably represents the simplest hypothesis explaining the data in Figure 3.2b.

Model B may fail to explain more sophisticated temperature data produced with

the alarm clock, and such more sophisticated data might force us to go into a more

detailed consideration of the alarm clock’s internals; for example, into a detailed

investigation of the temperature sensor, or into a more detailed modeling of the

alarm clock’s internal temperatures, which Model B summarizes into one single

quantity T

As long as we are concerned with the data only in Figure 3.2b, we can be

content with the rather rough and unsharp picture of the alarm clock’s internal

mechanics provided by Model B. The uncertainty remaining when mechanistic

models such as Model B are developed from data is a principal problem that

cannot be avoided. A mechanistic model represents a hypothesis about the internal

mechanics of a system, and it is well known that it is, as a matter of principle,

impossible to prove a scientiﬁc hypothesis based on data [95]. Data can be used

to show that a model is wrong, but they can never be used to prove its validity.

From a practical point of view, this is not a problem since we can be content with

a model as long as it explains the available data and can be used to solve our

problems.

130 3 Mechanistic Models I: ODEs

3.3

General Idea of ODE’s

3.3.1

Intrinsic Meaning of π

Simplicity is a characteristic of good mathematical models (Note 1.2.2), but also

a characteristic of good science in general (a fact that does not really surprise us

since good science must of course be based on good mathematical models ...).

A scientiﬁc result is likely to be well understood if it can be expressed in simple

terms. For example, if there is someone who wants to know about the meaning of

π, you could begin with a recitation like this

π = 3.14159265358979323846264338327950288419716939937510... (3.12)

If that someone is a clever guy, he might be able to do a similar recitation after

some time. But, of course, he would not understand the meaning of π ,evenifhe

could afford the time to recite every digit. Your next idea may be to use formulas

such as [96]

π = 4 ·

∞



n=0

(−1)

2n + 1

(3.13)

Probably, your listener would be happy that it is not so much effort to memorize

π this way. But, of course, he still would not understand. He would not understand

until you say a simple thing: ‘‘π is the surface area of a circle with radius 1.’’ He

would understand it this way because π is treated in the right setting here, namely

in terms of geometry. π can be treated in terms of algebra as above, but its intrinsic

meaning belongs to the realms of geometry, not algebra. This is reﬂected by the

fact that you have to tell comparatively long and complicated stories about π in

terms of algebra (Equations 3.12 and 3.13 above), while much more can be said by

a very simple statement in terms of geometry. The example may seem trivial, but

its practical importance can hardly be underestimated. To understand things, we

should always try to ﬁnd a natural setting for our objects of interest.

3.3.2

Solves an ODE

With this idea in mind, let us reconsider the body temperature example from

the last section. Remember that the data in Figure 3.1a were described using the

function

T(t) = T

− (T

− T

) ·e

−r · t

(3.14)

The main ingredient of this function is the exponential function. Now, sim-

ilar to above, let us ask a naive question: What is the exponential function? If

3.3 General Idea of ODE’s 131

you never have thought about this, you may start with similar answers as above.

For example, the exponential function is the power function e

having the Euler

number

e = 2.718281828459045235360287471352662497757247093699959... (3.15)

as its base. Or, you may use a formula similar to Equation 3.13 [17]:

∞



n=0

(3.16)

Undoubtedly, answers of this kind would prove that you are knowing a lot more

than those students who believe that ‘‘the exponential function is the e

key on

my pocket calculator’’. But you never would be able to understand the intrinsic

meaning of the exponential function this way. As above, algebra just is not the

rightsettingifyouwanttounderstande

The right setting for an understanding of the exponential function are ODEs,

just as geometry is the right setting to understand π. This is related with a fact

that most readers may remember from calculus: The derivative of e

is e

,thatis,e

anditsderivativecoincide.IntermsofODEs,thisisexpressedasfollows:e

solves

the ODE



= y (3.17)

with initial condition

y(0) = 1 (3.18)

3.3.3

Inﬁnitely Many Degrees of Freedom

Let us say a few words on the meaning of the last two equations (precise deﬁnitions

are given in Section 3.5). At a ﬁrst glance, Equation 3.17 looks similar to many

other algebraic equations that the reader will already have encountered during his

mathematical education. There is an unknown quantity y in this equation, and we

have to ﬁnd some particular value for y such that Equation 3.17 is satisﬁed. There

is, however, an important difference between Equation 3.17 and standard algebraic

equations: the unknown y in Equation 3.17 is a function. To solve Equation 3.17, y

must be replaced by a function y (t) which satisﬁes Equations 3.18 and 3.17 for every

t ∈ R.Notethaty



is the derivative of y with respect to the independent variable,

t in this case. Any of the common denotations of derivatives can be used when

writing down differential equations. In particular,

y is frequently used to denote a

time derivative.

The fact that functions serve as unknowns of differential equations makes their

solution a really challenging task. A simple analogy may explain this. Consider, for

132 3 Mechanistic Models I: ODEs

example, a simple linear equation such as

3x − 4 = 0 (3.19)

which everyone of us can easily solve for x. Now you know that the number of

unknowns and equations can be increased arbitrarily. For example,

2x − 8y = 4

7x − 5y = 2

(3.20)

is a system of two linear equations in the two unknowns x and y,andthis

generalizes to a system of n linear equations in n unknowns as follows:

+ a

+···+a

= b

+ a

+···+a

= b

+ a

+···+a

= b

(3.21)

To solve the last system of equations, the n unknowns x

, x

, ..., x

must be

determined. Although there are efﬁcient methods to solve equations of this kind

(e.g. Maxima’s

solve command could be used, see also [17]), it is obvious that the

computational effort increases beyond all bounds as n →∞. Now remember that

the unknown of a differential equation is a function. At a ﬁrst glance, it may seem

that we, thus, have ‘‘less unknowns’’ in an equation such as (3.17) compared to

Equation 3.21: y as the only unknown of Equation 3.17, and x

, x

, ..., x

as the n

unknowns of Equation 3.21. But now ask yourself how many numbers x

, x

, ..., x

you would need to describe a function y = f (x)onsomeinterval[a,b]. In fact, you

would need inﬁnitely many numbers, as you would have to recite all the values of

that function for those inﬁnitely many values of the independent variable x in the

interval [a,b]. This is related to the fact that a real function y = f (x)onsomeinterval

[a,b] is said to have inﬁnitely many degrees of freedom on its domain of deﬁnition. If we

solve a differential equation, we thus effectively solve an equation having inﬁnitely

many unknowns corresponding to the solution functions degrees of freedom. The

solution of differential equations is indeed one of the most challenging tasks in

mathematics, and a substantial amount of the scientiﬁc work in mathematics has

been devoted to differential equations since many years.

3.3.4

Intrinsic Meaning of the Exponential Function

It was said that to solve Equation 3.17, we need to replace the unknown y by

afunctionf (t)suchthatavalidequationresults.Choosingf (t) = e

and using



(t) = e

, Equation 3.17 turns into

= e

(3.22)

3.3 General Idea of ODE’s 133

which means that e

indeed solves Equation 3.17 in the sense described above.

Note that f (t) = c · e

would also solve Equation 3.17, where c ∈ R is some constant.

This kind of nonuniqueness is a general characteristic of differential equations,

usually overcome by imposing appropriate additional initial conditions or boundary

conditions.Equation3.18providesanexampleofaninitialcondition(boundary

conditions are discussed further below). It is not really hard to understand why this

condition is needed here to make the solution of Equation 3.17 unique. Remember

how we motivated the use of differential equations above (Note 3.1.1): they provide

us with a means to formulate equations in terms of the rates of change of the

quantities we are interested in. From this point of view, we can say that Equation

3.17 ﬁxes the rate of change for the quantity y in a way that it always equals the

current value of y. But it is of course obvious that we cannot derive the absolute

value of a quantity from its rate of change. Although your bank account may

have the same interest rates as the bank account of Bill Gates, there remains a

difference that can be made precise in terms of the concrete values of those bank

accounts at some particular time, which is exactly what is meant by an initial

condition.

Equation 3.18 is an initial condition for the ODE (Equation 3.17) since it

prescribes the initial value of y at time t = 0. f (t) = e

, our solution of Equation

3.17, obviously also satisﬁes Equation 3.18, and it thus solves the so-called initial

value problem comprising Equations 3.17 and 3.18. Thus, we can now answer the

above question ‘‘What is the exponential function?’’ as follows:

The exponential function is the real function f (t) which coincides everywhere

with its derivative and satisﬁes f (0) = 1.

This statement characterizes the exponential function in a simple and natural

way, similarly simple and natural as the characterization of π in Section 3.3.1.

Note that unique solvability of Equations 3.17 and 3.18 is assumed here – this is

further discussed in Section 3.5. There is indeed no simpler way to explain the

exponential function, and the complexity of algebraic statements such as Equations

3.15 and 3.16 emanates from that simple assertion above just as the algebraic

complexity of π emanates from the assertion that ‘‘π is the surface area of a circle

having radius 1’’. As an alternative, less mathematical formulation of the above

‘‘explanation of the exponential function’’, we could also say (neglecting the initial

condition)

The exponential function describes a process where the rate of change of a

quantity of interest always equals the actual value of that quantity.

This second formulation makes it understandable why the exponential function

appears so often in the description of technical processes. Technical processes

involve rates of changes of quantities of interest, and the exponential function

describes the simplest hypothesis that one can make regarding these rates of

134 3 Mechanistic Models I: ODEs

changes. This is the intrinsical meaning of the exponential function, and this

is why it is so popular. Note that the exponential function can also be used in

situations where the rate of change of a quantity of interest is a linear function of

the actual value of that quantity (Sections 3.4.1 and 3.7.2.1).

3.3.5

ODEs as a Function Generator

Ordinary differential equations also provide the right setting for an understanding

of many other transcendental functions. For example, the cosine function can be

characterized as the solution of



=−y (3.23)

with the initial conditions

y(0) = 1



(0) = 0

(3.24)

cos(t) satisﬁes Equation 3.23 in the sense described above since cos



(t) =−cos(t)

for t ∈ R. Equation 3.24 is also obviously valid since cos(0) = 1andcos



(0) =

−sin(0) = 0. Note that two initial conditions are required here, which is related

to the fact that Equation 3.23 is a so-called second-order equation since it involves

a second-order derivative. The order of a n ODE is always the order of its highest

derivative (precise deﬁnitions are given in Section 3.5).

Note 3.3.1 (Understanding transcendental functions) Ordinary differential

equations provide a natural framework for an understanding of the intrinsic

meaning of functions that are frequently used to describe technical processes

such as the exponential function and the trigonometric functions.

Beyond this, we will see that ODEs virtually serve as what might be called a

function generator in the sense that they can be used to produce functions that do not

belong to the classical zoo of functions such as the exponential and trigonometric

functions. An example is the function describing the time dependence of cell

biomass during wine fermentation (see Figure 3.20 in Section 3.10.2). This function

cannot by described by any of the classical functions, but it can be described using

asystemofODEsasitisexplainedbelow.

Note 3.3.2 (ODEs as function generator) ODEs can be used to generate func-

tions that do not belong to the classical zoo of functions such as the exponential

and trigonometric functions.