Yang X.-S. Introductory Mathematics for Earth Scientists

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

100 Chapter 7. Integration

7.2 Integration by Parts

In diﬀerentiation, we can easily ﬁnd the gradient of x sin x if we know

the derivatives of x and sin x by using the rule of products. In inte-

gration, ca n we ﬁnd the integral of x sin x if we know the integrals of x

and sin x? The answer is yes, this is the integration by parts.

From the diﬀerentiation rule

d(uv)

= u

+ v

, (7.11)

we integrate it with respect to x, we have

uv =

dx +

dx. (7.12)

By rearra nging, we have

dx = uv −

dx. (7.13)

This is the well-known formula for the technique known as integration

by parts. You may wonder where the consta nt of integration is? You

are right, there is a constant of integration for uv, but as we know

they exist for indeﬁnite integrals, we simply omit to write it out in the

formula, but we have to remember to put it back in the end. Now let

us look at a simple example.

Example 7.3: A unit hydrograph is an important tool in hydrology.

There are many diﬀerent methods a n d formulae to approximate the unit

hydrographs such as the triangular representations. From the mathe m a ti-

cal po int of view, the following form can ﬁt well

f(t) = At

−t/τ

where A is a sc a ling constant, and τ is the time c onstant. The exponent

n > 0 and τ as well as A can be obtained by ﬁtting to the experimental

data.

The total discharge from 0 to t = T can be calcul a te d by the integral

Q =

f(t)dt.

Let us consider a s pecial case of n = 1 for simplicity. Now we have

Q = A

−t/τ

dt.

7.2 Integration by Parts 101

Setting u = t and v

= e

−t/τ

and using integration by part, we have

Q = A{

t(−τ)e

−t/τ

−

−t/τ

(−τ)dt}

= A{−τT e

−T/τ

+ τ

(−τ)e

−t/τ

} = Aτ

− Aτ(τ + T )e

−T/τ

We can see that

Q → Aτ

, as T → ∞.

So for a given discharge Q, the characteristic discharge time is

τ =

When applying the technique of integration by parts, it might be

helpful to paus e and look at diﬀerent perspectives, and the solutions

sometimes come out more naturally. For example, we try to ﬁnd the

integral

I =

sin xdx. (7.14)

If we us e u = e

and dv/dx = sin x, we have du/dx = e

and v =

−cos x. By using integratio n by parts, we then have

I =

sin xdx = e

(−cos x) −

(−cos x)

= −cos xe

cos xdx. (7.15)

It seems that we are stuck here. The integration by parts does not seem

to help much, as it only transfers

sin xdx to

cos xdx, which is

not an improvement. In this case, we have to pause and think, what

happens if we use the integration by parts a gain for

cos xdx? We

let u = e

and dv/dx = cos x, which gives du/dx = e

and v = sin x,

and we have

cos xdx = e

sin x−

sin x

dx = e

sin x−

sin xdx. (7.16)

Then, you may think, what is the use, as we have come back to the

original integral? But what happens if we combine the above two equa-

tions? We have

I =

sin xdx = −e

cos x −

cos xdx

= −e

cos x + e

sin x −

sin x = e

(sin x − cos x) − I. (7.17)

102 Chapter 7. Integration

We want to ﬁnd I, so we can solve the above equation by simple rear-

rangement

2I = e

(sin x − cos x), (7.18)

I =

(sin x − cos x). (7.19)

Here we have found the integral as

I =

sin xdx =

(sin x − cos x) + C. (7.20)

Such a proc e dure is very useful in deriving r e ductio n formulae for some

integrals . Let us take the Wallis’ formula as an ex ample.

Example 7.4: In order to derive Wallis’ formula for the integral

π/2

sin

xdx, (n ≥ 2). (7.21)

We use u = sin

n−1

x and dv/dx = sin x; we get

= (n − 1) sin

n−2

x cos x, v =

sin x = −cos x.

We have

− sin

n−1

x cos x

π/2

−

π/2

(−cos x)(n − 1) sin

n−2

cos xdx

= 0 + (n − 1)

π/2

cos

x sin

n−2

dx.

Since sin

x + cos

x = 1 or cos

x = 1 − sin

x, we have

= (n − 1)

π/2

(1 − sin

x) sin

n−2

= (n−1)

π/2

sin

n−2

dx−(n−1)

π/2

sin

xdx = (n−1)I

n−2

−(n−1)I

whose rearrangement leads to

n − 1

n−2

which is Wallis’ reductio n formula for the i n te gral I

. If we continue to

use the reduction formula, when n i s odd, we have

(n − 1)(n − 3)...4 × 2

n(n − 2)...5 × 3

7.3 Integration by Substitution 103

When n is even, we have

(n − 1)(n − 3)...3 × 1

n(n − 2)...4 × 2

It is left as an exercise to prove that the integral

π/2

cos

xdx

follows the same formulae.

When using integration by parts, sometimes we have to try to use

diﬀerent combinations for u and v. In the following example, we have

to use u = x and v

= xe

−x

. If we try to use u = x

, it will not work.

Example 7.5: Let us now try to calcula te the integral

J =

∞

−∞

−x

dx,

using integration by parts or

udv = uv −

vdu. By setting u = x and

= xe

−x

, we have

= 1, v =

−x

dx =

−x

d(x

) = −

−x

We have

J =

∞

−∞

−x

dx =

−

−x

∞

−∞

−

∞

−∞

(

−1

−x

)dx

= 0 +

∞

−∞

−x

dx =

√

erf(x)

∞

−∞

√

× [1 − (−1)] =

√

where we have used the error function

erf(x) =

√

−u

du,

and erf(∞) = 1 and erf(−∞) = −1. We wil l use this result in deriving

the variance of the Gaussian distribution.

7.3 Integration by Substitution

Sometimes, it is not easy to carry out the integration directly. However,

it might become easier if we use the change of variables or integration

by substitution. For exa mple, we want to calculate the integral

I =

f(x)dx. (7.22)

104 Chapter 7. Integration

We c an change the variable x into another variable u = g(x) where

g(x) is a known function of x. This mea ns

= g

(x), (7.23)

du = g

(x)dx, dx =

(x)

du. (7.24)

This means that

I =

f(x)dx =

f[g

−1

(u)]

du, (7.25)

where it is usually not necessary to calculate g

−1

(u) as it is relatively

obvious. For example, in order to do the integration

I =

dx, (7.26)

we let u = x

+ 5 and we have

= (x

+ 5)

= 2x, (7.2 7)

du = 2xdx, (7.28)

which means dx =

du. Therefore, we have

I =

dx =

du =

+ A = e

+ A, (7.29)

where A is the constant of integration. Here we have substituted u =

+ 5 back in the last step.

Example 7.6: Compaction occurs when loos e ly packaged materials

(such as soils and snow) deform under their own weight, resulting in the

reduction of porosity. Compacted sediments will form sedimentary rocks,

and compacted snow w ill form glaci e rs. In the simplest case, the porosity

φ will decrease a lmost expon e n tially in terms of

φ = φ

−αz

where z is the depth of the column of the porous materials, and φ

is the

initial porosity at z = 0 . The parameter α can be considered as a constant,

typically α = 0.001/m. In sedimentary basins, the pores can be considered

fully saturated with ﬂuid, and the reduction of pore space is often linked

7.4 Bouguer Gravity 105

gure 7.4: Derivation of the Bouguer correction formula.

with the generation of hydrocarbon under the right conditions. Suppose a

fraction (say, β = 0.01 or about 1%) of such ﬂuids has been turned into

oil, then the total volume c a n be estimated by the total reduction of the

pore space

V = Aβ

(φ

− φ)dz = βA

(1 − e

−αz

)dz,

where A is the area of interest. Setting u = αz, we have

V = βAφ

dz − βAφ

αh

−u

du = βAφ

h −

βA

−αh

− 1].

When changing variables, the integration limits should also change

accordingly. Otherwise, the results will be diﬀerent and incorrect. This

will be demonstrated by the application in the next section.

7.4 Bouguer Gravity

The Boug uer formula for gravity correction is derived using an inﬁnite

plate of thickness t with a uniform density ρ. As the plate is inﬁnite,

we can choose any point as the origin r = 0 and cut a small ring, and

we ar e only interested in the gravity at point P outside the plate at a

perpendicula r distance of h away (see Fig. 7.4).

Due to symmetry, the gravitational force will point towards p oint

O. According to Newton’s law of gravitation, the gravity caused by the

small ring is

dg = −

Gdm

cos θ, (7.30)

where G is the universal gr avitational c onstant, θ is the angle. R is the

distance from P to the ring. The ma ss of the ring is dm = ρtdA and

106 Chapter 7. Integration

dA = 2πrdr is the area of the ring with the thickness of dr. Thus, we

have

dm = 2πρtrdr. (7.31)

In addition, from the basic trigonometry, we know that

cos θ =

, R =

+ r

. (7.32)

Therefore, the gravity due to the ring b e c omes

dg = −2πGρht

rdr

= −2πGρht

rdr

+ r

)

3/2

. (7.33)

So the total gravity due to the whole plate is the integration

∆g = −2πGρht

∞

rdr

+ r

)

3/2

. (7.34)

Using the change of variables by setting u = (r

) so that du = 2rdr

and u = h

for r = 0, we have

∆g = −2πGρht

∞

3/2

. (7.35)

Since

−3/2

du = −2u

−1/2

, we have

∆g = −2πGρht

∞

−

1/2

∞

= −2πGρht[−(0 −

)] = −2πGρt. (7.36)

This is exactly the well-known formula for Bouguer reduction of a layer

with a density of ρ and thickness t.

Using the typical values of G = 6.67 ×10

−11

kg, a nd dens ity

ρ = 2670 kg/m

for rocks, we have

∆g = −2π × 6.67 × 10

−11

× 2670 t

≈ 1.1189 × 10

−6

t (m/s

) ≈ 0.112 t (mGal). (7.37)

This is to say, for e very 1-metre layer removed, the Bouguer gravity

anomaly is corrected by −0.112 mGal. The Bouguer gravity variations

reﬂect changes in mass distr ibution below the surface after removing

most of the e ﬀect of mass excess above the reference level, o ften the

sea level.

Chapter 8

Fourier Transforms

Mathematical transform is a method of changing one kind of function

and equation into another, often simpler and easier to solve. Fourier

transform maps a function in the time domain such as a sig nal into

another function in the frequency domain, which is commonly used in

signal process ing. In earth sciences, the processing of seismic signals is

pivotal to oil and gas exploratio n and the understanding of almost all

underground geological structures.

As an application, we will analyse the Milankovitch cycles using

simple Fourier transform in the last s e c tion o f this chapter.

8.1 Fourier Series

From earlier discussions, we know that the function e

can be ex panded

into a series in terms of a polynomial with terms of x, x

, ..., x

In this case, we ar e in fact trying to expand it in terms of the basis

functions 1, x, x

, ..., and x

. There are many other basis functions.

For example, the basis functions sin(nπt) and cos(nπt) are more widely

used in seismic signal processing. In general, this is essentially about

the Fourier series. For a function f(t) on an interval t ∈ [−T, T ] where

T > 0 is a ﬁnite constant or half period, the Fourier series is deﬁned

f(t) =

∞

n=1

cos(

nπt

) + b

sin(

nπt

)

, (8.1)

where

−T

f(t)dt, a

−T

f(t) cos(

nπt

)dt, (8.2)

and

−T

f(t) sin(

nπt

)dt, (n = 1, 2, ...). (8.3)

107

108 Chapter 8. Fourier Transforms

Here a

and b

are the Fourier coeﬃcients of f (t) on [−T, T ]. The

function f(t) can be continuous or piecewise continuous with a ﬁnite

number of jump discontinuities. For a jump disc ontinuity at t = t

, if

−) and f

+) both exist with f(t

−) 6= f(t

+), then the Fourier

series converge to the average value. Tha t is

f(t

) =

[f(t

−) + f(t

+)]. (8.4)

You may wonder how to calculate the coeﬃcient a

and b

? Before

we proceed, let us prove the orthogonality relation

J =

−T

sin(

nπt

) sin(

mπt

)dt =







0 (n 6= m)

T (n = m)

, (8.5)

where n and m are integers.

From the trigonometrical functions, we know that

cos(A + B) = cos A cos B − sin A sin B, (8.6)

and

cos(A − B) = cos A cos B + sin A sin B. (8.7)

By subtracting, we have

cos(A − B) −cos(A + B) = 2 sin A sin B. (8.8)

Now the orthogonality integral becomes

J =

−T

sin(

nπt

) sin(

mπt

)dt

−T

cos[

(n − m)πt

] − cos[

(n + m)πt

]

dt. (8.9)

If n 6= m, we have

J =

(n − m)π

sin[

(n − m)πt

]



−T

−

(n + m)π

sin[

(n + m)πt

]



−T

[

(n − m)π

× (0 − 0) −

(n + m)π

× (0 − 0)] = 0. (8.10)

If n = m, we have

J =

−T

1 − cos[

2nπt

]

8.1 Fourier Series 109



−T

−

2nπ

sin[

2nπt

]



−T

[2T −

2nπ

× 0] = T, (8.11)

which proves the relation (8.5). Using similar calculations, we can easily

prove the following ortho gonality relations

−T

cos(

nπt

) cos(

mπt

)dt =







0 (n 6= m)

T (n = m)

, (8.12)

and

−T

sin(

nπt

) cos(

mπt

)dt = 0, for all n and m. (8.13)

Now we can try to derive the expression for the coeﬃcients a

Multiplying b oth sides of the Fourier series (8.1) by cos(mπt/T ) and

taking the integration from −T to T , we have

−T

f(t) cos(

mπt

)dt =

−T

cos(

mπt

)dt

∞

n=1

−T

cos(

nπt

) cos(

mπt

)dt + b

−T

sin(

nπt

) cos(

mπt

)dt

Using the relations (8.12) and (8.13) as well as

−T

cos(mπt/T )dt = 0,

we know that the only non-zero integral on the right-hand side is when

n = m. Therefore, we get

−T

f(t) cos(

nπt

)dt = 0 + [a

T + b

× 0], (8.14)

which gives

−T

f(t) cos(

nπt

)dt, (8.15)

where n = 1, 2, 3, ···. Interestingly, when n = 0, it is still valid and

becomes a

as cos 0 = 1. That is

−T

f(t)dt. (8.16)

In fact, a

/2 is the average of f (t) over the period 2T . The coeﬃcients

can be obtained by multiplying sin(mπt/T ) and following similar

calculations.

Fourier series in general tends to converge slowly. In order for a

function f(x) to be expanded, it must satisfy the Dirichlet conditions:

f(x) must be periodic with a t most a ﬁnite number of discontinuities,