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

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

produced by the numerical method will grow increasingly inaccurate with time.

The

ODE solver is said to be numerically unstable for the problem at hand. An ODE

solver that successfully solves one ODE problem may not necessarily produce a

convergent solution for another ODE problem, when using the same step size. The

numerical stability of an ODE solver is deﬁned by a range of step sizes within which a

non-divergent solution (solution does not blow up at long times) is guaranteed. This

step size range is not only ODE solver speciﬁc, but also problem speciﬁc, and thus

varies from problem to problem for the same solver. Some numerical integration

schemes have be tter stability properties than others (i.e. wider step size range for

producing a stable solution), and selection of the appropriate ODE solver will

depend on the nature of the ODE problem.

Some ODE problems are more “difﬁcult” to solve. Solutions that have rapidly

changing slopes can cause difﬁculties during numerical integration. A solution that

has a fast decaying component, i.e. an exponential term, e

αt

, where α  0, will

change rapidly over a short period. The solution is said to have a large transient.

ODEs that have large transients are called stiff equations. In certain processes, often

multiple phenomena with different time scales will contribute to the behavior of a

system. Some of the contributing factors will rapidly decay to zero, while others will

evolve slowly with time. The transient parts of the solution eventually disappear and

give way to a steady state solution, which is unchanging with time. The single ODE

initial-value problem given by

¼ ft; yðÞ; y 0ðÞ¼y

;

can be linearized to the form dy=dt ¼ αy, such that α ¼ df=dy. Note that αðtÞ varies

in magnitude as a function of time. If at any point α becomes large and negative, the

ODE is considered stiff.

Unless the time step is chosen to be vanishingly small throughout the entire

integration interval, the decay of a transient part of the solution will not occur

when using certain ODE solvers. Large errors can accumulate due to the magniﬁ-

cation of errors at each time step when, at any time during the integration, the step

size jumps outside the narrow range that guarantees a stable, convergent solution.

Solvers that vary the step size solely based on local truncation error are particularly

prone to this problem when solving stiff equations.

Ordinary ODE solvers must use prohibitively small step sizes to solve stiff prob-

lems. Therefore, we turn to special solution techniques that cater to this category of

differential equations. Stiff ODE solvers have better stability properties and can

handle stiff problems using a step size that is practical. Very small step sizes increase

the solution time considerably, and may reach the limit of round-off error. In the

latter case, an accurate solution cannot be found. Stiffness is mo re commonly

observed in systems of ODEs rather than in single ODE problems.

For a system of linearized, ﬁrst-order ordinary differential equations, the set of

equations can be compactly represented as

¼ Jy; y 0ðÞ¼y

;

Although we have referred to time as the independent variable in an initial-value problem, the methods

developed for solving initial-value problems apply equally well to any continuous independent variable,

such as distance x, speed s, or pressure P.

457

7.6 Stability and stiff equations

where

J ¼

∂f

∂y

∂f

∂y

...

∂f

∂y

∂f

∂y

∂f

∂y

...

∂f

∂y

∂f

∂y

∂f

∂y

...

∂f

∂y

J is called the Jacobian, and was introduced in Chapter 5. The exact solution to this

linear problem is

y ¼ e

: (7:52)

Before we can simplify Equat ion (7.52), we need to introduce some important linear

algebra concepts. The equation Ax ¼ λx comprises an eigenvalue problem. This

linear matrix problem represents a linear transformation of x that simply changes

the size of x (i.e. elongation or shrinkage) but not its direction. To obtain non-zero

solutions of x, A  λIjjis set as equal to zero, and the values of λ for which the

determinant A  λI

equals zero are called the eigenvalues of A. The solution x for

each of the n eigenvalues are called the eigenvectors of A.

If J (a square matrix of dimensions n × n)hasn distinct eigenvalues λ

and n distinct

eigenvectors x

, then it is called a perfect matrix, and the above exponential equation

can be simpliﬁed by expressing J in terms of its eigenvectors and eigenvalues. How this

can be done is not shown here since it requires the development of linear algebra

concepts that are beyond the scope of this book. (The proof is explained in Chapter 5

of Davis and Thomson (2000).) The above equation simpliﬁes to

y ¼

i¼1

: (7:53)

If an eigenva lue is complex, then the real part needs to be negative for stability. If the

real part of any of the n eigenvalues λ

is positive, the exact solution will blow up. If

all λ

are negative, the exact solution is bounded for all time t. Thus, the sign of the

eigenvalues of the Jacobian reveal the stabi lity of the true solution, i.e. whether the

solution remains stable for all times, or becomes singular at some point in time. If the

eigenvalues of the Jacobian are all negative and differ by several orders of magni-

tude, then the ODE problem is termed stiff. The smallest step size for integration is

set by the eigenvalue of largest magnitude.

Using MATLAB

MATLAB offers four ODE solvers for stiff equations, ode15s, ode23s, ode23t,

and ode23tb .

ode15s is a mu ltistep solver. It is suggested that ode15s be tried when ode45 or

ode113 fails and the problem appears to be sti ff.

ode23s is a one-step solver and can solve some stiff problems for which ode15s fails.

ode23t and ode23tb both use the modiﬁed Euler method, also call ed the trape-

zoidal rule or the Adams–Moulton second-order method; ode23tb also imple-

ments an implicit RK formula.

458

Numerical integration of ODEs

For stiff ODE solvers, you can supply a function that evaluates the Jacobian matrix

to the solver (using options) to speed up the calculations. If a Jacobian matr ix is

not provided for the function, the components of the Jacobian are estimated by the

stiff ODE solver using ﬁnite difference formulas. See MATLAB help for further

details.

Box 7.4C Microbial population dynamics

We re-solve the coupled ODEs presented below that describe predator–prey population dynamics in a

mixed culture microbial system:

1;0

 y





;max

þ y

;

¼

;max

þ y



;max

þ y

;

¼

;max

þ y

;

using the stiff ODE solver ode15s. For this parameter set, we get an oscillatory response of the

system, shown in Figure 7.13. Changes in the initial values of the variables N

and N

do not affect the

long-term oscillatory properties of the system. On the other hand, the nature of the oscillations that are

characteristic of the population dynamics predicted by the Lotka–Volterra predator–prey model are

affected by changes in initial predator or prey numbers.

There are a number of operating variables that can be changed, such as the feed-to-volume ratio and

the substrate concentration feed. At time t = 50 days, we introduce a step change in the feed glucose

Figure 7.13

Evolution of number densities of E.coli and Dictyostelium discoideum, and substrate concentration with time in a

chemostat for S

= 0.5 mg/ml.

0 5 10 15 20 25 30 35 40 45 50

0.5

(mg/ml)

0 5 10 15 20 25 30 35 40 45 50

× 10

0 5 10 15 20 25 30 35 40 45 50

× 10

t (days)

(cells/ml)

459

7.6 Stability and stiff equations

MATLAB program 7.14

% Use ode15s to solve the microbial population dynamics problem.

clear all

% Variables

S0 = 0.5; % glucose concentration in feed mg/ml

% Initial conditions

y0(1) = 0.5; % mg glucose/ml

y0(2) = 13e8; % bacteria/ml

y0(3) = 4e5; % amoeba/ml

% Time interval

tf = 1200; % hr

% Call ODE Solver

[t, y] = ode15s(‘popdynamics’, tf, y0, [], S0);

concentration from 0.5 mg/ml to 0.1 mg/ml. The behavior of the system for the next 50 days is shown in

Figure 7.14.

The system quickly stabilizes to steady state values. The long-term system behavior is a function of

the F/V ratio and the S

concentration. Operating at F=V ¼ð1=5Þ/s and S

¼ 0:5 mg/ml leads to

washout of amoebae (predator).

In the program listing, the plotting commands have been removed for brevity.

Figure 7.14

Evolution of number densities of E.coli and Dictyostelium discoideum amoebae, and substrate concentration with

time, in a chemostat for S

= 0.1 mg/ml.

0 5 10 15 20 25 30 35 40 45 50

0.05

0.1

(mg/ml)

0 5 10 15 20 25 30 35 40 45 50

× 10

(cells/ml)

0 5 10 15 20 25 30 35 40 45 50

× 10

t (days)

(cells/ml)

460

Numerical integration of ODEs

% Step change in substrate concentration in feed

S0 = 0.1; %mg glucose/ml

y0 = y(end, :); % Values of variables at last time point

% Call ODE Solver

[t, y] = ode15s(‘popdynamics’, tf, y0, [], S0);

MATLAB program 7.15

function f = popdynamics(t, y, ﬂag, S0)

% Evaluate slope f(t,y) of coupled ODEs for microbial population dynamics

% Constants

FbyV = 1/16; % feed to chemostat volume ratio (/hr)

muN1max = 0.25; % max speciﬁc growth rate for E coli (/hr)

muN2max = 0.24; % max speciﬁc growth rate for amoeba (/hr)

KN1 = 5e-4; % saturation constant (mg glucose/ml)

KN2 = 4e8; % saturation constant (bacteria/ml)

invYN1S = 3.3e-10; % reciprocal yield factor (mg glucose/bacteria)

invYN2N1 = 1.4e3; % reciprocal yield factor (bacteria/amoeba)

f = [FbyV*(S0 - y(1)) - invYN1S*muN1max*y(1)*y(2)/(KN1 + y(1)); ...

-FbyV*y(2) + muN1max*y(1)*y(2)/(KN1 + y(1)) - ...

invYN2N1*muN2max*y(2)*y(3)/(KN2 + y(2)); ...

-FbyV*y(3) + muN2max*y(2)*y(3)/(KN2 + y(2))];

7.7 Shooting method for boundary-value problems

A unique solution for an ordinary differential equation of order n is pos sible only

when n constraints are speciﬁed along with the equation. This is because, if a solution

exists for the ODE, the constraints are needed to assign values to the n integrati on

constants of the general solution, thus making the solution unique. For an initial-

value problem, the n constraints are located at the initial point of integration. Such

problems are often encountered when integration is performed over the time

domain. When the interval of integration lies in the physical space domain, it is

likely that some of the n constraints are avail able at the starting point of integration,

i.e. at one end of the boundary, while the remaining constraints are speciﬁed at the

endpoint of integration at the other end of the boundary. Constraints that are

speciﬁed at two or more boundary points of the interval of integration are called

boundary conditions. ODEs of second order and higher that are supplied with

boundary conditions are called boundary-val ue problems (BVPs).

Consider the second-order ODE



þ y ¼ ρ:

To solve this non-homogeneous ODE problem, y is integrated with respect to x. The

interval of integration is a  x  b, and ρ is some constant.

There are several types of boundary conditions that can be speciﬁed for a second-

order ODE. At any boundary endpoint, three types of boundary conditions are

possible.

461

7.7 Shooting method for boundary-value problems

(1) y ¼ α A boundary condition that speciﬁes the value of the dependent variable at an

endpoint is called a Dirichlet boundary condition.

(2) y

¼ α A boundary condition that speciﬁes the value of the derivative of the

dependent variable at an endpoint is called a Neum ann boundary condition.

(3) c

þ c

y ¼ c

A boundary condition that linearly combines the dependent vari-

able and its derivative at an endpoint is called a mixed boundary condition.

For any second-order boundary-value ODE problem, one boundary condition is

speciﬁed at x ¼ a and another condition at x ¼ b. The boundary-value problem

 y

þ y ¼ ρ; yaðÞ¼α; ybðÞ¼β

is said to be have two Dirichlet boundary conditions.

Several numerical methods have been designed speciﬁcally to solve boundary-

value problems. In one popular technique, called the ﬁnite difference method, the

integration interval is divided into many equally spaced discrete points called

nodes. The derivatives of y in the equation are substituted by ﬁnite difference

approximations. This creates one algebraic equation for each node in the inter val.

If the ODE system is linear, the numerically equivalent system of algebraic equa-

tions is linear, and the methods discussed in Chapter 2 for solving linear equations

are used to solve the problem. If the ODE system is nonlinear in y, the resulting

algebraic equations are also nonlinear, and iterative methods such as Newton’s

method must be used.

Another method for solving boundary-value problems is the shooting method.

This method usually provides better accuracy than ﬁnite difference methods.

However, this method can only be applied to solve ODEs and not partial differential

equations (PDEs). On the other hand, ﬁnite difference methods are easier to use for

the higher-order ODEs and are commonly used to solve PDEs. In the shooting

method, the nth-order ODE is ﬁrst converted to a system of ﬁrst-order ODEs. Any

unspeciﬁed initial conditions for the set of ﬁrst-order equations are guessed, and the

ODE system is treated as an initial-value pro blem. The ODE integration schemes

described earlier in this chapter are used to solve the problem for the set of known

and assumed initial conditions. The ﬁnal endpoint values of the variables obtained

from the solution are compared with the actual boundary conditions supplied. The

magnitude of deviation of the numerical solution at x ¼ b from the desired solution

is used to alter the assumed initial conditions, and the system of ODEs is solved

again. Thus, we “shoot” from x ¼ a (by assuming values for unspeciﬁed boundary

conditions at the ﬁrst endpoint) for a desired solution at x ¼ b, our ﬁxed target.

Every time we shoot from the starting point, we reﬁne our approximation of the

unknown boundary conditions at x ¼ a by taking into account the difference

between the numerically obtained and the known boundary conditions at x ¼ b.

The shooting method can be combined with a method of nonlinear root ﬁnding (see

Chapter 5) for improving the future iterations of the initial boundary condition.

Boundary-value problems for ODEs can be classiﬁed as linear or nonlinear, based

on the dependence of the differential equations on y and its derivatives. Solution

techniques for linear boundary-value problems are different from the methods used

for nonlinear problems. Linear ODEs are solved in a straightforward manner with-

out iteration, while nonlinear BVPs require an iterative sequence to converge upon a

solution within the tolerance provided. In this section, we discuss shooting methods

for second-order ODE boundary-value problems. Higher-order boundary-value

problems can be solved by extending the methods discussed here.

462

Numerical integration of ODEs

7.7.1 Linear ODEs

A linear second-order boundary value problem is of the form

¼ pxðÞy

þ qxðÞy þ rxðÞ; x 2 a; b½; (7:54a)

where the independent variable is x. It is assumed that a unique solution exists for the

boundary conditions supplied. Suppos e the ODE is subjected to the following

boundary conditions:

yaðÞ¼α; ybðÞ¼β: (7:54b)

A linear problem can be reformulated as the sum of two simpler linear problems,

each of which can be determined separately. The superposition of the solutions of

the two simpler linear problems yields the solution to the original problem . This is

the strategy sought for solving a linear boundary-value problem. We express

Equation (7.54) as the sum of two linear initial-value ODE problems. The ﬁrst linear

problem has the same ODE as given in Equation (7.54a):

¼ pxðÞu

þ qxðÞu þ rxðÞ; uaðÞ¼α; u

aðÞ¼0: (7:55)

One of the two initial conditions is taken from the actual boundary conditions

speciﬁed in Equation (7.54b). The second boundary condition assumes that the ﬁrst

derivative is zero at the left boundary of the interval. The second linear ODE is the

homogeneous form of Equat ion (7.54a):

¼ pxðÞv

þ qxðÞv; vaðÞ¼0; u

aðÞ¼1: (7:56)

The initial conditions for the second linear problem are chosen so that the solution of

the origina l ODE (Equation (7.54)) is

yxðÞ¼uxðÞþcv xðÞ; (7:57)

where c is a constant that needs to be determined. Before we proceed further, let us

conﬁrm that Equation (7.57) is indeed the solution to Equation (7.54a).

Differentiating Equation (7.57) twice, we obtain

¼ u

þ cv

¼ pxðÞu

þ qxðÞu þ rxðÞþcpxðÞv

þ qxðÞvðÞ:

Rearranging,

¼ pu

þ cv

ðÞ

þ quþ cv

ðÞ

þ r:

We recover Equation (7.54a),

¼ pxðÞy

þ qxðÞy þ rxðÞ:

Thus, the solution given by Equation (7.57) satisﬁes Equation (7.54a).

Equation (7.57) must also satisfy the boundary conditions given by Equation

(7.54b) at x ¼ a:

uaðÞþcv aðÞ¼α þ c  0 ¼ α:

We can use the boundary condition at x ¼ b to determine c. Since

ðÞ

þ cv b

ðÞ

¼ β;

c ¼

β  u ðbÞ

vðbÞ

463

7.7 Shooting method for boundary-value problems

If vbðÞ¼0, the solution of the homogeneous equation can be either v ¼ 0 (in which

case y ¼ u is the solution) or a non-trivial solut ion that satisﬁes homogeneous

boundary conditions. When vbðÞ¼0, the solution to Equation (7.54) may not be

unique.

For vbðÞ6¼ 0, the solution to the second-order linear ODE with boundary con-

ditions given by Equation (7.54b) is given by

yxðÞ¼uxðÞþ

β  uðbÞ

vðbÞ

vxðÞ: (7:58)

To obtain a solution for a linear second-order boundary-value problem subject

to two Dirichlet boundary conditions, we solve the initial-value problems

given by Equations (7.55) and ( 7.56), and then combine the solutions using

Equation (7.58).

If the boundary condition at x ¼ b is of Neumann type ( y

bðÞ¼β) or of mixed

type, and the boundary condition at x ¼ a is of Dirichlet type, then the two initial-

value problems (Equations (7.55) and (7.56)) that must be solved to obtain u and v

remain the same. The method to calculate c is, for a Neumann boundary condition at

x ¼ b , given by

bðÞ¼u

bðÞþcv

bðÞ¼β

c ¼

β  u

bðÞ

The initial conditions of the two initial-value ODEs in Equations (7.55) and (7.56)

will depend on the type of boundary condition speciﬁed at x ¼ a for the original

boundary-value problem.

7.7.2 Non-linear ODEs

Our goal is to ﬁnd a solution to the nonlinear second-order boundary-value problem

¼ fx; y; y

ðÞ; a  x  b; yaðÞ¼α; ybðÞ¼β: (7:59)

Because the ODE has a nonlinear dependence on y and/or y

, we cannot simplify the

problem by converting Equation (7.59) into a linear combination of simpler ODEs,

i.e. two initial-value problems. Equation (7.59) must be solved iteratively by starting

with a guessed value s for the unspeciﬁed initial condition, y

aðÞ, and improving our

estimate of s until the error in our numerical estimat e of yðb Þ is reduced to below the

tolerance limit E. The deviation of yðbÞ from its actual value β is a nonlinear function

of the guessed slope s at x ¼ a. We wish to minimize the error, ybðÞβ, which can be

written as a function of s:

ybðÞβ ¼ zsðÞ: (7:60)

We seek a value of s such that zsðÞ¼0. To solve Equation (7.60), we use a nonlinear

root-ﬁnding algorithm. We begin by solving the initial-value problem,

¼ fx; y; y

ðÞ; a  x  b; yaðÞ¼α; y

aðÞ¼s

;

464

Numerical integration of ODEs

using the methods discussed in Sections 7.2– 7.6, to obtain a solution y

ðxÞ.If

¼ y

bðÞβ  E, the solution y

ðxÞ is accepted.

Otherwise, the initial-value problem is solved for another guess value, s

, for the

slope at x ¼ a,

¼ fx; y; y

ðÞ; a  x  b; yaðÞ¼α; y

aðÞ¼s

;

to yield another solution y

ðxÞ corresponding to the guessed value y

aðÞ¼s

.If

¼ y

bðÞβ  E, then the solution y

ðxÞ is accepted. If neither guess for the

slope s

or s

produce a solution that meets the boundary condition criterion at

x ¼ b , we must generate an impr oved guess s

. The secant method, which approx-

imates the function z within an interval Δs as a straight line, can be used to

produce a next guess for y

aðÞ. If the two points on the functi on curve, zsðÞ, that

deﬁne the interval, Δs, are s

; z

ðÞand s

; z

ðÞ, then the slope of the straight line

joining these two points is

m ¼

bðÞz

bðÞ

 s

A third point lies on the straight line with slope m such that it crosses the x-axis

at the point s

; 0ðÞ, where s

is our approximation for the actual slope that

satisﬁes zsðÞ¼0. The straight-line equation passing through these three points is

given by

0  z

bðÞ¼ms

 s

ðÞ

¼ s



 s

ðÞz

ðbÞ

bðÞz

bðÞ

We can generalize the formula above as

iþ1

¼ s



 s

i1

ðÞz

ðbÞ

bðÞz

i1

bðÞ

; (7:61)

which is solved iteratively until zs

iþ1

ðÞE. A nonlinear equation usually admits

several solutions. Therefore, the choice of the starting guess value is critical for

ﬁnding the correct solution.

Newton’s method can also be used to iteratively solve for s. This root-ﬁnding

algorithm requires an analytical expression for the derivative of yðbÞ with respect

to s. Let’s write the function y at x ¼ b as y ðb; sÞ, since it is a function of the initial

slope s. Since the analytical form of zðsÞ is unknown, we cannot directly take the

derivative of yðb; sÞ to obtain dy=dsðb; sÞ. Instead, we create another second-order

initial-value problem whose solution gives us ∂y=∂s ðx; sÞ. Construction of this

companion second-order ODE requires obtaining the analytical form of the deriv-

atives of fx; y; y

ðÞwith respect to y and y

. While Newton’s formula results in

faster convergence, i.e. fewer iterations to obtain the desired solution, two second-

order initial-value problems must be solved simultaneously for the slope s

obtain an improved guess value s

iþ1

. An explanation of the derivation and appli-

cation of the shooting method using Newton’s formula can be found in Burden

and Faires (2005).

465

7.7 Shooting method for boundary-value problems

Box 7.5 Controlled-release drug delivery using biodegradable polymeric microspheres

Conventional dosage methods, such as oral delivery and injection, are not ideally suited for sustained

drug delivery to diseased tissues or organs of the body. Ingesting a tablet or injecting a bolus of drug

dissolved in aqueous solution causes the drug to be rapidly released into the body. Drug concentrations

in blood or tissue can rise to nearly toxic levels followed by a drop to ineffective levels. The duration over

which an optimum drug concentration range is maintained in the body, i.e. the period of deriving

maximum therapeutic benefit, may be too short. To maintain an effective drug concentration in the

blood, high drug dosages, as well as frequent administration, become necessary.

When a drug is orally administered, in order to be available to different tissues, the drug must first

enter the bloodstream after absorption in the gastrointestinal tract. The rate and extent of absorption may

vary greatly, depending on the physical and chemical form of the drug, presence or absence of food,

posture, pH of gastrointestinal fluids, duration of time spent in the esophagus and stomach, and drug

interactions. Uncontrolled rapid release of drug can cause local gastrointestinal toxicity, while slow or

incomplete absorption may prevent realization of therapeutic benefit.

In recent years, more sophisticated drugs in the form of protein-based and DNA-based com-

pounds have been introduced. However, oral delivery of proteins to the systemic circulation is

particularly challenging. For the drug to be of any use, the proteins must be able to pass through the

gastrointestinal tract without being enzymatically degraded. Many of these drugs have a small

therapeutic concentration range, with the toxic concentration range close to the therapeutic range.

Research interest is focused on the development of controlled-release drug-delivery systems that can

maintain the therapeutic efficacy of such drugs. Controlling the precise level of drug in the body

reduces side-effects, lowers dosage requirements and frequency, and enables a predictable and

extended duration of action. Current controlled-release drug-delivery technologies include trans-

dermal patches, implants, microencapsulation, and inhaled and injectable sustained-release peptide/

protein drugs.

An example of a controlled-release drug-delivery system is the polymer microsphere. Injection of

particulate suspensions of drug-loaded biodegradable polymeric spheres is a convenient method to

deliver hormonal proteins or peptides in a controlled manner. The drug dissolves into the

surrounding medium at a pre-determined rate governed by the diffusion of drug out of the polymer

and/or by degradation of the polymer. Examples of commercially available devices for sustained

drug release are Gliadel (implantable polyanhydride wafers that release drug at a constant rate as

the polymer degrades) for treating brain cancer, Lupron Depot (injectable polymer microspheres) for

endometriosis and prostate cancer, and Nutropin Depot for pituitary dwarfism. Lupron Depot

(leuprolide) is a suspension of microspheres made of poly-lactic-glycolic acid (PLGA) polymer

containing leuprolide acetate, and is administered once a month as an intramuscular injection. The

drug is slowly released into the blood to maintain a steady plasma concentration of leuprolide for

one month.

Since controlled-release drug-loaded microparticles are usually injected into the body, thereby

bypassing intestinal digestion and first-pass liver metabolism, one must ensure that the polymeric

drug carrier material is non-toxic, degrades within a reasonable time frame, and is excreted from the

body without any accumulation in the tissues. Biodegradable polymers gradually dissolve in body

fluids either due to enzymatic cleavage of polymer chains or hydrolytic breakdown (hydrolysis) of the

polymer. PLGA is an example of a biocompatible polymer (polyester) that is completely biodegrad-

able. As sections of the polymer chains in the outer layers of the drug–polymeric system are cleaved

and undergo dissolution into the surrounding medium, drug located in the interior becomes exposed

to the medium. The surrounding fluid may penetrate into the polymer mass, thereby expediting drug

diffusion out of the particle. Transport of drug molecules from within the polymeric mass to the

surrounding tissues or fluid is governed by several mechanisms that occur either serially or in

parallel. Mathematical formulation of the diffusional processes that govern controlled release is far

from trivial.

Here, we avoid dealing with the par tial differential equations used to model the diffusion of drug in

time and space within the particle and in the surrounding medium. We make the simplifying assumption

466

Numerical integration of ODEs