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

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

3.8 Numerical Solutions 185

1: dfn <-

2: function(t, y, p)

3: {

4: ### User: names of state variables

5: ##################################

6: T=y[1]

7: ### User: differential equation

8: ##################################

9: dTdt=r*(Tb-T)

10: ### User: list of derivatives

11: ##################################

12: list(c(dTdt))

13: }

(3.241)

As a whole, the code 3.241 deﬁnes a function dfn that gives the right-hand side

of the ODE and is used later in

ODEEx1.r to solve the ODE. In the general case,

this function assumes a system of ODEs of the form (compare Deﬁnition 3.5.3)



(t) = F(t, y(t), p) (3.242)

where y(t) = (y

(t), y

(t), ..., y

(t))

expresses the n dependent variables of the

ODE and p = (p

, p

, ..., p

)

is a vector of parameters on which the right-hand

side of the ODE may depend (n, s ∈ N). We do not use p here, but it is

important in the estimation of parameters of the ODE from experimental

data (Section 3.9). The arguments of

dfn in program 3.241 are in obvious

correspondence with the arguments of F(t, y(t), p) in Equation 3.242. The com-

ponents of

y and p are accessed as y[1], y[2], ... and p[1], p[2], ... respec-

tively, in R. Square brackets are used generally in R to access array compo-

nents [45].

In lines 4–6 of program 3.241, you set names for the dependent variables

that are then used in lines 7–9 to deﬁne the ODE. Lines 4–6 can be omitted if

you formulate the ODE using

y[1], y[2] etc. as the dependent variables. Note,

however, that your code will be better readable if you are using the original names

of your dependent variables instead of

y[1], y[2] etc. Line 9 deﬁnes the ODE

in obvious correspondence with Equation 3.231. Line 12 of program 3.241 is

the last statement within the function body of

dfn (all statements enclosed in

brackets between lines 3 and 13), and it is, thus, the function value produced by

dfn. dfn is expected to produce a list of all right-hand sides of the ODE system,

which is a mere

list(c(dTdt)) in this case and which would read, for example,

list(c(dTdt,dSdt)) if our ODE would involve a second state variable S (see the

examples below). Note that the concatenate operator

c(...) is used in R to build

vectors from its arguments [45].

186 3 Mechanistic Models I: ODEs

3.8.3.2 Deﬁning Model and Program Control Parameters

Now you know how to deﬁne an ODE in R – let us go on looking through ODEEx1.r.

The next thing you can do is the deﬁnition of a closed form solution,ifavailable,and

if you would like to compare the numerical result with a closed form solution. In

the body temperature example, we have Equation 3.230 as a closed form solution of

Equations 3.231 and 3.232. Program 3.243 shows the part of

ODEEx1.r that deﬁnes

the closed form solution, Equation 3.230. As you see, the closed form solution func-

tion is denoted as

AnSol in ODEEx1.r. AnSol corresponds to T(t) in Equation 3.230,

and in R it can be invoked in the form

AnSol(t). The right-hand side of Equation

3.230 is in immediate correspondence with the expression you see in line 3 of 3.243.

1: ### User: closed form solution (if available)

2: #############################################

3: AnSol=function(t) Tb-(Tb-T0)*exp(-r*t)

(3.243)

After this, there is a section Parameters of the model in ODEEx1.r where you

do exactly what this section title suggests. In the body temperature model, there

are three parameters to be deﬁned (r, T

,andT

) and looking into ODEEx1.r you

willseethatexactlythesamevaluesasintheMaxima program

FeverODE.mac are

used here.

In the Program control parameters section of

ODEEx1.r,youbeginwitha

selection of the plots you want to see by setting the appropriate logical variables.

After this, the start and end times

aa and bb are deﬁned, which limit the interval

in which you want to see the numerical solution. If applicable, you can set the

number of points of the analytical solution in the plot using

nAnSol,andyou

can deﬁne a data ﬁle in csv-format in the variable

DataFile (e.g. if you want

to compare the numerical solution with data). Using

nTime, you can prescribe

the number of intermediate points between

aa and bb where you want lsoda to

compute numerical approximations (this corresponds to n in the Euler method,

see Section 3.8.1.1). Basically,

nTime controls the smoothness of the numerical

solution curve that you see in the plots, since this curve is a polygon connecting

the points of the solution that have been computed by the numerical algorithm.

Note that

nTime assumes an equidistant distribution of the intermediate points

between

aa and bb. You can easily change this if you like by looking at the way in

which

nTime is used later in ODEEx1.r.Itisimportantthatyoudonotconfuse

the distance between numerical approximation points controlled by

nTime with

the stepsize h discussed above. h is automatically controlled by

lsoda so as to

satisfy given local error bounds (see below), while

nTime just controls the amount

of output you need.

3.8.3.3 Local Error Control in lsoda

After this you set the absolute and relative local error tolerances, atolDef and

rtolDef. As it was mentioned above, the fact that these error tolerances can be

prescribed using

lsoda is an important reason why you should better use lsoda

3.8 Numerical Solutions 187

instead of Maxima’s rk function in your everyday work, although the use of

Maxima is simpler as we have seen. Suppose we are in n dimensions, that is,

we have y(t) = (y

(t), y

(t), ..., y

(t))

in Equation 3.242. Assume that the solver

has just performed one step of its numerical algorithm leading to a numerical

approximation ˜y(t) = (

(t),

(t), ...,

(t))

at some time t. You may imagine one

step of the Euler method here that advances the numerical solution from a previous

time t − h to the current approximation at time t. Above it was demonstrated how

the error in one step of the Euler method can be estimated, see Equation 3.224.

This kind of per-step-error is called the local error of the numerical algorithm. It

must be distinguished from the global error of the numerical method, which is the

difference between the numerical approximation of the solution of the ODE and

the exact solution of the ODE, both considered over the entire interval of interest

(the difference being measured in some appropriate norm that we do not need to

discuss here).

Using methods appropriate for the numerical algorithms used in

lsoda,the

local error in one step of this algorithm can be estimated similar to Equation 3.224.

For each of the n components of ˜y(t),

lsoda computes such an estimated local

error, which we denote e

. Then, the local errors e

are required to satisfy

≤ rtol

·|

(t)|+atol

, i = 1 ...n (3.244)

Here, rtol

and atol

are the relative and absolute local error tolerances that you

specify in

ODEEx1.r using atolDef and rtolDef. For example, assume n = 2,

and assume you want to have atol

= 10

−2

and atol

= 10

−4

. Then, you can do this

by setting

1: # absolute tolerance

2: atolDef=c(1e-2,1e-4)

(3.245)

in the ‘‘program control parameters’’ section of ODEEx1.r. On the other hand,

1: # absolute tolerance

2: atolDef=1e-3

(3.246)

would give atol

= 10

−3

and atol

= 10

−3

.rtol

canbecontrolledinthesameway.

Basically,

lsoda uses Equation 3.244 to optimize the stepsize, that is, the stepsize

is chosen as large as possible such that Equation 3.244 is still satisﬁed.

3.8.3.4 Effect of the Local Error Tolerances

Various strategies of local error control can be realized by an appropriate choice of

atolDef and rtolDef. Consider one particular solution component i ∈

{

1 ...n

}

Setting rtol

= 0 would transform Equation 3.244 into

≤ atol

, i = 1 ...n (3.247)

188 3 Mechanistic Models I: ODEs

Choosing atol

= 10

−3

would then mean that a step of the numerical algorithm

inside

lsoda is accepted only if the absolute value of the estimated local error of

solution component i is below 10

−3

. Setting atol

= 0, on the other hand, would

transform Equation 3.244 into

≤ rtol

·|

(t)|, i = 1 ...n (3.248)

In this case, choosing rtol

= 10

−3

would mean that a step of the numerical

algorithm inside

lsoda is accepted only if the local error of solution component

i is smaller than 10

−3

when expressed relative to the size of the numerical

approximation of solution component i, |

|. Equations 3.247 and 3.248 correspond

to (purely) absolute error control and (purely) relative error control, respectively, while

Equation 3.244 with rtol

= 0andatol

= 0 expresses a mixed absolute and relative

error control. Such a mixed error control basically corresponds to a relative error test

when the solution component |

(t)|is much larger than atol

, while it approximates

an absolute error test when |

(t)| is smaller than atol

3.8.3.5 A Rule of Thumb to Set the Tolerances

This is the theory – but how should atol

and rtol

be chosen in a practical situation?

Naively, one might conjecture that it is a good idea to choose these tolerances as

small as possible. But this is a bad idea since Equation 3.244 would then force

the numerical algorithms inside

lsoda to use extremely small stepsizes, which

would increase the computational time beyond all bounds and, more seriously,

which might increase the global error in an uncontrollable way since much more

computational steps would be needed until the numerical approximation of the

solution of the ODE would be obtained.

If you want to know more about the various kinds of errors in numerical

algorithms, we recommend [111] for a nice discussion of ‘‘errors and uncertainties

in computations’’ or the comments on ‘‘error, accuracy and stability’’ in [92]. Two

important kinds of error discussed there are the roundoff error and the truncation

error or approximation error. Basically, the truncation or approximation error is

related to the fact that any numerical algorithm can only perform a ﬁnite number

of steps. Remember the discussion of the Euler method in Section 3.8.1.1, which

evaluates the right-hand side of the ODE only at a ﬁnite number of x values, and

which is based on formula (3.208), a truncated version of the originally inﬁnite

Taylor series. The roundoff error, on the other hand, arises because computers

store ﬂoating-point numbers using only a ﬁnite number of digits, which means,

for example, that a number such as 1/3 is stored as 0.3333 if we assume a computer

storing four decimal digits. If you then use 0.3333 instead of the exact 1/3 in your

computations, it is obvious that this and any other arithmetical operation on such

a computer may introduce a roundoff error in the order of 10

−4

. These errors add

up with the number of steps performed in your numerical algorithm, and this is

the reason why atol

and rtol

should be chosen with care (not too small). As a

guideline, we recommend to use the following ‘‘rule of thumb’’ from [112]:

3.8 Numerical Solutions 189

Note 3.8.2 (Rule of thumb to set the tolerances) If m is the number of signiﬁcant

digits required for solution component y

,setrtol

= 10

−m+1

.Setatol

to the value

at which |y

| is essentially insigniﬁcant.

Letusapplythistothebody temperature example. In this example, there is only

one dependent variable (T in Equation 3.231), and hence we have n = 1, which

means we can omit the index at atol

and rtol

and simply write atol and rtol instead.

The body temperature example refers to a clinical thermometer that works in a range

roughly between 30 and 42

◦

C, and which displays one digit behind the decimal

point. Therefore, it is clear that we need at least three signiﬁcant digits here. In

terms of the above ‘‘rule of thumb’’, this means m = 3andrtol= 10

−2

.Regarding

atol, values below 10

−1

can be considered as ‘‘essentially insigniﬁcant’’ since only

one digit behind the decimal point is displayed by the clinical thermometer. In

terms of the ‘‘rule of thumb’’ this means we should set atol = 10

−2

. These settings

of atol and rtol can be viewed as the largest values of these tolerances that may give

us the desired result with sufﬁcient accuracy.

However, you should remember that atol and rtol limit the local error of the

algorithm. As explained above, what really counts is the global error, that is, the ﬁnal

difference between the exact solution of the ODE and the numerical approximation.

In each step of the numerical algorithm, the local errors add up to the global error.

If we are lucky, the signs of the local errors produced in each step of the numerical

algorithms may be randomly distributed and thus cancel out, giving a global error

in a similar order of magnitude as the local errors. Fortunately, this is frequently the

case, that is, we may expect the global errors to be in a similar order of magnitude

as the local errors in many cases. But it may also happen that the local errors add

up to give a much larger global error. You get a fairly reliable idea about the size of

the global error by using the heuristic procedure explained above in our discussion

of Maxima’s

rk command (Note 3.8.1). There it was recommended to compute

numerical solutions for two stepsizes h and h/2, and then to compare the results.

The same can be done here by reducing the tolerances, for example, by an order

of magnitude. If the differences between the numerical solutions obtained in this

way are negligible, you have a good reason to assume a negligible global error, and

hence that your choice of the tolerances was acceptable.

3.8.3.6 The Call of lsoda

Remember that we are still discussing the R code ODEEx1.r in the book software.

The error tolerances

rtolDef and atolDef are the last settings in the ‘‘Program

control parameters’’ section of

ODEEx1.r. If you go through the rest of the code, you

will ﬁnd several other places where the user needs to do some smaller adjustments,

such as the speciﬁcation of column names in the data ﬁle, setting initial values

for

lsoda using appropriate variable names, and adjusting the plots. Those places

are indicated by commentaries beginning with ‘‘

# User:...’’, such as lines 1–2

in program 3.243. As a last point in our discussion of

ODEEx1.r, let us have a look

190 3 Mechanistic Models I: ODEs

at the call of lsoda within that code, since this will be important, for example, in

our discussion of parameter estimation in ODEs in Section 3.9:

1: out = as.data.frame(

2: lsoda(

3: c(T=T0)

4: ,seq(aa,bb,length=nTime)

5: ,dfn

6: ,c()

7: , rtol= rtolDef

8: , atol= atolDef

9: ))

(3.249)

Note that we have left out the user comments here for brevity. In line 1 of 3.249,

the output of

lsoda is stored as a ‘‘data frame’’ since this is the form needed

further below in

ODEEx1.r to produce the plots (see [45] for more on the data frame

concept). Line 3 corresponds to Equation 3.232 of the body temperature model, that

is, here you supply the initial values of the state variables. Line 4 deﬁnes a vector of

times where

lsoda is required to generate numerical approximations of the state

variable(s), that is, of the temperature T(t) in the case of the body temperature

model, Equations 3.231 and 3.232. The

seq command in line 4 generates an

array of

nTime equally spaced times between aa and bb. You can replace this

by any other vector of times that you might want to use. In Section 3.10.2, for

example, a code

Fermentation.r is discussed where a sophisticated, nonequally

spaced time vector is used in the call of

lsoda. Line 5 tells lsoda about the

function that generates the right-hand sides of the ODE, which is

dfn in our case

(Section 3.8.3.1). In line 6 of 3.249 you can add parameters required by

dfn,which

will be used in the parameter estimation problems treated in Section 3.9. In that

section, other optional parameters of

lsoda such as hmax are also discussed, which

canbesetbyaddingalinesuchas

,hmax=0.001

e.g. between lines 8 and 9 of 3.249.

3.8.3.7 Example Applications

If you now set PlotStateAna=TRUE and PlotStateData=TRUE in the ‘‘Pro-

gram control parameters’’ section of

ODEEx1.r and then execute the code e.g. using

the

source command as described in Appendix B, you get the result shown in

Figure 3.11. Figure 3.11a shows a comparison of the numerical solution of the body

temperature model (Equations 3.231 and 3.232) with the appropriate closed form

solution (Equation 3.230). The numerical and closed form solutions match perfectly

well, as it was achieved previously using Maxima’s

rk command (Figure 3.8).

Looking at

ODEEx1.r you will ﬁnd that rtolDef=1e-4 and atolDef=1e-4 have

been used for this result. Remember that the ‘‘rule of thumb’’ discussed above

led us to

rtolDef=1e-2 and atolDef=1e-2. The actual, smaller values of the

3.8 Numerical Solutions 191

Temperature (°C)

Numerical solution

Analytical solution

Numerical solution

Data

Temperature (°C)

0 5 10 15 20 0 5 10 15 20

t (min)

(a) (b)

t (min)

= 0.998

Fig. 3.11 (a) Line: numerical solution of Equations 3.231

and 3.232 computed with

ODEEx1.r.Circles:exactsolution

(Equation 3.230). (b) Line: numerical solution of Equations

3.231 and 3.232 computed with

ODEEx1.r.Circles:data

Fever.csv.

tolerances have been chosen based on an application of the heuristic procedure

discussed above, that is, by taking

rtolDef=1e-2 and atolDef=1e-2 as a starting

point, and then comparing the result with the result obtained using tolerance

values one order of magnitude smaller, and so on, until a further reduction of

the tolerances did not affect the result. Figure 3.11b shows a comparison of the

numerical solution of Equations 3.231 and 3.232 with the data in

Fever.csv

similar to Figure 3.1b above which was produced using the Maxima program

FeverExp.mac. Figure 3.11b reports an R

value of 0.998, which coincides with

the value computed by

FeverExp.mac.

As a second example of using R to solve ODEs, consider the alarm clock model

(Equations 3.234–3.237). To solve these equations based on R,weusetheprogram

ODEEx2.r in the book software. The general structure of this program is identical

with

ODEEx1.r as discussed above, the only difference being the fact that we are

concerned with a system of two ODEs here. In terms of the

dfn function, this is

expressed as follows (compare with 3.241):

1: dfn <-

2: -function(t, y, p)

3: {

4: ### User: names of state variables

5: ##################################

6: Ts=y[1]

7: Ti=y[2]

8: ### User: differential equation

9: ##################################

10: dTsdt=rsi*(Ti-Ts)

11: dTidt=ria*(Ta-Ti)

12: ### User: list of derivatives

(3.250)

192 3 Mechanistic Models I: ODEs

13: ##################################

14: list(c(dTsdt,dTidt))

15: }

As in 3.241, the code begins with a deﬁnition of the state variables in lines 6–7.

In this case, there are two state variables, which are denoted as

y[1] and y[2]

(y[3], y[4], ... in cases with more than two variables). In lines 6–7, these state

variables are named as Ts and Ti corresponding to T

and T

in Equations 3.234

and 3.235. Lines 10–11 deﬁne the ODEs (Equations 3.234 and 3.235). In line 14,

ﬁnally, the list of derivatives computed in lines 10–11 is returned as the result of

the

dfn function. Note that in line 12 of 3.241 this was a list of length 1 since only

one ODE was solved there. This list must always have as many elements as the

number of ODEs in the ODE system to be solved. As before, the deﬁnition of

dfn

is followed by the deﬁnition of the closed form solution (Equations 3.238 and 3.239

in this case). Since the parameters used in

ODEEx2.r satisfy r

= r

,Equations

3.238 and 3.239 lead us to

1: ### User: closed form solution (if available)

2: ##################################

3: AnSolTs=function(t) Ta+(Ts0-Ta)*exp(-rsi*t)

+(Ti0-Ta)*rsi*(exp(-rsi*t)

-exp(-ria*t))/(ria-rsi)

4: AnSolTi=function(t) Ta+(Ti0-Ta)*exp(-ria*t)

(3.251)

Compared with Equation 3.243 where there was only one closed form solution

AnSol, you see that two closed form solutions AnSolTs and AnSolTi are deﬁned

here, corresponding to Equations 3.238 and 3.239. In the ‘‘Parameters of the

model’’ section of

ODEEx2.r, the parameters are set to the same values as in the

Maxima code

RoomODE.mac discussed above. In the ‘‘Solution of the ODE’’ section

ODEEx2.r, you have to set the initial values as follows:

1: ### User: set initial values

2: ##################################

3: c(Ts=Ts0,Ti=Ti0)

(3.252)

In ODEEx1.r, the corresponding code is

1: ### User: set initial values

2: ##################################

3: c(T=T0)

(3.253)

So you see that you have to provide a list of initial values here with as many

elements as the number of ODEs in your ODE system. In the last part of

ODEEx2.r

where the plots are deﬁned, we need two plots here comparing a numerical solution

with a closed form solution, while only one such plot was needed in

ODEEx1.r.

Comparing

ODEEx1.r and ODEEx2.r, you will see that this was achieved just by

3.8 Numerical Solutions 193

20.5

19.5

Sensor temperature (°C)

Internal temperature (°C)

18.5

0 102030

t (min)

(a) (b)

t (min)

40 50 60 0 10

Numerical solution

Analytical solution

20 30 40 50 60

Fig. 3.12 (a) Line: numerical approximation of T

based on

Equations 3.234–3.237. Circles: exact solution (Equation

3.238). (b) Line: numerical approximation of T

based on

Equations 3.234–3.237. Circles: exact solution (Equation

3.239). Plots generated using

ODEEx2.r.

18.5

19.5

20.5

= 0.994

Numerical solution

Data

(

min

)

Sensor temperature (°C)

0 102030405060

Fig. 3.13 Line: Numerical approximation of T

based on

Equations 3.234–3.237, obtained using

ODEEx2.r.Circles:

data from

room.csv.

using two copies of the relevant part of ODEEx1.r, one of these copies plotting Ts

compared with AnSolTs and the other one plotting Ti compared with AnSolTi.

Figures 3.12 and 3.13 show the result produced by

ODEEx2.r. Compare Figure 3.12

with 3.10, which shows the same result as it was obtained using Maxima’s

function. As before, there is a perfect coincidence between the numerical solution

and the closed form solution. Figure 3.13 shows a comparison of the numerical

approximation with the data in

room.csv similar to Figure 3.5 (see the discussion

of these data in Section 3.2.3).

We have seen now how ODEs can be solved using Maxima

rk or R’s lsoda

commands. Although the procedure in Maxima is simpler, we have seen that

it is better to use R’s

lsoda in your everyday work since this is a much more

professional software that provides, for example, a local error control and a number

194 3 Mechanistic Models I: ODEs

of other options which you do not have when using Maxima.Ifyouwantto

apply R’s

lsoda to your own ODE systems, you can use the codes ODEEx1.r

and ODEEx2.r discussed in this section as a starting point, editing these codes

appropriately (for example, you will then have to insert your ODE into the

dfn

function as discussed above). All examples in the following sections are treated

based on R’s

lsoda command. Some of these examples use extended version of

ODEEx1.r and ODEEx2.r, using new options of R and lsoda. For example, in the

discussion of the wine fermentation model (Section 3.10.2), a maximum stepsize

hmax is prescribed in lsoda to make sure that external supplies of nitrogen are

correctly used in the computations.

3.9

Fitting ODE’s to Data

InSection3.8,itwasshownhowODEscanbesolvedusingMaxima or R.Oneof

the most important advantages of using R liesinthefactthatwithinR the solution

of ODEs can be coupled with R’s unsurpassed statistical capabilities. This is of

particular relevance when you are facing the important problem of ﬁtting an ODE to

experimental data. Consider, for example, Figure 3.13, which shows a comparison

of the alarm clock model (Equations 3.234–3.237) with the data in

room.csv.This

ﬁgure was obtained by hand tuning of the parameters in Equations 3.234–3.237,

that is, by hand tuning of r

, r

, T

,andT

(using the initial sensor temperature

room.csv for T

) until a good coincidence between the measurement data in

room.csv and T

as computed from Equations 3.234–3.237 was obtained. Besides

the fact that such a hand tuning of the parameters of on ODE can be tedious and

time consuming, it has the disadvantage that you do not get any information on the

statistical quality of the parameters estimated in this way. For example, it would be

nice to have conﬁdence intervals around the parameter estimates that would tell

us about the precision of these estimates, similar to the conﬁdence intervals in our

above discussion of nonlinear regression (Section 2.4).

3.9.1

Parameter Estimation in the Alarm Clock Model

Reconsidering the alarm clock model and data (Section 3.8.2.2), let us try to make

it better now. First of all, note that we have already discussed the parameter

estimation problem in Section 2.4. It was shown there how R’s

nls function can be

used to statistically estimate the parameters p = (p

, ..., p

)(s ∈ N) of a nonlinear

function y = f (x , p ) in a way that minimizes the distance between this function

and a given data set (x

, y

), ...,(x

, y

).Allwehavetodohereistoobservethatthe

discussion in Section 2.4 applies to our situation. In the alarm clock model, T

(t)is

the nonlinear function we want to ﬁt to the data in

room.csv using the parameter

vector p = (r

, r

, T

). Using the notation of Section 2.4, we can write T

(t)as

(t, p). The only difference between our situation and the discussion in Section