King M.R., Mody N.A. Numerical and Statistical Methods for Bioengineering: Applications in MATLAB

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

options = odeset(‘parameter1’, ‘value1’, ‘parameter2’, ‘value2’, ...)

For example, you can set the relative tolerance by doing the following. Create an

options structure,

options = odeset(‘RelTol’, 0.0001);

Now pass the newly created options structure to the ODE solver using the syntax

[t, y]= ode45(odefunc, [t0 tf], y0, options)

There are a variety of features of the ODE solvers that you can control using the

options structure. For more information, see MATLAB help.

Values of constants or parameters can be passed to the function that calculates the

right-hand side of Equation (7.26). These parameters should be listed after the

options argument. If you do not wish to supply an options structure to the ODE

solver, you should use [] in its place. When creating the deﬁning line of the function

odefunc, you must pass an additional argument called ﬂag, as shown below:

f = odefunc(t, y, ﬂag, p1, p2)

where p1 and p2 are additional parameters that are passed from the calling pro-

gram. See Box 7.3B for an example that uses this syntax .

Box 7.1C HIV–1 dynamics in the blood stream

We re-solve the coupled ODE problem using the MATLAB ode45 function:

¼ 0:05y

 0:5y

;

¼3:0y

;

¼ 30y

 3:0y

The initial concentrations of T*, V

, and V

are y

0ðÞ¼10; y

0ðÞ¼100; y

0ðÞ¼0.

The ode45 function can be called either in a script file or in the Command Window:

[t, y] = ode45(‘HIVODE’, 5.0, [10; 100; 0]);

The first five and last five time points for integration are

t(1:5)

t(end-4:end)

ans = ans =

1.0e-006 * 4.8440

0 4.8830

0.1675 4.9220

0.3349 4.9610

0.5024 5.0000

0.6698

The time step increases by five orders of magnitude from the beginning to the end of the interval. The

number of infected cells, infectious viral particles, and non-infectious viral particles at t = 5 days is

y(end, :)

ans –

0.9850 0.0000 11.8201

447

7.4 Adaptive step size methods

Box 7.2B Enzyme deactivation in industrial processes

Defining y

 S and y

 E

total

, we rewrite the first-order differential equations as follows:

¼

þ y

;

¼

1 þ y

;

where

¼ 21 M/s/e.u. (e.u.  enzyme units)

¼ 0:004 M

¼ 0:03/s.

The initial conditions are y

¼ 0:05 M and y

¼ 1  10

6

units.

enyzmedeact is an ODE function written to evaluate the right-hand side of the coupled ODE

equations. A script file is written to call ode45 and plot the result.

MATLAB program 7.10

% Use ode45 to solve the enzyme deactivation problem.

clear all

% Variables

y0(1) = 0.05; % M

y0(2) = 1e-6; % units

tf = 3600; % s

% Call ODE Solver

[t, y] = ode45(‘enzymedeact’, tf, y0);

% Plot Results

ﬁgure;

subplot(1,2,1)

plot(t,y(:,1),‘k-’,‘LineWidth’,2)

set(gca,‘FontSize’,16,‘LineWidth’,2)

xlabel(‘{\itt} in sec’)

ylabel(‘{\itS} in M’)

subplot(1,2,2)

plot(t,y(:,2),‘k-’,‘LineWidth’,2)

set(gca,‘FontSize’,16,‘LineWidth’,2)

xlabel(‘{\itt} in sec’)

ylabel(‘{\itE}_t_o_t_a_l in units’)

MATLAB program 7.11

function f = enzymedeact(t, y)

% Evaluate slope f(t,y) of coupled ODEs for enzyme deactivation

% Constants

k2 = 21; % M/s/e.u.

Km = 0.004; % M

kd = 0.03; % s-1

f = [-k2*y(1)*y(2)/(Km + y(1));-kd*y(2)/(1 + y(1)/Km)];

448

Numerical integration of ODEs

Box 7.3B Dynamics of an epidemic outbreak

Two sets of parameters are chosen to study the dynamics of the spread of disease in a population of

100 000.

Case (1) k =1× 10

−5

and γ = 4. This set of values represents a disease that spreads easily and

recovery is rapid, e.g. flu.

Case (2) k =5× 10

−8

and γ = 4000. This set of values represents a disease that spreads with

difficulty yet a person remains infectious for a longer period of time, e.g. diseases that are

spread by contact of bodily fluids.

The birth and death rates are β ¼ 0:0001; and μ ¼ 0:00001, respectively. Defining y

 S,

 D, and y

 I, we rewrite the first-order differential equations as

¼ky

þ β y

þ y

ðÞμy

;

¼ ky



;

 μy

The MATLAB ODE solver ode23 is used to solve this problem. For case 1, we observe from

Figure 7.11 that the disease becomes an epidemic and spreads through the population within a

month. As large sums of people become sick, recover, and then become immune, the disease

disappears because there is a lack of people left to spread the disease.

For the second case, because the recovery period is so long, we must modify the equations to

account for people who die while infected. From Figure 7.12, we see that, in 7 years, most of the

Figure 7.10 plots the drop in concentration with time. At t = 440 s, the active enzyme concentration

reduces to half its starting value. Note the quick deactivation of the enzyme in comparison to

consumption of the substrate. The enzymatic process needs to be modified to increase the life of the

enzyme.

Figure 7.10

Decrease in concentration of substrate and active enzyme with time.

0 2000 4000

0.038

0.04

0.042

0.044

0.046

0.048

0.05

0.052

t (s)

S (M)

0 2000 4000

0.2

0.4

0.6

0.8

× 10

−6

t (s)

total

(units)

449

7.4 Adaptive step size methods

Figure 7.11

Spread of disease among individuals in an isolated population for case 1.

0 50 100

× 10

t (days)

0 50 100

0.5

1.5

2.5

3.5

4.5

× 10

t (days) t (days)

0 50 100

Figure 7.12

Spread of disease among individuals in an isolated population for case 2.

0 2000 4000

× 10

t (days)

0 2000 4000

× 10

t (days) t (days)

0 2000 4000

0.5

1.5

2.5

3.5

4.5

× 10

450

Numerical integration of ODEs

MATLAB program 7.12

% Use ode23 to solve the epidemic outbreak problem.

clear all

% Variables

k = 1e-5; % infectiousness of disease, per person per contact per day

gamma = 4; % recovery period, days

% Initial conditions

y0(1) = 100000;

y0(2) = 100;

y0(3) = 0;

% Time interval

tf = 100; % days

% Call ODE Solver

[t, y] = ode23(‘epidemic’, tf, y0, [], k, gamma);

MATLAB program 7.13

function f = epidemic(t, y, ﬂag, k, gamma)

% Evaluate slope f(t,y) of coupled ODEs for dynamics of epidemic outbreak

% Constants

beta = 0.0001; % birth rate, per person per day

mu = 0.00001; % death rate, per person per day

f = [-k*y(1)*y(2) + beta*(y(1)+y(3)) - mu*y(1); ...

k*y(1)*y(2) - 1/gamma*y(2); ...

1/gamma*y(2) - mu*y(3)];

7.5 Multistep ODE solvers

The Euler and the RK methods of numerical ODE integ ration are called one-step

methods because they use the solution y

at the current time step t

to estimate the

solution at the next time step t

kþ1

. In RK methods, the formula to obtain the

numerical approximation to yt

kþ1

ðÞrequires evaluation of the slope at one or

more time points within the subinterval t

; t

kþ1

½over which integration is currently

being performed. Information on the solution prior to t

is not used to determine

kþ1

in these methods. Thus, in these methods the integration moves forward by

stepping from the current step to the next step, i.e. by a single step.

population succumbs to the disease. The disease however does not disappear, but remains embedded

in the population, even after 10 years.

Note that these models are highly simplified. The equations assume that each person has an equal

probability of coming into contact with another person. It has been found that some epidemic models

are better explained using complex network theory that uses a power law distribution to model the

number of people that an individual of a population comes into contact with. Complex network theory

assigns some people the status of “hubs.” Disease spreads most easily through “hubs” since they are

the most well-connected (Barabasi, 2003).

451

7.5 Multistep ODE solvers

Much useful information is contained in the solution at time points earlier than t

The change in the values of y

; y

k1

; y

k2

; ... at previous time points provides a

measure of the rate of change in the slope f

ðt; yÞ and higher derivatives of y. The

knowledge of prior solutions over a larger time interval consisting of severa l time

steps can be used to obtain more accurate app roximations to the solution at the next

time point. ODE integration methods that use not only the so lution at the current

step t

to approximate the solution at the next step, but also previous solution values

at t

k1

; t

k2

; ...are called multistep methods.

A general-purpose numerical integration formula for a three-step method is

kþ1

¼p

þ p

k1

þ p

k2

þ hq

kþ1

; y

kþ1

ðÞþq

; y

ðÞþq

k1

; y

k1

ðÞþq

k2

; y

k2

ðÞ½;

where p

; p

; q

, and q

are constants that are derived using the Taylor

series.

One adv antage of multistep methods for ODE integration is that only one

evaluation of fðt; yÞ is required for each step regardless of the order of the method.

Function values at previous time points are already available from previous inte-

gration steps in the multistep method. However, multistep methods are not self-

starting, because several previous solution points are required by the formula to

calculate y

kþ1

. An initial condition, in most cases, consists of the value of y at one

initial time point only. To begin the integrati on, a one-step method of order equal to

that of the multistep method is used to generate the solution at the ﬁrst few time

points, and after that the multistep method can take over.

There are two kinds of multistep methods: explicit methods and implicit methods.

Explicit multistep methods do not require an evaluation of the slope at t

kþ1

to ﬁnd

kþ1

. Implicit multistep methods include a term for the slope at t

kþ1

on the right-

hand side of the integration formula. The implicit methods have better stability

properties, but are harder to execute than explicit methods. Unless the slope function

is linear in y, a nonlinear root-ﬁnding algori thm must be used to search for y

kþ1

variety of multistep algorithms have been developed. We discuss the Adams methods

of explicit and implicit type for multistep ODE integration. MATLAB offers an

ODE solver based on an Adams routine.

7.5.1 Adams methods

The Adams–Bashforth methods are explicit multistep ODE integration schemes in

which the number of steps taken during the integration at each time step is equal to

the order of the method. Keep in mind that the order of a method is equal to the

order of the global truncation error produced by that method. The integration

formula for an m-step Adams–Bashforth method is given by

kþ1

¼ y

þ h

j¼1

kjþ1

; y

kjþ1



: (7:44)

The formula for an m-step Adams–Bashforth method can be derived using the

Taylor seri es expansion about yt

ðÞwith step size h, as follows:

kþ1

ðÞ¼yt

ðÞþft

; yt

ðÞðÞh þ

; yt

ðÞðÞh

; yt

ðÞðÞh

þ Oh



452

Numerical integration of ODEs

Substituting numerical differentiation formulas for derivatives of f generates a

multistep difference formula. This is demonstrated for a two-step Adams–

Bashforth method.

We can approximate the ﬁrst derivative of f using the ﬁrst-order backward differ-

ence formula (Equation (1.22)),

; yt

ðÞðÞ

; yt

ðÞðÞft

k1

; yt

k1

ðÞðÞ

; yt

ðÞðÞþOh



Substituting the backward difference formula above into the Taylor expansion,

we get

kþ1

ðÞ¼yt

ðÞþft

; yt

ðÞðÞh þ

; yt

ðÞðÞft

k1

; yt

k1

ðÞðÞðÞ

; yt

ðÞðÞþOh



This simpliﬁes to

kþ1

ðÞ¼yt

ðÞþ

3ft

; yt

ðÞðÞft

k1

; yt

k1

ðÞðÞðÞþ

ξ; y ξðÞðÞ;

where ξ 2 t

k1

; t

kþ1

½. Writing ft

; yt

ðÞðÞas f

, the second-order Adams–Bashforth

formula is given by

kþ1

¼ y

 f

k1

ðÞ: (7:45)

Upon substituting a second-order backward difference approximation for the ﬁrst

derivative of f and a ﬁrst-order backward difference formula for the second deriv-

ative of f, we obtain the third-order Adams–Bashforth formula:

kþ1

¼ y

23f

 16f

k1

þ 5f

k2

ðÞ: (7:46)

In this way, higher-order formulas for the Adams–Bashforth ODE solver can be

derived.

The Adams–Moulton methods are implicit multistep ODE integration schemes.

For these methods, the number of steps m taken to advance to the next time point is

one less than the order of the method. The integration formula for an m-step

Adams–Moulton method is given by

kþ1

¼ y

þ h

j¼0

kþ1j

; y

kþ1j



: (7 :47)

The form ula for an m-step Adams–Moulton method can be derived using a Taylor

series expansion about yt

kþ1

ðÞwith step size −h, as follows:

ðÞ¼yt

kþ1

ðÞft

kþ1

; yt

kþ1

ðÞðÞh þ

kþ1

; yt

kþ1

ðÞðÞh



kþ1

; yt

kþ1

ðÞðÞh

þ Oh



Substituting numerical differentiation formulas for derivatives of f generates a

multistep difference formula. Retaining only the ﬁrst two terms on the right-hand

side, we recover Euler’s implicit method, a ﬁrst-order method. Approximating the

ﬁrst derivative of f in the Taylor series expansion with a ﬁrst-order backward

453

7.5 Multistep ODE solvers

difference formula (Equation (1.22)), and dropping the Oh



terms, we obtain the

numerical integration formula for the modiﬁed Euler method, a secon d-order

method (see Section 7.2.3):

kþ1

¼ y

þ f

kþ1

ðÞ: (7:48)

The associated Oh



error term is ðh

=12Þf

ξ; y ξðÞðÞ. Try to derive the error term

yourself.

Upon substituting a second-order backward difference approximation for the

ﬁrst derivative of f and a ﬁrst-order backward difference formula for the second

derivative of f, we obtain the third-order Adams–Moulton formula :

kþ1

¼ y

kþ1

þ 8f

 f

k1

ðÞ: (7:49)

In this manner, higher-order formulas for the Adams–Moulton ODE solver can be

derived.

Upon comparing the magnitude of the error term of the second -order Adams–

Bashforth method and the second-or der Adams–Moulton method, we notice that

the local truncation error is smaller for the implicit Adams method. This is also the

case for higher-order explicit and implicit Adams methods. The smaller error lends

itself to smaller global truncation error and impro ved numerical stability.

Remember that instability of a numerical solution occurs when the error generated

at each step is magniﬁed over successive time steps by the integration formula. The

second-order Adams–Moulton method is stable for all step sizes when the system of

ODEs is well-posed.

However, implicit multistep methods of order greater than 2

are not stable for all step sizes, and can produce truncation errors that grow without

bound when large step sizes are used (Grasselli and Pelinovsky, 2008).

7.5.2 Predictor–corrector methods

Implicit numerical ODE integration methods have superior stability characteristics

and smaller local truncation errors compared to explicit methods. However, execu-

tion of implicit methods can be cumbersome since an iterative root-ﬁnding algo-

rithm must be employed in general to ﬁnd the solution at each time step. Also, to

begin the root- ﬁnding process, one must provide one or two guesses that lie close

to the true solution of the implicit ODE integration formula. Rather than having to

resort to iterative numerical methods for solving implicit multistep formulas, one

can obtain a prediction of the solution using the explicit multistep formula. The

explicit multistep ODE solver is called the predictor. The prediction y

ð0Þ

kþ1

is plugged

into the right-hand side of the implicit multistep ODE solver to calculate the slope at

ðt

kþ1

; y

ð0Þ

kþ1

Þ. The implicit formula imparts a correction on the prediction of the

solution made by the predictor. Therefore, the implicit ODE integration formula

is called the corrector.

Predictor and corrector formulas of the same order are paired with each other. An

example of a predictor–corrector pair is the second-order Adams–Bashforth method

paired with the second-order Adams–Moulton method shown below:

Well-posed ODEs have a unique solution that changes only slightly when the initial condition is

perturbed. An ODE whose solution ﬂuctuates widely with small changes in the initial condition

constitutes an ill-posed problem.

454

Numerical integration of ODEs

predictor:

ð0Þ

kþ1

¼ y

 f

k1

ðÞ;

corrector:

ð1Þ

kþ1

¼ y

ðft

; y

ðÞþfðt

kþ1

; y

ð0Þ

kþ1

ÞÞ.

A numerical integration method that uses a predictor equation with error Oh

ðÞto

predict the solution and a corrector equation with error Oh

ðÞto improve upon the

solution is called a predictor–corrector method. In some predictor–corrector algo-

rithms, the corrector formula is applied once per time step and y

ð1Þ

kþ1

is used as the

ﬁnal approximation of the solution at t

kþ1

. Alternatively, the approximation y

ð1Þ

kþ1

can be supplied back to the right-hand side of the corrector equation to obtain an

improved corrector approximation of the solution y

ð2Þ

kþ1

. This is analogous to per-

forming ﬁxed-point iteration (see Chapter 5). Convergence of the corrector equation

is obtained when

ðiÞ

kþ1

 y

ði1Þ

kþ1



ðiÞ

kþ1

 E:

An advantage of the predictor–corrector algorithm is that it also provides an

estimate of the local truncation error at each time step. The exact solution at t

kþ1

can be expressed in terms of the predictor solut ion,

kþ1

ðÞ¼y

ð0Þ

kþ1

þ c

nþ1

ðÞ;

and the corrector solution as

kþ1

ðÞ¼y

ðiÞ

kþ1

þ c

nþ1

ðÞ:

By assuming that ξ

 ξ

¼ ξ, we can subtract the two expressions above to get

ðiÞ

kþ1

 y

ð0Þ

kþ1

¼ c

 c

ðÞh

nþ1

ξðÞ:

Simple algebraic manipulations yield

 c

ðiÞ

kþ1

 y

ð0Þ

kþ1



¼ c

nþ1

ξðÞ: (7:50)

Thus, the approximations yielded by the predictor and corrector formulas allow us

to estimate the local truncation error at each time step. Speciﬁcally for the second-

order Adams–Bashforth–Moulton method:

; c

¼

Thus, for the second-order predictor–corrector method based on Adams formulas,

the local truncation error at t

kþ1

can be estimated with the formula



ðiÞ

kþ1

 y

ð0Þ

kþ1



As discussed in Section 7.4, local error estimates come in handy for optimizing the

time step such that the local error is always well within the tolerance limit. The local

truncation error given by Equation (7.50) can be compared with the tolerance limit

to determine if a change in step size is warranted. The condition imposed by the

tolerance limit E is given by c

nþ1

opt

nþ1

ξðÞE. Substituting y

nþ1

ξðÞE=c

nþ1

opt

into

Equation (7.50), we obtain, after some rearrangement,

455

7.5 Multistep ODE solvers

opt



nþ1



E c

 c

ðÞ

ðiÞ

kþ1

 y

ð0Þ

kþ1



: (7:51)

An adjustment to the step size requires recalculation of the solution at new, equally

spaced time points. For example, if the solutions at t

; t

k1

; and t

k2

are used to

calculate y

kþ1

, and if the time step is halved during the interval t

; t

kþ1

½, then, to

generate the approximation for y

kþ1=2

, the solutions at t

; t

k1=2

; and t

k1

are

required. Also, y

k1=2

must be determined using interpolation formulas of the

same order as the predictor–corrector method. Step size optimization for multistep

methods entails more computational work than for one-step methods.

Using MATLAB

The MATLAB ODE solver ode113 is based on the Adams–Bashforth–Moulton

method. This solver is preferred for problems with strict tolerance limits since the

multistep method on which the solver is based produces smaller local truncation

errors that the one-step RK method of corresponding order.

7.6 Stability and stiff equations

In earlier sections of this chapter, we discussed the stability properties of one-step

methods, speciﬁcally the Euler and explicit RK methods (implicit RK methods have

been devised but are not discussed in this book.) When solving a well-posed system

of ODEs, the ODE solver that yields a solution whose relative error is kept in control

even after a large number of steps is said to be numerically stable for that particular

problem. If the relative error of the solution at each time step is magniﬁed at the next

time step, the numerical solution will diverge from the exact solution. The solution

Box 7.4B Microbial population dynamics

The coupled ODEs that describe predator–prey population dynamics in a mixed culture microbial

system are presented below after making the following substitutions: y

¼ S; y

¼ N

; y

¼ N

1;0

 y





;max

þ y

;

¼

;max

þ y



;max

þ y

;

¼

;max

þ y

Using the kinetic parameters that define the predator–prey system (Tsuchiya et al., 1972), we attempt to

solve this system of ODEs using ode113. However, the numerical algorithm is unstable and the

solution produced by the ODE solver diverges with time to infinity. For the given set of parameters, this

system of equations is difficult to integrate because of sudden changes in the substrate concentration

and microbial density with time. Such ODEs are called stiff differential equations, and are the topic of

Section 7.6.

456

Numerical integration of ODEs