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

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

severity of motor symptoms is different among patients receiving neurostimulation

versus patients without medication.

4.15. Efﬁcacy of antivenom for severe cases of neurotoxicity caused by scorpion

stings Some varieties of scorpions produce neurotoxic venom that causes muscle

hyperactivity, respiration difﬁculties, and visual problems when injected into the

human body. A scorpion sting often produces severe illness among children but

milder symptoms in adults. Treatment of severe cases requires heavy sedation to

control violent involuntary movements, and assisted ventilation. A double-blinded

study was conducted to test the efﬁcacy of a scorpion antivenom (Boyer et al., 2009).

Eight children were randomly assigned to the treatment group that received the

antivenom dug, while seven children were randomly placed in the placebo group.

The placebo group received a higher sedative dose than the treatment group. At

baseline, the mean venom concentrations in blood were comparable in both groups.

The plasma levels of venom were measured at 1 hour from treatment time, (values

are close to, but are not exactly, the true measurements). These are presented in

Table P4.8. At 1 hour from treatment, are the median plasma venom levels of the

two groups signiﬁcantly different? Construct the null hypothesis, choose the appro-

priate hypothesis test, and state your conclusions based on the value of the calculated

test statistic.

4.16. Treatment of familial adenomatous polyposis Formation of adenomatous polyps

in the colon is often followed by development of colon cancer. Non-steroidal anti-

inﬂammatory drugs (NSAIDs), such as aspirin, are cyclooxygenase (cox) inhibitors

and have been linked in multiple studies to regression of colorectal adenomas and a

lowering of colon cancer incidence. Mutation of the adenomatous polyposis coli

(APC) gene predisposes a person to colon adenoma. Those who inherit the mutated

APC gen e develop familial adenomatous polyps and are at greater risk for develop-

ing colon cancer.

In a randomized, double-blinded study (Steinbach et al., 2000), a cox-2 inhibitor

called celecoxib was administered to patients who had developed familial adenom-

atous polyps. Placebo was given to another group of patients as a control case. The

percent changes from baseline after 6 months of treatment in the number of color-

ectal polyps in patients treated with celecoxib twice per day and with placebo are

listed below.

Table P4.8.

Plasma venom levels (ng/ml)

Antivenom treatment Placebo treatment

0.007 1.212

0.001 0.026

0.022 8.676

0.032 0.334

0.003 3.634

0.041 0.885

0.009 2.419

1.795

307

4.12 Problems

Placebo (n = 15)

20, 17, 14, 0, 0, 0, 0, 0, −2, −8, −10, −13, −17, −20, −45

400 mg of celecoxib administered 2 × daily (n = 30)

20, 4, 0, 0, 0, 0, 0, −10, −14, −16, −16.5, −18, −19, −19.5, −28.5, −30, −30.5, −38, −38,

−40, −40, −43, −45, −47, −50, −52, −55, −58, −70, −80.

Using a non-parametric method, test the hypothesis that celecoxib reduces the

occurrence of adenomatous polyps in patients with familial adenomatous polyposis.

What assumptions are necessary to apply the appropriate hypothesis test?

References

Bassett, M. T., Dumanovsky, T., Huang, C., Silver, L. D., Young, C., Nonas, C.,

Matte, T. D., Chideya, S. and Frieden, T. R. (2008). Purchasing Behavior and Calorie

Information at Fast-Food Chains in New York City, 2007. Am J Public Health , 98,

1457–9.

Bedair, H. S., Karthikeyan, T., Quintero, A., Li, Y., and Huard, J. (2008) Angiotensin II

Receptor Blockade Administered after Injury Improves Muscle Regeneration and

Decreases Fibrosis in Normal Skeletal Muscle. Am. J. Sports Med., 36, 1548–54.

Boyer, L. V., Theodorou, A. A., Berg, R. A., Mallie, J., Chavez-Mendez, A., Garcia-

Ubbelohde, W., Hardiman, S., and Alagon, A. (2009) Antivenom for Critically Ill

Children with Neurotoxicity from Scorpion Stings. N. Engl. J. Med., 360, 2090–8.

Brooke, B. S., Habashi, J. P., Judge, D. P., Patel, N., Loeys, B., and Dietz, H. C., 3rd. (2008)

Angiotensin Ii Blockade and Aortic-Root Dilation in Marfan’s Syndrome. N. Engl.

J. Med., 358, 2787–95.

Czeisler, C. A., Walsh, J. K., Roth, T. et al. (2005) Modaﬁnil for Excessive Sleepiness

Associated with Shift-Work Sleep Disorder. N. Engl. J. Med., 353, 476–86.

Deuschl, G., Schade-Brittinger, C., Krack, P. et al. (2006) A Randomized Trial of Deep-Brain

Stimulation for Parkinson’s Disease. N. Engl. J. Med., 355, 896–908.

Doll, R. and Kinlen, L. (1977) Fluoridation of Water and Cancer Mortality in the U.S.A.

Lancet, 1, 1300–2.

Grocott, M. P., Martin, D. S., Levett, D. Z., McMorrow, R., Windsor, J., and

Montgomery, H. E. (2009) Arterial Blood Gases and Oxygen Content in Climbers on

Mount Everest. N. Engl. J. Med., 360, 140–9.

Malizia, B. A., Hacker, M. R., and Penzias, A. S. (2009) Cumulative Live-Birth Rates After In

Vitro Fertilization. N. Engl. J. Med., 360, 236–43.

Moser, K., Stanﬁeld, K. M., and Leon, D. A. (2008) Birthweight and Gestational Age by

Ethnic Group, England and Wales 2005: Introducing New Data on Births. Health Stat.

Q., issue no. 39, 22–31, 34–55.

Mruk, J. S., Webster, M. W., Heras, M., Reid, J. M., Grill, D. E., and Chesebro, J. H.

(2000) Flavone-8-Acetic Acid (Flavonoid) Profoundly Reduces Platelet-Dependent

Thrombosis and Vasoconstriction after Deep Arterial Injury in Vivo. Circulation,

101, 324–8.

Newburger, J. W., Sleeper, L. A., McCrindle, B. W. et al. (2007) Randomized Trial of Pulsed

Corticosteroid Therapy for Primary Treatment of Kawasaki Disease. N. Engl. J. Med.,

356, 663–75.

Nooyens, A. C., van Gelder, B. M., and Verschuren, W. M. (2008) Smoking and Cognitive

Decline among Middle-Aged Men and Women: The Doetinchem Cohort Study. Am.

J. Public Health, 98, 2244–50.

Piot, C., Croisille, P., Staat, P., et al. (2008) Effect of Cyclosporine on Reperfusion Injury in

Acute Myocardial Infarction. N. Engl. J. Med., 359, 473–81.

308

Hypothesis testing

Samuels, M. L. and Witmer, J. A. (2003) Statistics for the Life Sciences (Upper Saddle River;

NJ: Pearson Education, Inc.).

Steinbach, G., Lynch, P. M., Phillips, R. K. et al. (2000) The Effect of Celecoxib, a

Cyclooxygenase-2 Inhibitor, in Familial Adenomatous Polyposis. N. Engl. J. Med.,

342, 1946–52.

Visentin, M., Stea, S., Squarzoni, S., Reggiani, M., Fagnano, C., Antonietti, B., and Toni, A.

(2005) Isolation and Characterization of Wear Debris Generated in Patients Wearing

Polyethylene Hylamer Inserts, Gamma Irradiated in Air. J. Biomater. Appl., 20, 103–21.

Whitney, C. G., Farley, M. M., Hadler, J. et al. (2000) Increasing Prevalence of Multidrug-

Resistant Streptococcus Pneumoniae in the United States. New Engl. J. Med., 343,

1917–24.

Yiamouyiannis, J. and Burk, D. (1977) Fluoridation and Cancer: Age-Dependence of Cancer

Mortality Related to Artiﬁcial Fluoridation. Fluoride, 10, 102–125.

309

References

5 Root-finding techniques for

nonlinear equations

5.1 Introduction

As previously discussed in Chapter 2, a system of linear equations can be expressed as

Ax = b, where A and b are not functions of x. On the other hand, a system of

equations written as Ax = b, in which A and/or b is a function of x, contains one

or more nonlinear equations. A polynomial in x of order  2 has nonlinear depend-

ency on the variable x and is an example of a nonlinear equation. Equations

involving trigonometric and/or transcedental functions such as sin(x), cos(x),

log(x), e

, and coth(x) are also nonlinear in x. The solution x of a nonlinear equation

f (x) = 0 is termed a root of the equation or a zero of the function f (x). Application of

numerical techniques to obtain the roots of nonlinear equations is the focus of this

chapter.

Nonlinear equations that are not amenable to an analytical solution may be

solved using iterative numerical schemes. A system of consistent linear equations

will always yield a solution. In fact, linear equations that are consistent but not

underdetermined will always have a unique solution, regardl ess of the technique

used to solve the system of linear equations. On the other hand, a numerical method

used to solve a nonlinear equation may not necessarily converge to a solution. Often,

certain conditions must be met before a particular numerical technique can success-

fully produce a solution. For example, convergence of a method may be contingent

on the magnitude of the derivative of a function lying within bounded limits, or may

depend on how close the initial guess value is to the root. In many cases, nonlinear

equations do not have unique solutions, i.e. an equation may admit several roots.

The root that is obtained will usually depend on the numerical method chosen as well

as the starting guess value.

Numerical schemes for solving nonlinear equations in one variable can be broadly

classiﬁed into two categories: bracketing methods and open interval methods.

Sections 5.2 and 5.3 discuss the bisection method and the regula-falsi method,

which are numerical root-ﬁnding methods of the ﬁrst category. Bracketing methods

of solution require that a closed interval on the number line containing the desired

root be speciﬁed at the beginning of the ﬁrst iteration. With each iteration, the

interval that brackets the root shrinks, allowing the location of the root to be

known with increasing accuracy. When employing bracketing methods, convergence

is usually guaranteed unless the initially speciﬁed bracket does not contain the roo t,

or contains a singularity point

instead of a zero. The latter can be avoided by

plotting the function graphically to study the behavior of the function and location

of the roots. Bracketing methods of solution are reliable but converge slowly. Faster

methods of convergence are usually preferred. This explains the popularity of open

A singularity point of a function f (x) is the x value at which the magnitude of the function blows up, i.e.

tends towards plus and/or minus inﬁnity.

Box 5.1A Rheological properties of blood

Blood is responsible for the convective transport of various vital compounds in the body such as oxygen,

nutrients, and proteins. Blood consists of plasma and formed elements, i.e. red blood cells (RBCs) or

erythrocytes, white blood cells (WBCs) or leukocytes, and platelets. The distinct size, shape, and

deformability of each of the three blood cell types – erythrocytes, leukocytes, and platelets strongly

influence the flow characteristics of these cells within the vasculature.

The volume fraction of RBCs in blood is called the hematocrit and typically ranges from 40–50% in

humans. RBCs have a biconcave disc shape with a diameter of 6–8 μm and thickness of 2 μm, and are

highly deformable. The distinctive shape of RBCs facilitates aggregation of red blood cells during flow.

Aggregates of RBCs flow along the centerline of blood vessels, pushing other blood cells (leukocytes

and platelets) to the wall periphery. Thus, flowing blood segregates into two layers: (1) an RBC-rich

core surrounded by (2) a concentric cell-depleted plasma layer near the wall. Figure 5.1 illustrates

this. The overall hematocrit of blood is less than the hematocrit exclusively observed in the central flow

(core) region where the RBC concentration is highest.

The tendency of RBCs to flow in the center of a vessel gives rise to a flow phenomenon, known as the

Fah raeus effect, in which the apparent hematocrit of blood that discharges from a vessel is gr eater than the

actual hematocrit of the flowing blood. This results from the gr eater velocity of the core region of the flowing

blood and therefore a greater concentration of RBCs in the collected discharge from the vessel.

The viscosity

of blood in the core region is much greater than that of the concentric plasma layer at

the periphery.

The viscosity of plasma is given by

plasma

¼ 1 :2cP

(cP = centiPoise, where 1 Poise (g/(cm s)) = 100 cP). The viscosity of the core region can be

calculated as follows:

core

¼ μ

plasma

1  αH

core

; (5:1)

where H

core

is the core region hematocrit,

α ¼ 0:070 exp 2:49H

core

1107

expð1:69H

core



;

Viscosity is a physical property of a ﬂuid that describes the ease with which a ﬂuid can be sheared. A

larger viscosity means that more energy is required to induce and maintain ﬂow.

Figure 5.1

Cross-sectional and longitudinal schematic view of blood flow in a blood vessel.

Blood flow

Blood vessel wall

Core region (RBC-rich)

Cell-depleted plasma layer (RBC-poor)

311

5.1 Introduction

interval methods such as Newton’s method and the secant method discussed in

Sections 5.5 and 5.6, which have higher orders of convergence but do not always

guarantee a solution.

Often, a biomedical engineer will face the challenge of solving engineering prob-

lems consisting of a system of two or more simultaneous nonlinear equations.

Newton’s method for solution of multiple nonlinear equations is discussed in

Section 5.7. MATLAB provides a built-in function called fzero that uses a hybrid

numerical method to determine the zeros of any nonlinear function. The MATLAB

fzero function to solve for the roots of an equation is discussed in Section 5.8.

5.2 Bisection method

The bisection method is perhaps the simplest method devised for nonlinear root

ﬁnding. The ﬁrst step in this numerical procedure involves determining an interval

[a, b] that contains the desired root x* of a nonlinear function f(x). The interval [a, b]

is certain to contain at least one zero of the function if

(1) fðxÞ is continuous within the interval, and

(2) faðÞfbðÞ

If the ﬁrst condition is true, then we know that the function has no discontinuities

(e.g. a singularity or a step function) within the bracketing interval and thus the

function must cross the x-axis if the second condition is also true.

How do you choose a bracketing interval(s)? Before tackling any nonlinear root-

ﬁnding problem, regardless of the numerical method of solution chosen, you should

ﬁrst plot the function f (x). The graphical plot supplies a host of information, such as

the number of roots, approxim ate location of each root, presence of any disconti-

nuities or singularities, behavior of the function in proximity of each root, and

and T is the temperature in Kelvin.

As an investigator of blood rheology, you want to calculate the hematocrit of the core region in a

blood vessel. You have obtained the following expression for the core hematocrit (Fournier, 2007):

core

inlet

¼ 1 þ

ð1  σ

21 σ

ðÞþσ

plasma

core

for H

core

0:6; (5:2)

where

σ ¼ 1 

;

δ is the thickness of the plasma layer, and R is the radius of the blood vessel.

If H

inlet

(hematocrit at inlet of blood vessel) = 0.44, δ =2.92μm, R =15μm, and T =310Karethe

estimated values, then σ is calculated as 0.8053. On substituting Equation (5.1) into (5.2), we obtain

core

¼ H

inlet

1 þ

ð1  σ

21 σ

ðÞþσ

ð1  αH

core

Þ½

: (5:3)

Nonlinear root-finding techniques are required to solve Equation (5.3).

312

Root-finding techniques for nonlinear equations

existence of multiple roots. A root that occurs more than once in the solution is

called a multiple root. For example, the equation

fxðÞ¼x

 6x

þ 9x  4 ¼ 0(5:4)

can be factored as

x  1ðÞx  1ðÞx  4ðÞ¼0:

This equation has a double root at x =1.

Once a plot of f(x) has revealed the approximate location of the desired root, we

can easily obtain the interval [a, b] encapsulating the root. Next, we determine the

value of the function at the endpoints of the interval. If the interval contains one

root, then f(a)·f(b) must be less than zero, i.e. f(x) evaluated at a and b will lie on

opposite sides of the x-axis. This may not be the case if the interval contains more

than one root. For example, on both sides of a double root, the value of the function

remains on the same side of the x-axis. Equation (5.4) is plotted in Figure 5.2 to

illustrate the behavior of a function at a double root.

The behavior of a function at a multiple root is discussed in greater detail in

Section 5.5. The bisection method will fail to uncover a root if the function does not

cross the x-axis at the root, and thus cannot be used to solve for a double root.

The bisection method searches for the root by depleting the size of the bracketing

interval by half so that the location of the root is known with an accuracy that

doubles with every iteration. If our initial interval is [a, b], then we proceed as follows

(Figure 5.3 illustrates these steps).

(1) Bisect the inter val such that the midpoint of the interval is x

¼ða þbÞ=2: Unless x

is precisely the root, i.e. unless f(x

) = 0, the root must lie either in the interval [a, x

]

or in [ x

, b].

(2) To determine which interval contains the root, f(x

) is calculated and then compared

with f(a) and f(b).

Figure 5.2

Plot of a function that contains a double root. At the double root, the function does not cross the x-axis, but is tangent

to the axis.

−1 0 1 2 3 4 5

−20

−10

f(x) = x

− 6x

+ 9x − 4

f(x)

Double root at x = 1

Single root at x = 4

313

5.2 Bisection method

(a) If fðaÞfðx

0, then the root is located in the interval [a, x

], and the current

value of b is replaced by x

, i.e. b = x

. The new interval is half the size of the

previous starting interval.

(b) If fðaÞfðx

0, the root lies within [x

, b], and a is set equal to x

(3) The bisect ion step is repeated to ﬁnd the midpoint, x

, of the new bracketing interval,

which is eithe r [a, x

]or[x

, b]. Then, the next smaller interv al that encloses the root

is selected. In Figure 5.3, the interval [x

, x

] is found to bracket the root. This new

interval is one-quarter the size of the original interval [a, b].

(4) The iterative process is stopped when the algorithm has converged upon a solution.

The criteria for attaining convergence are discussed below.

We must specify certain criteria that establish whether we have arrived at a

solution. After every iteration, the solution should be tested for convergence to

determine whether to continue the bisection method or terminate the process.

Tests for convergence estimate the error in the numerical solution and check if this

is within the speciﬁed tolerance.

One method to test for convergence calculates the difference between the numeri-

cally calculated root and the previous root estimate. If the numerical procedure is

convergent such that |x

n+1

– x

|<|x

– x

n–1

| for all n, then the difference between

successive estimates of the root indicates whether the solution has converged. The

convergence criterion may be speciﬁed as follows:

n þ1

 x

TolðxÞðabsolut e differenceÞ (5:5)

n þ1

 x

nþ1

TolðxÞðrelative differenceÞ: (5:6)

Here, Tol(

) is the maximum allowable error tolerance for the root or the

function whose zeros we seek. For the bisection method, there is an alternative

Figure 5.3

Illustration of the bisection method.

f(x)

f(b)

f(x

)

f(a)

f(x

)

3rd interval

]

2nd interval [a x

]

1st interval [a b]

314

Root-finding techniques for nonlinear equations

way of calculating the error in the root estimation. The bracketing interval

speciﬁes the maximum error with which the actual root is known. Within the

bracketing interval [a, b] the actual root lies within a distance of jb  aj=2 from

the midpoint of the interval. At the ﬁrst iteration i = 1, the midpoint of the

interval ða þ bÞ=2 is our approximated root, and the maximum possible error in

root estimation is jb  aj=2, in which case the root lies at either of the interval

endpoints. Thus, the magnitude of the error E

¼ x

 x



is between 0 and

jb  aj=2.

At the second iteration, the interval size jb  aj=2 is again halved. For

i ¼ 2 the error bounds in root estimation is given by

jb  aj

;

and for

i ¼ 3 the error bounds in root estimation is given by

b  a

;

i ¼ n the error bounds in root estimation is given by

b  a

At every iteration, the error in root estimation reduces linearly by half.

If a tolerance for the error in estimation of the root is speciﬁed, we can deﬁne a

convergence criterion

jb  aj

 TolðxÞ: (5:7)

The number of iterations n that will result in a bracketing interval small enough to

permit tolerable error is calculated by taking the log of both sides of Equation (5.7),

so that

n 

log jb  aj=TolðxÞðÞ

log 2

: (5:8)

The rate of convergence determines how quickly a numerical procedure approaches

the solution or root of the equation, or, in other words, how many iterations it takes

to converge upon a solution. The rate of convergence differs from one numerical

root-ﬁnding method to another.

If the function has a gradual slope at the location of the root, then Equations

(5.5),(5.6), or (5.8) are usually sufﬁcient to establish a stop criterion. However, if the

function has a steep slope when it crosses the x-axis, then even small deviations from

the actual root can lead to substantial de partures of the functio n value from zero. In

such cases, another important criterion for convergence speciﬁes the tolerance for

f(x), i.e.

fxðÞ

TolðfxðÞÞ; (5:9)

TolðfxðÞÞmay be speciﬁed as a fraction of the function values at the endpoints of the

initial interval [a, b].

315

5.2 Bisection method

The rate of convergence or order of convergence r for a numerical technique is defined as

lim

i!∞

iþ1

¼ C; (5:10)

where

iþ1

is the error in root estimation for the (i + 1)th iteration,

is the error in root estimation for the ith iteration; E

¼ x

 x



, where x* is the root.

The efficiency of a root-finding method depends on the values of C and r. It is desirable to have a large

value of r and a small value for C. For the bisection method, C = 0.5 and r =1.

Box 5.1B Rheological properties of blood

As described in Box 5.1A, we would like to calculate the core region hematocrit (the volume fraction of

red blood cells within the core region). Equation (5.3) is rearranged to yield

core

ðÞ

¼ 1 þ

ð1  σ

21 σ

ðÞþσ

ð1  αH

core

Þ½



core

inlet

;

where α is given by Equation (5.1).Let’s denote H

core

by x. On combining the above equation with

Equation (5.1) and substituting the known constants into the resulting equation, we obtain

fxðÞ¼1 þ

ð1  σ

21 σ

ðÞþσ

ð1  0:070 exp 2:49x þ

1107

expð1:69xÞ



xÞ





0:44

;

(5:11)

and f(x) is plotted in Figure 5.4 to determine the behavior of the function, the number of roots, and the

location of the desired root. We seek the value of x at f(x)=0.

The plot in Figure 5.4 shows that the function is well behaved near the root. Equation (5.3) is only

valid for H

core

≤ 0.6, and we do not expect H

core

< H

inlet

. Hence we choose the plotting interval 0.4 < x

< 0.6. On inspection of Figure 5.4 we choose the bracketing interval [0.5, 0.55].

Figure 5.4

Plot of the function specified by Equation (5.11).

0.4 0.45 0.5 0.55 0.6

−0.2

−0.1

0.1

0.2

0.3

f(x)

Location of

desired root

316

Root-finding techniques for nonlinear equations