Yang X.-S. Introductory Mathematics for Earth Scientists

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

200 Chapter 14. Numerical Methods

= f (x

i+1/2

) +

i+1/2

). (14.17)

where O(h

) means that the value is about the order of h

, or

O(h

) = Kh

where K is a constant.

The trapezium rule is just one of the simple and popular schemes

for numerical integration with the error of O(h

). If we want higher

accuracy, we can either reduce h or use a better approximation for f (x).

A small h means a large n, which implies that we have to do the sum

of many small sections, and it may increase the computational time.

On the other hand, we can use higher-order approximations for the

curve. Instead of using straight lines or linear approximations for curve

segments, we can use parabolas or quadratic approximations. For any

consecutive three points x

i−1

, x

and x

i+1

, we can construct a parabola

in the form

f(x

+ t) = f

+ αt + βt

, t ∈ [−h, h]. (14.18)

As this parabola must go through the three k nown points (x

i−1

, f

i−1

)

at t = −h, (x

, f

) at t = 0 and x

i+1

, f

i+1

at t = h, we have the

following equations for α and β

i−1

= f

− αh + βh

, (14.19)

and

i+1

= f

+ αh + βh

, (14.20)

which lead to

α =

i+1

− f

i−1

, β =

i−1

− 2f

+ f

i+1

. (14.21)

In fact, α is the centred approximation for the ﬁrst derivative f

and β

is related to the central diﬀerence scheme for the second derivative f

Therefore, the integral from x

i−1

to x

i+1

can be approximated by

i+1

i−1

f(x)dx ≈

−h

+ αt + βt

]dt =

i−1

+ 4f

+ f

i+1

where we have substituted the expressions for α and β. To e nsure

the whole interval [a, b] ca n be divided up to form three-point approx-

imations without any point left out, n must be even. Therefore, the

estimate of the integral becomes

I ≈

+4(f

+...+ f

n−1

)+2(f

+...+ f

n−2

)+f

], (14.22)

which is the standard Simpson’s rule.

As the approximation for the function f (x) is quadratic, an order

higher than the linear form, the error estimate of Simpson’s rule is thus

14.4 Numerical Solutions of ODEs 201

O(h

) or O(h

000

) to be more speciﬁc. There a re many variations of

Simpson’s rule with higher-or der accuracies such as O(h

(4)

Example 14.3: We know the value of the integral

I = erf(1) =

√

−x

dx = 0.8427007929.

Let us now estimate

I =

√

J, J =

−x

dx,

using the Simpson rule with n = 8 and h = (1 − 0)/8 = 0.125. We have

J ≈

+ 4(f

+ f

) + 2(f

+ f

) + f

Since f

= e

−x

= e

−(i∗h)

, we h a ve f

= 1, f

= 0.984496, f

0.939413, f

= 0.868815, f

= 0.778801, f

= 0.676634, f

= 0.569783,

= 0.465043, and f

= 0.367879. Now the integral estimate of J is

J ≈

0.125

[1 + 4 × 2.9 949885 + 2 × 2.2879 967 + 0.367879]

≈

0.125

× 17.92 3827 ≈ 0.746826.

Finally, the integral estimate of I is

I =

√

J =

√

3.1415926

× 0.746 826 ≈ 0.84270 3.

We can see that this estimate is acc u rate to the 5th decimal place.

There are other even better ways for evaluating the integral more

accurately using fewer points of evaluation. Such numerical integration

is called the Gaus sian integration or Gaussian quadrature. Interested

readers can refer to more advanced books on numerical methods.

14.4 Numerical Solutions of ODEs

The simplest ﬁrst-order diﬀerential equation can be written as

= f(x, y), (14.23)

where f(x, y) is a known function of x and y. In principle, the solution

can be obtained by direct integration

y(x) = y(x

) +

f(x, y(x))dx, (14.24 )

202 Chapter 14. Numerical Methods

but in practice it is usually impossible to do the integration analytically,

as it requires the solution y(x) to evaluate the right-hand s ide. Numer-

ical integration in this case is the most common technique for obtaining

approximate solutions. A naive approach is to use the standard tech-

niques such as Simpson’s rule for numerical integration to evaluate the

integral numerically. However, since we have to use some approxima-

tion for y(x), such techniques hardly work here without modiﬁcation.

A better approach is to start from the known initial value y(x

) at x

and try to march to the next point at x = x

+ h where h is a small

increment. In this manner, the solution at all other values of x can be

estimated numerically. This is essentially the Euler scheme.

14.4.1 Euler Scheme

The basic idea of the Euler method is to approximate the integral

(14.24) using an iterative procedure

n+1

= y

n+1

f(x, y)dx = y

+ hf (x

, y

), (14.25)

where h = ∆x = x

n+1

− x

is a small increment. Here we use the

notations x

= x

+ nh(n = 0, 1, 2, ..., N ) and y

= y(x

). This is

essentially to divide the interval from x

to x = Nh into N small

interval of width h, so the value of f(x, y) is approximated by its value

in the interval as f (x

, y

The above equation can also be viewed as the approximation of the

ﬁrst derivative

= y

n+1

− y

. (14.26)

This is a forward diﬀerence scheme for approximating the ﬁrst deriva-

tive. From the diﬀerentiation and integration discussed earlier in this

book, we know that such an approximation is very crude, and the accu-

racy is at most O(h

). A potentially better way is to us e the so-called

central diﬀerence to approximate the ﬁrst derivative

n+1

− y

n−1

, (14.27)

which uses two steps. This does not always solve the problem. In order

to get a reasonable accuracy, we have to use a very small value of h.

There is also an issue called numerical stability. If the step size h is to o

large, there is a possibility that any error will be ampliﬁed during each

iteration, and after some iterations, the value y

might be unstable and

become meaningless. A fa r better and more stable numerical scheme

is the well-known Runge-Kutta method.

14.4 Numerical Solutions of ODEs 203

14.4.2 Runge-Kutta Method

The Runge-Kutta method uses a trial step to march the solution to the

midpoint of the interval by the central diﬀerence

n+1/2

= y

f(x

, y

), (14.28)

then it combines with the standard forward diﬀerence as used in Euler

scheme

n+1

= y

+ hf(x

n+1/2

, y

n+1/2

). (14.29)

This scheme can be viewed as a predictor -corrector method a s it ﬁrst

predicts a value in the midpoint, and then corrects it to an even better

estimate at the next step.

If we follow the same idea, we can devise hig her-order methods with

multiple steps, which will give higher accuracy. The popular classical

Runge-Kutta method is a fourth-order method, which involves four

steps to advance from n to n + 1. We have

a = hf(x

, y

), b = hf(x

+ h/2, y

+ a/2),

c = hf(x

+ h/2, y

+ b/2), d = hf(x

+ h, y

+ c),

n+1

= y

a + 2(b + c) + d

, (14.30)

which is fourth-order accura te.

Example 14.4: Let us solve the nonlinear equ a tion numerically

+ y

= −1, x ∈ [0, 2]

with the initial c ondition y(0) = 1. We know that

= −(1 + y

), or

−

1 + y

dy = −tan

−1

y =

dx = x + K,

where K is the constant of the integration. This gives

y = −tan(x + K).

Using the initial condition y = 1 and x = 0, we have 1 = −tan(K), or

K = −tan

−1

(1) = −

The analytical solution becomes

y(x) = −tan(x −

204 Chapter 14. Numerical Methods

On the interval [0, 2], let us ﬁrst solve the equation using the Euler

scheme for h = 0.5. There are ﬁve points x

= ih (where i = 0, 1, 2, 3, 4).

As dy/dx = f(y) = −1 − y

, we have the Euler scheme

n+1

= y

+ hf(y

) = y

− h − hy

From the initial condition y

= 1, we now have

= y

− h − hy

= 1 − 0.5 − 0.5 × 1

= 0,

≈ −0.5, y

≈ −1.125, y

= −2.2578.

These are signiﬁcantly diﬀerent (about 30%) from the exact solution

∗ = 1, y

∗

= 0.2934079, y

∗

= −0.21795809,

∗

= −0.86756212, y

∗

= −2.68770693.

Now let us use the Ru n ge-Kutta method to so lve the same equation

to see if it is better. Since f(x

, y

) = −1 − y

, we have

a = hf(x

, y

) = −h(1 + y

), b = −h[1 + (y

)

c = −h[1 + (y

)

], d = −h[1 + (y

+ c)

and

n+1

= y

a + 2(b + c) + d

From y

= 1 and h = 0.5, we have

≈ 0.29043, y

≈ −0.22062, y

= −0.87185, y

≈ −2.67667.

These values are wi th in about 1% of the an a lytical solution y

∗

. We can

see that even with the same step size h, the Runge-K u tta method is more

eﬃcient and accurate than the Euler sche m e .

Generally speaking, higher-order schemes are better than lower-

order schemes, but not always. Let us look at an example.

Example 14 .5: We know from the example in an earlier c h a p te r on

ODEs that the solution to the diﬀerential e q u a tion

= (x

+ e

−y

with the initial c ondition y(0) = 0 at x

= 0 is

y = ln(e

14.4 Numerical Solutions of ODEs 205

Now let us try to solve it numerically to estimate y(1).

First, we divide the interval [0, 1 ] into two steps so that h = (1−0)/2 =

0.5. We ha ve x

= 0 and y

= 0. Then, we use the Euler scheme, and we

have y

at x

= 0.5

= y

+hf(x

, y

) = y

+h(x

−y

= 0+0.5×(0

−0

= 0.5.

The next step is to evaluate y

at x

= 2h = 1, and we have

= y

+ h(x

+ e

−y

= 0.5 + 0.5 × (0.5

+ e

0.5

−0.5

≈ 1.0189.

We know the exact value at x = 1 is

∗

= y(1 ) = ln(e

) = 1.07099503.

We can see that the value y

obtained by the Euler scheme is about 5%

from the true solution. Now let us see whether the Runge-Kutta method

can do better.

For the ﬁrst step to get y

, we have

a = hf(x

, y

) = h(x

+ e

)

−y

= 0.5 × (0

+ e

) × e

−0

= 0.5,

b= hf(x

+h/2, y

+a/2)=0.5×[(0+h/2)

0+h/2

]×e

−(0+a/2)

=0.50152,

and

c = hf(x

+h/2, y

+b/2) ≈ 0.5 0114, d = hf (x

+h, y

+c) ≈ 0.518363.

Now we get

= y

a + 2(b + c) + d

= 0.503947.

The next ste p to get y

from x

= 0.5 and y

= 0.503947, we have

a = hf(x

, y

) = 0.5 × (0.5

+ e

0.5

) × e

−0.503947

≈ 0.51691,

b= hf(x

, y

)≈0.56765, c =hf (x

, y

)≈0.55343,

and d = hf(x

+ h, y

+ c) ≈ 0.645 80. Finally, we have

= y

a + 2(b + c) + d

≈ 1.07142.

We can see that Runge-Kutta is m ore accu rate than the Euler method,

and the solution is about 0.04% from the true value. However, this Runge-

Kutta method is of O(h

) accurate.

This example demonstrates that higher-order methods are usually

the best choice, and they generally work better for most problems.

206 Chapter 14. Numerical Methods

That is why they are more widely used and have been implemented in

many software packages. The g e neral guidance is that you should try to

use these methods ﬁrst until you ﬁnd a better method for a particular

problem.

So far, we have only introduced the numerical method to the ﬁrst-

order equations . Higher-order or dinary diﬀerential equations can al-

ways be converted to a ﬁrst-or der system. For example, the following

second-order equation

(x) + p(x)y

(x) + q(x)y(x) = r(x), (14.31)

can be rewritten as a system if we let u = y

(x)

= u, u

= r(x) − q(x)y(x) −p(x)u. (14.32)

The above system of two equations can be written as







r(x) −q(x)y − p(x)u



. (14.33)

This is essentially the same as (14.23) in the context of vectors and

matrices, and it can essentially be solved using the similar method

with little modiﬁcatio ns.

Numerical solutions of partial diﬀerential equations are even more

diﬃcult, especia lly when the domain is multidimensional and irregular.

There are various metho ds for solving PDEs, including ﬁnite diﬀer-

ent methods, ﬁnite volume methods, ﬁnite element methods, boundar y

element methods, spectral methods, and meshless methods.

The literature about numerical methods is vast, especially about

the solution of partial diﬀerential eq uations with complex boundary

conditions and irregular geometry. Interested readers c an refer to more

advanced textb ooks. For a good star t, the book Mathematical Mod-

elling for Earth Sciences by the same a utho r is recommended as it cov-

ers all the main numerical methods and more advanced mathematical

tools in earth sciences.

Bibliography

References

[1] Amundson J. M., Truﬀer M., and Luthi M. P., Time-dependent basal

stress conditions beneath Black Rapids Glacier, Alaska, USA, inferred

from measurements of ice deformation and surface motion, J. of Glaciol-

ogy, 52, 347-357 (2006).

[2] Beard K. V., Terminal velocity and shape of cloud and precipitation

drops aloft, J. Atmosp. Sci., 33, 851-864 (1976).

[3] Bercovic D., Ricard Y., and Richards M. A., The relation between mantle

dynamics and plate tectonics: A primer, in: The History and Dynamics

of Global Plate Motions (Eds Richards M. A., Gordon R. G., van de Hilst

R. O.), AGU Geophysical Monograph 121, AGU Publications ( 2000).

[4] Berger A. L., Long term variations of the Earth’s orbital elements, Ce-

lestial Mechanics, 15, 53-74 (1977).

[5] Brown P., Spalding R. E., ReVelle D. O., Tagliaferri E., Worden S. P.,

The ﬂux of small near-Earth objects colliding with the Earth, Nature,

420, 294-296 (2002).

[6] Committee on Grand Research Questions in the Solid-Earth Sciences,

National Research Council, Origin and Evolution of Earth: Research

Questions for a Changing Planet, National Academies Press, (2008).

[7] Doin M. P. and Fleitout, Thermal evolution of t he oceanic lithosphere:

An alternate view, Earth Planet. Sci. Lett., 142, 121-136 (1996).

[8] Drew D. A., Mathematical modelling of two-phase ﬂow, A. Rev. Fluid

Mech., 15, 261-291 (1983).

[9] England P., Molnar P., and Richter F., John Perry’s neglected critique of

Kelvin’s age for the Earth: A missed opportunity in geodyn amics, GSA

Today, 17, 4-9 (2007).

[10] Fowler A. C., A mathematical model of magma transport in the as-

thenosphere, Geophys. Astrophys. Fluid Dyn., 33, 63-96 (1985).

[11] Fowler A. C. and Yang X. S ., Fast and slow compaction in sedimentary

basins, SIAM J. Appl. Math., 59, 365-385 (1998).

[12] Fowler A. C. and Yang X. S., Pressure solution and viscous compaction

in sedimentary basins, J. Geophys. Res., 104, 12989-12997 (1999).

[13] Fowler A. C. and Yang X. S ., Dissolution-precipitation mecha-

nism for diagenesis in sedimentary basins, J. Geophys. Res., 108,

10.1029/2002JB002269 (2003).

207

208 BIBLIOGRAPHY

[14] Holter O., Ingebretsen F., and Kanestrom I ., Analytic CO

model cal-

culations and global temperature, Eur. J. Phys., 20, 483-494 (1999).

[15] Krige D. G., A statistical approach to some basic mine valuation prob-

lems on the Witwatersrand, J. of the Chem., Metal. and Mining Soc. of

South Africa, 52, 119-139 (1951).

[16] Lemons D. S., Gythiel A., Paul Langevin’s 1908 paper ‘On the Theory

of Brownian Motion’, Am. J. Phys., 65, 1079-1081 (1997).

[17] Mason B. J., Physics of a raindrop, Phys. Educ., 13, 414-419 (1978).

[18] Melosh H. J., Impact Cratering, Oxford University Press, (1989).

[19] Milankovitch M., Canon of isolation and the ice-age problem, K. Serb.

Adad. Geogr., special publication, 132, 484 (1941).

[20] NASA’s Near-Earth Object Science Deﬁnition Team, Study to Deter-

mine the Feasibility of Extending the Search for Near-Earth Objects to

Small er Limiting Diameters, (2003).

[21] Parsons B. and Sclater J. G., Analysis of the variation of ocean ﬂoor

bathymetry and heat ﬂow with age, J. Geophys. Res., 82, 803-827(1977).

[22] Peltier W. R ., Mantle viscosity and ice-age ice sheet topography, Sci-

ence, 273, 1359-1364 (1996).

[23] Price, R., The mechanical paradox of large overthrusts, GSA Bull., 100,

1989-1908 (1988).

[24] Revil A., Pervasive pressure-solution transfer: a poro-visco-plastic

model, Geophys. Res. Lett., 26, 255-258 (1999).

[25] Richter F., Dynamical models for sea ﬂoor spreading, Rev. Geophys.

Space Phys., 11, 223-287 (1973).

[26] Southern California Earthquake Data Centre, www.data.scec.org/

[27] Turcotte D. L. and Schubert G., Geodynamics: Applications of Contin-

uum Physics to Geologi cal Problems, Wiley & S ons, New York, (1982).

[28] USGS Hawaiian Volcano Observatory, hvo.wr.usgs.gov /volcanoes/

[29] Waltham D., Mathematics: a Simple Tool for Geologists, Routeldge

(1997).

[30] Weisstein E. W., http://mathworld.wolfram.com

[31] Wikipedia, http://en.wikipedia.com

[32] Yang X. S. and Young Y., Cellular automata, PDEs and patt ern for-

mation (Chapter 18), in Handbook of Bioinspired Algorithms, edited by

Olarius S. and Zomaya A., Chapman & Hall/CRC, (2005).