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

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

5.5 Newton’s method

Newton’s method, also called the Newton–Raphson method, is an open interval root-

ﬁnding method. It is a popular choice because of its second-order or quadratic rate of

convergence (r =2inEquation (5.10)). Newton’s method is applicable only when both

f (x)andf

(x) are continuous and differentiable. Only a single initial guess x

is required

to begin the root-ﬁnding algorithm. The x-intercept of the tangent to f (x)atx

yields

the next approximation, x

, to the root. This process is repeated until the solution

converges to within tolerable error. Unlike the bracketing methods discussed in

Sections 5.2 and 5.3, an open interval method does not always converge to a solution.

The solution technique is derived from the Taylor series expansion of f (x) about

an initial guess point x

located close to the root:

þ ΔxðÞ¼fx

ðÞþf

ðÞΔx þ

ζðÞΔx

; ζ E x

; x

þ Δx½: (5:17)

If x = x

+ Δx is the root of the function then fx

þ ΔxðÞ¼0 and the Tay lor series

expansion can be rewritten as follows:

0 ¼ fx

ðÞþf

ðÞx  x

ðÞþ0:5f

ζðÞx  x

ðÞ

: (5 :18)

If only the ﬁrst two terms on the right-hand side of the Taylor expansion series in

Equation (5.18) are retained, then

x ¼ x



ðÞ

: (5:19)

The validity of Equation (5.19) is contingent on our initial guess value x

being close

to the root enough that the error term squared (x – x

)

is much smaller than (x – x

Thus, the last term will be much smaller than the ﬁrst two terms in Equation (5.18),

which justiﬁes neglecting the f ″ term in the series expansion. Equation (5.19) is a

ﬁrst-order approximation of the function as a straight line through point (x

, f (x

))

with slope equal to f

). Because f (x) is not linear in x, x obtained from Equation

(5.19) will not be the root, but an approximation.

Equation (5.19) can be graphically interpreted as follows. If the function f (x)is

approximated by the tangent to the curve at an initial guess value of the root x

, then

the x-intercept of this tangent line produces the next approximation x

of the actual

root (Figure 5.9). The slope of this tangent line is calculated as

ðÞ¼

ðÞ0

 x

;

which, on rearrangement, yields Equation (5.19) where x,orx

, is the ﬁrst approx-

imation to the root. Equation (5.19) can be generalized as Equation (5.20).

Newton’s method

iþ1

¼ x



fðx

ðx

: (5:20)

Equation (5.20) deﬁnes the step-wise iterative technique of Newton’s method.

Figure 5.9 demonstrates graphically the step-wise algorithm to approach the root.

The algorithm for calculating a zero using Newton’s method is as follows.

327

5.5 Newton’s method

(1) Provide an initial guess x

of the root of f (x) = 0 either by plotting the function or by

performing an initial bisection or regula-falsi root-approximation step.

(2) Derive the analytical function f

(x).

(3) Calculate f (x

) and f

)atx

. For the ﬁrst iteration, x

= x

(4) Evaluate x

iþ1

¼ x

 fx

ðÞ=f

ðÞ.

(5) Repeat steps 3 and 4 until the method has converged upon a solution.

Example 5.2

Let’s use Newton’s method to solve the simple problem of determining the square root of 8, i.e.

fxðÞ¼x

 8 ¼ 0:

Then, f

xðÞ¼2x:

The iterative formula for this problem obtained from Equation (5.20) is

iþ1

¼ x



 8

Let’s start with an initial guess value x

= 3. Then

¼ 3 

 8

2  3

;



 8

2 

¼ 2:828431373;

¼ 2:828431373 

2:828431373

 8

2  2:828431373

¼ 2:828427125:

The actual solution is

ﬃﬃﬃ

¼ 2:828427125.

The error at each stage of calculation is given by

¼0:1715729;

¼0:0049062;

¼4:24780  10

6

;

¼3:18945  10

12

Figure 5.9

Illustration of Newton’s method.

f(x)

slope

f′(x

)

Actual root

slope

f′(x

)

328

Root-finding techniques for nonlinear equations

If Newton’s method is quadratically convergent then jε

iþ1



jε

for r = 2 should yield a number C that

converges to a constant value as i increases; i.e.

jε

¼ 0:1667;

jε

¼ 0:1765;

jε

¼ 0:1768:

Note the fast rate of convergence when r =2.

The last term on the right-hand side of Equation (5.18) provides valuable informa-

tion regarding the rate of convergence of Newton’s method. Equation (5.18) may be

rewritten a s follows:



¼ x



ðÞ



ζðÞ

ðÞ



 x

ðÞ

: (5:21)

After substituting Equation (5.20) into Equation (5.21) we obtain



 x

iþ1

ðÞ¼

ζðÞ

ðÞ



 x

ðÞ

j x



 x

iþ1

ðÞj

j x



 x

ðÞj

iþ1

ζðÞj

2jf

ðÞj

: (5:22)

Equation (5.22) proves that the order of convergence is r = 2 for Newton’s method.

This does not hold true if f

(x) = 0. Here

C ¼

ζðÞj

2jf

ðÞj

and

lim

i ! ∞

! x

root

C ¼



ðÞj

2jf



ðÞj

For Example 5.2, x

 8 ¼ 0, f

ﬃﬃﬃ



¼ 2; f

ﬃﬃﬃ



¼ 5:65685, and C = 0.1768.

5.5.1 Convergence issues

Convergence is not guaranteed

Newton’s method converges if the initial approximation of the solution is sufﬁciently

close to the root. How close the initial guess value must be to the root to guarantee

convergence depends upon the shape of the function, i.e. the values of f

(x) and f

(x)

at the guess values, at subsequent root approximations, and at the actual root. If the

value of f

(x) becomes zero during the iterative procedure, then the solution

329

5.5 Newton’s method

technique will diverg e. In fact, a near-zero slope or inﬂection point anywhere

between the initial guess and the root can cause persistent oscillations or divergence

of the solution. Figure 5.10(a) and (b) graphically depict two situations where

convergence cannot be attained.

Note that Equation (5.20) is of the form x = g(x), where

gxðÞ¼x 

fxðÞ

xðÞ

: (5:23)

It is necessary that g

xðÞ

1 near the root to establish convergence in any x = g(x)

ﬁxed-point iterative scheme (see Section 5.4), where

xðÞ¼

fxðÞf

xðÞ

ðf

xðÞÞ

: (5:24)

At the root, f (x) = 0 and g

(x) = 0. Therefore, there must be a range of values of x

near the root such that |g

(x)| < 1. Convergence of the method is certain as long as

the initial guess lies within this range.

Convergence rate is not always quadratic

The convergence rate was shown to be quadratic for this method; however, this

applies only at simple roots, i.e. roots of order p =1.

A root x



of f(x)issaidtobeoforder p if the function’s derivatives are continuous and the following hold

true:

fðx



Þ¼0; f

ðx



Þ¼0; f

ðx



Þ¼0;  f

ðp1Þ

ðx



Þ¼0; f

ðpÞ

ðx



Þ 6¼ 0:

For root x



of order p =1,

fðx



Þ¼0; f

ðx



Þ 6¼ 0; f

ðx



Þ 6¼ 0;

and f(x) can be written as x  x



ðÞhxðÞ, where hxðÞ6¼ 0 when x = x



For root x



of order p =2,

fðx



Þ¼0; f

ðx



Þ¼0; f

ðx



Þ 6¼ 0;

and f(x) can be written as x  x



ðÞ

hxðÞ.

Figure 5.10

Two situations in which Newton’s method does not converge.

f(x)

Location

of root

Divergence

Location

of root

Cyclic

behavior

(a) (b)

f(x)

330

Root-finding techniques for nonlinear equations

For simple roots (p = 1), Equation (5.24) yields g

xðÞ¼0, and from ﬁxed-

point iteration theory New ton’s method will have second-order convergence.

For mult iple roots, i.e. p >1, g

xðÞtakes on a 0/0 form and converges to a

number 0 < g

xðÞ< 1, which signiﬁes linear convergence (see Section 5.4). If

fxðÞ¼x  x



ðÞ

hxðÞ, is substituted into Equation (5.24), where p is the multi-

plicity of the root, it can be shown that g

xðÞ¼ðp  1Þ=p. For a double root,

p = 2 and g

xðÞ¼1=2: Quadratic convergence can be enforced by modifying

Equation (5.23) as follows:

gxðÞ¼x  p

fxðÞ

xðÞ

; (5:25)

where p is the order of the root x*.

Box 5.2A Equilibrium binding of multivalent ligand in solution to cell surfaces

Macromolecules that bind to receptors located on the surface of cells are called ligands.

Receptors are spec ial proteins, like enzymes, that have specific binding properties and usually

perform a function, such as transport of a ligand across a cell membrane or activation of a

signal to turn on certain genes, when specific ligands are bound. When a ligand contains

multiple binding sites such that it can bind to more than one receptor at any time, the ligand is

said t o crosslink or aggregate multiple receptors. When ligands such as antibody–antigen

complexes, ho rmones, or other extracellular signaling molec ules induce the aggregatio n of

receptors on a cell surface, various important cellular responses are triggered, such as the

expression of new protein.

A model for the equilibrium binding of multivalent ligand to cell surface receptors was proposed by

Perelson (1981). In this model, the multivalent ligand is assumed to have a total of ν active sites

available for binding to receptors on the surfaces of cells suspended in solution. A ligand molecule in

solution can initially attach to a cell surface by forming a bond between any of its ν binding sites and a

receptor on the cell surface (see Figure 5.11). However, once a single bond has formed between the

ligand and a cell receptor, the ligand orientation may allow only a limited number of binding sites to

interact with that cell. The total number of binding sites that a ligand may present to a cell surface for

binding is f.

Figure 5.11

Binding of a multivalent ligand in solution to monovalent receptors on a cell surface.

Cell

Cell surface

receptor

Multivalent ligand

ν ligand

binding sites

Extracellular space

Only f out of

ν binding sites

are presented by a ligand

to a cell at a time

Intracellular space

331

5.5 Newton’s method

The extent of binding of a ligand containing multiple identical binding sites to receptors on a cell

surface can be measured by determining the concentration of unbound receptors present on the cell

surface at equilibrium, which is given by the following equation (Lauffenburger and Linderman,

1993):

¼ R

1 þ ν



ð1 þ K

f1



; (5:26)

where

ν is the total number of binding sites present on the ligand surface,

f is the total number of binding sites available for binding to a single cell,

is the two-dimensional crosslinking equilibrium constant (1/(# ligands/ cell)),

is the ligand concentration in solution (M),

is the three-dimensional dissociation constant for binding of ligand molecules in solution to

cells (M),

is the equilibrium concentration of unbound receptors present on the cell surface (# receptors/

cell),

is the total number of receptors present on the platelet surface.

V on Willebrand factor, a plasma protein, binds platelet cells via multiple receptor–ligand bonds. For the

von Willebrand factor–platelet system, the parameters are estimated to be as follows (Mody and King,

2008):

(1) ν = 18,

(2) f =9,

(3) K

= 7.73 × 10

–5

(4) K

= 5.8 × 10

−5

cell/# ligands

(5) L

=2× 10

−9

(6) R

= 10 688 # receptors/cell.

We wish to determine the equilibrium concentration R

of unbound or free receptors on a platelet.

Replacing R

with x, we obtain

fxðÞ¼x 1 þ ν



ð1 þ xK

f1



 R

; (5:27)

xðÞ¼1 þ ν



1 þ xK

ðÞ

f1

þxK

ν f  1ðÞ



ð1 þ xK

f2

: (5:28)

A function called ReceptorLigandEqbmBinding was written in MATLAB to calculate f(x)

and f

(x) as given by Equations (5.27) and (5.28). A general function called Newtonsmethod was

created to solve a nonlinear problem in a single variable x using Newton’s method. When calling this

function at the command prompt, we specify the function ReceptorLigandEqbmBinding as

an input parameter to Newtonsmethod, which supplies the values of f(x) and f

(x) for successive

values of x. Note that the feval function requires the function name to be specified in string format,

and therefore the function name is included within single quotation marks when input as a parameter.

newtonsmethod(‘Receptor Ligand Eqbm Binding’,8000,0.5,0.5)

i x(i) f(x(i)) f’(x(i))

1 10521.7575 55.262449 1.084879

2 10470.8187 0.036442 1.083452

3 10470.7850 0.000000 1.083451

Note the rapid convergence of f(x) → 0. The number of receptors on the platelet surface bound to ligand

when equilibrium conditions are established is R

– 10 471 = 217 receptors or 2% of the total available

332

Root-finding techniques for nonlinear equations

receptors. The number of multivalent ligand molecules (containing multiple identical binding sites)

bound to receptors on a cell surface at equilibrium is given by the following equation:

i;eq

i!ðf  iÞ!



i1



; i ¼ 1; 2; ...; f;

where C

i,eq

is the number of ligand molecules bound to the cell by i bonds. Now that you have solved

Equation (5.26), determine the total number of ligand molecules

i¼1

that bind a platelet cell at

equilibrium. Use sum(ﬂoor(Ci)) to calculate the integral number of bound ligand molecules.

MATLAB program 5.6

function [fx, fxderiv] = ReceptorLigandEqbmBinding (Req)

% Equilibrium binding of multivalent ligand to cell surface

% Constants

RT = 10688; % Total number of receptors on a cell

L0 = 2.0e-9; % Ligand concentration in medium (M)

Kd = 7.73e-5; % Equilibrium dissociation constant (M)

v = 18; % Total number of binding sites on ligand

f = v/2; % Total number of binding sites that a single cell can bind

Kx = 5.8e-5; % Crosslinking equilibrium constant (cell/#)

fx = Req*(1+v*(L0/Kd)*(1+Kx*Req)^(f -1)) - RT;

fxderiv = 1+v*(L0/Kd)*(1+Kx*Req)^(f-1) + . . .

Req*v*(L0/Kd)*(f-1)*Kx*(1+Kx*Req)^(f-2);

MATLAB program 5.7

function newtonsmethod(func, x0, tolx, tolfx)

% Newton’s method used to solve a nonlinear equation in x

% Input variables

% func : nonlinear function

% x0 : initial guess value

% tolx : tolerance for error in estimating root

% tolfx: tolerance for error in function value at solution

% Other variables

maxloops = 50;

[fx, fxderiv] = feval(func,x0);

fprintf(‘ i x(i) f(x(i)) f’’(x(i)) \n’);

% Iterative solution scheme

for i = 1:maxloops

x1 = x0 - fx/fxderiv;

[fx, fxderiv] = feval(func,x1);

fprintf(‘%2d %5.4f %7.6f %7.6f \n’,i,x1,fx,fxderiv);

if (abs(x1 - x0) <=tolx && abs(fx) < tolfx)

break % Jump out of the for loop

end

x0 = x1;

end

333

5.5 Newton’s method

Box 5.3 Immobilized enzyme packed bed reactor

Enzymes are proteins in the body whose main functions are to catalyze a host of life-sustaining

chemical reactions, such as breaking down food into its primary components, e.g. proteins into amino

acids, converting glucose to CO

and H

O, synthesizing molecular components of cells, enabling

transcription and translation of genes, and so on. Enzymes have tremendous specificity for the substrate

acted upon and the reaction catalyzed. The rate of reaction r for an enzyme-catalyzed reaction such as

E þ S ! E þ P

where E is the enzyme, S is the substrate, and P is the product, is commonly described using

Michaelis–Menten kinetics:

R ¼

max

þ S

; (5:29)

where

S is the substrate concentration,

max

is the maximum reaction rate for a particular enzyme concentration E,

is the Michealis constant, equal to the substrate concentration at which the reaction rate is half of the

maximum rate.

Equation (5.29) can be integrated to yield an analytical relationship between the substrate concen-

tration and time as follows:



þ S

 S

ðÞ

¼ R

max

t: (5:30)

Enzymes can be produced (using genetic engineering techniques), purified, and then used ex vivo in a

suitably designed contacting device to catalyze similar reactions when placed in a conducive

physiological environment, such as at body temperature and favorable pH. Often enzymes are

immobilized on a support structure within the contacting equipment rather than kept free in solution.

This is advantageous for several reasons, such as (1) ease of separation of the contact stream from the

enzyme (the exiting stream must be cleared of enzyme if it is introduced into the body) and (2) reduced

loss of enzymes, which is especially important when the process is expensive.

In certain diseased conditions, the inefficient breakdown of a specific accumulating molecular entity

can pose a threat to the body since abnormally high concentration levels of the non-decomposed

substrate may be toxic to certain body tissues. In such cases the use of an immobilized enzyme reactor

to catalyze the breakdown reaction extracorporeally may be of therapeutic benefit. There are several

options available for the type of contacting device and contacting scheme used to carry out the reaction

process. A packed bed is a commonly used reactor configuration in which the catalyst particles are

dispersed on a support structure throughout the reactor and an aqueous stream containing the substrate

continuously passes through the reactor (McCabe et al., (2004) provides an introduction to packed bed

columns).

A prototype immobilized enzyme packed bed reactor has been designed whose function is to

degrade a metabolite X that is normally processed by the liver except in cases of liver disease or liver

malfunction. Accumulation of X in the body prevents functioning of certain cell types in the body.

Enzyme AZX is a patented enzyme that is immobilized within iner t pellet particles contained in the

reactor (Figure 5.12) and degrades X according to Michaelis–Menten kinetics (Equation (5.29)). A

simple mass balance over the reactor based on certain simplifying assumptions produces the following

design equation for the reactor (Fournier, 2007):

feed

out



þ KS

feed

 S

out

ðÞ¼

ηR

max

1  "ðÞL

; (5:31)

334

Root-finding techniques for nonlinear equations

where

feed

is the concentration of substrate in the feed stream,

out

is the concentration of substrate in the exiting stream,

K is the partition coefficient and is the ratio of the substrate concentration inside the pellet pores at the

surface to the substrate concentration in the external fluid at the pellet surface,

is the Michaelis constant (see Equation (5.29)),

is the reactor column cross-sectional area,

L is the reactor length,

ε is the porosity of the column (fractional volume not occupied by pellets),

max

is the maximum reaction velocity (see Equation (5.29)),

Q is the flowrate of the feed stream into the plug flow packed bed reactor,

η is the reaction efficiency that measures the effect of mass transfer resistance on reaction rate.

Suppose you are a bioengineer working for a company and you are asked to estimate the extent of

degradation of substrate X in a prototype reactor designed by the company. The parameters available to

you are as follows.

Packed bed column:

d = 2 cm,

L = 10 cm,

ε = 0.4,

Q = 80 ml/min,

K = 0.5,

feed

= 0.073 mg/ml.

Reaction kinetics:

= 0.022 mg/ml,

max

= 31 mg/(ml hr)

η = 0.6.

We rewrite Equation (5.31) as follows:

out

ðÞ¼K

feed

out



þ KS

feed

 S

out

ðÞ

ηR

max

1  "ðÞL

¼ 0(5:32)

and solve using Newton’s method. The derivative of fS

out

ðÞspecified in Equation (5.32) is given by

out

ðÞ¼

out

þ K



A MATLAB program packedbedreactor.m has been written to evaluate fS

out

ðÞand its derivative at a

supplied value of S

out

. Here, we run our function newtonsmethod using relevant arguments and an

initial guess of S

out

= 0.05:

Figure 5.12

Immobilized enzyme packed bed reactor.

Q, S

feed

Inlet or feed

stream

Pellet particle containing

immobilized enzymes

Exiting

stream

Q, S

out

335

5.5 Newton’s method

5.6 Secant method

The secant method is a combination of two root approximation techniques – the

regula-falsi method and Newton’s method – and accordingly has a convergence rate

that lies between linea r and quadratic. You may have already noted that f(x)in

Equation (5.3) of Example 5.1 has a complicated and messy derivative. In many

situations, the calculation of derivatives of a function may be prohibitively difﬁcult

or time-consuming. We seek an alternative to evaluating the analytical form of the

derivative for such situations. In the secant method, the analytical form of the

derivative of a function is substituted by an express ion for the slope of the secant

line connecting two points (x

and x

) that lie on f (x) close to the root, i.e.

xðÞ¼

ðÞfðx

ðx

 x

: (5:33)

If two initial estimates x

and x

of the root are speciﬁed, the x-intercept of the

line that joins the points (x

, f(x

)) and (x

, f(x

)) provides the next improved

approximation, x

, of the root being sought. Equation (5.12), reproduced below, is

newtonsmethod(‘packedbedreactor’,0.05,0.0001,0.001)

i x(i) f(x(i)) f’(x(i))

1 0.032239 0.001840 -1.182405

2 0.033795 0.000025 -1.150986

3 0.033816 0.000000 -1.150571

The extent of conversion of the substrate is (0.073 – 0.0338)/0.073 = 0.537 or 53.7%.

MATLAB program 5.8

function [f fderiv] = packedbedreactor(sout)

% Calculation of mass balance equation for immobilized enzyme ﬁxed bed

% reactor

% Input Variable

% sout : substrate concentration in outlet stream

% Constants

sfeed = 0.073; % mg/ml : substrate concentration at inlet

q = 80; % ml/min : ﬂowrate

dia = 2; % cm : diameter of reactor tube

Ac = pi/4*dia^2; % cross-sectional area of reactor

porosity = 0.4; % fraction of unoccupied volume in reactor

l = 10; % cm : length of reactor

K = 0.5; % partition coefﬁcient

vmax = 31/60; % mg/ml/min : maximum reaction velocity

Km = 0.022; % mg/ml : Michaelis constant

eff = 0.6; % reaction efﬁciency

f = Km*log(sfeed/sout) + K*(sfeed - sout) - ...

(1-porosity)*Ac*l*K*vmax*eff/q;

fderiv = -Km/sout - K;

336

Root-finding techniques for nonlinear equations