Partial Differential Equations

Partial Differential Equations Herman Jaramillo May 10, 2016

2

Contents 1 First Order 1.1 Finding the characteristics . . . . . . . . . . . 1.1.1 Initial conditions . . . . . . . . . . . . 1.1.2 Description of the General Algorithm: book . . . . . . . . . . . . . . . . . . .

. . . . . . . . . . . . Following . . . . . .

. . . . . . . . . . . . Bleistein’s . . . . . .

. . . . [2] . .

7 7 8 8

2 second order 2.1 constant coefficients . . . . . . . . . . . . . . . . . . . . . . . . . . . . . 2.1.1 varying coefficients: general hyperbolic equation . . . . . . . . . 2.2 Problems . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . 2.3 Semi–infinite string . . . . . . . . . . . . . . . . . . . . . . . . . . . . . 2.3.1 Mixed problem for the d’Alembert equation . . . . . . . . . . . 2.3.2 Other boundary conditions . . . . . . . . . . . . . . . . . . . . . 2.3.3 Discontinuities of a solution along a principal characteristic curve. Continuity conditions . . . . . . . . . . . . . . . . . . . . . . . .

11 11 15 18 31 31 32

3 Green Functions 3.1 Problems . . . . . . . . . . . . . . . . 3.1.1 The sifting property of a delta 3.1.2 Problem 23.7 c . . . . . . . . 3.1.3 Problem 23.7 d . . . . . . . .

39 39 39 40 43

. . . . . . . distribution . . . . . . . . . . . . . .

4 Higher dimensions 4.1 Problems . . . . . . . . . . . . . . . . . . . . . . 4.1.1 Five different derivations of the volume hyper–sphere . . . . . . . . . . . . . . . 4.1.2 Fundamental solutions of the Laplacian . 4.1.3 The Helmholtz equation . . . . . . . . . 4.1.4 The Cauchy–Riemann equation . . . . . Appendices

. . . .

. . . .

. . . .

. . . .

. . . .

. . . .

. . . . . . . and surface . . . . . . . . . . . . . . . . . . . . . . . . . . . .

. . . .

. . . .

. . . .

. . . area . . . . . . . . . . . .

. . . .

. . . .

. . . .

. . . of a . . . . . . . . . . . .

33

47 47 47 65 70 71 72

4 A The generalized polar–spherical coordinates A.1 Derivation . . . . . . . . . . . . . . . . . . . A.2 The Jacobian . . . . . . . . . . . . . . . . . A.2.1 Partial Derivatives . . . . . . . . . . A.3 The Volume and Surface Area . . . . . . . .

CONTENTS

. . . .

B Derivations of wave equations B.1 One dimensional compressional wave equation B.1.1 From a system of springs to a bar . . . B.1.2 From continuum mechanics tools . . . B.2 One dimensional transverse wave equation . . B.3 Elastic Waves in 3D Anisotropic Media . . . . B.3.1 Cauchy’s Lemma . . . . . . . . . . . . B.3.2 Cauchy’s Law . . . . . . . . . . . . . .

. . . .

. . . . . . .

. . . .

. . . . . . .

. . . .

. . . . . . .

. . . .

. . . . . . .

. . . .

. . . . . . .

. . . .

. . . . . . .

. . . .

. . . . . . .

. . . .

. . . . . . .

. . . .

. . . . . . .

. . . .

. . . . . . .

. . . .

. . . . . . .

. . . .

. . . . . . .

. . . .

. . . . . . .

. . . .

73 73 78 78 82

. . . . . . .

85 86 86 87 88 90 91 93

Introduction In chapter 1 I show the method of characteristics by a simple example for first order, quasi–linear Partial Differential Equations (PDEs). The second order PDEs are also studied using the method of characteristics in chapter 2. Chapter 3 deals with Fourier methods, and chapter 4 with Green functions.

6

CONTENTS

Chapter 1 First Order Solve ∂U ∂U + = 0. ∂t ∂x

1.1

(1.0.1)

Finding the characteristics

Let us define a vector of the coefficients of the equation as:   1 a= . 1 Now it is clear that equation 2.1.9 is: a · ∇U = 0. So the vector a is perpendicular to the gradient ∇U , and therefore parallel to contour levels of U defined by U = (x, t) = constant. If C is a contour level of the field U along its arc length s (assuming it is rectifiable) then ∂U dx ∂U dt dU = + = 0. (1.1.2) ds ∂x ds ∂t ds So we see that the vector b defined as:   dx  ds    b= .  dt  ds is also perpendicular to the gradient ∇U and so it is parallel to a. Therefore we have: dx = 1.F (s) ds (1.1.3) dt = 1.F (s) ds

8

First Order

Given that  t·t=

dx ds

2

 +

dt ds

2

then dx/ds and dt/ds cannot vanish simultaneously. We use equation 1.1.3 to find the curve C. From 1.1.3 dx dx/ds = = 1/1 = 1. dt dt/ds This is the equation of the characteristics. That is (integrating), the path C satisfies x=t+c where c is a constant. So we can say that U is constant along all the 45 degree lines that cross the t axis. The next step is to solve the equation along the characteristic.

1.1.1

Initial conditions

Because the PDF 2.1.9, is equivalent to an ordinary differential equation along the characteristic curve, Eq. is is not possible to prescribe the dependence of U along a characteristic curve. Let N be any curve such that the slope of N at any point (x, t) is not equal to the slope of the characteristic curve through (x, t). Prescribing the initial values of U on N is enough to determine a unique solution of the original PDF 2.1.9 on each characteristic curve that issues from a point on N . Therefore the values of U on N uniquely determine U , and conversely. Since U is constant along the characteristic we have that the value of U is only determined by the value c in the equation x = t + c, that is by the value x − t. That is, U (x, t) = f (x − t) = f (c) Let as assume that the initial condition is f (x) = sin 4πx at the time t = 0. So the solution is U (x, t) = f (x − t) = sin 4π(x − t).

1.1.2

Description of the General Algorithm: Following Bleistein’s [2] book

Given the firs–order quasi–linear equation a(x, y, u)ux + b(x, y, u)uy = c(x, y, u).

1.1 Finding the characteristics

9

Find the equations of the characteristics: b c a = = dx dy du Note that I inverted the fractions since I want c upstairs, to consider the homogeneous case c = 0. Bleistein’s idea here is to parametrize the solution surface U (x, y) by two parameters τ and σ such that x = x(σ, τ ) ,

y = y(σ, τ ) ,

U = U (σ, τ ).

(1.1.4)

In such a way that one parameter is along the characteristics (say σ) and the other (τ ) in a trajectory that cuts through the characteristics (nowhere kissing them). The equations of the characteristics in the new parameterization are: dx = λa , dσ

dy = λb , dσ

du = λ c. dσ

(1.1.5)

We can pick λ= √

a2

1 + b2 + c 2

to make σ be the arc length along the characteristics, or λ= √

1 a2 + b 2

to make σ be the arc length along the base characteristics, or simply λ = 1. Data specification For example if σ = 0, we can have a curve parameterized by τ x = x(0, τ ) ,

y = y(0, τ ) ,

U = U (0, τ ).

For each fixed τ we solve equations 1.1.5 and get to the solution in the form of 1.1.4. In order, for the parameterization 1.1.4 to work we need a non–vanishing Jacobian. That is: xσ xτ 6= 0. J = yσ yτ Bleistein show that the Jacobian vanishes along tangential (kissing) points of the initial curve with the characteristics. If the Jacobian vanishes everywhere then the initial curve is a characteristic. There could be infinite solutions.

10

First Order

The method of envelopes We describe the method following Bleistein’s example ux + uy = u The characteristics are dx dy du = = 1 1 u The integration of the characteristics in terms of “x” is: y = x + c1

,

u = c2 ex

(1.1.6)

Each choice of the constants c1 , c2 produces a characteristic curve. That is, 1.1.6 is a two–parameter family of characteristic curves. A solution surface could be determined by prescribing a functional relationship between c1 and c2 , that i,s by prescribing, for example, c2 = f (c1 ) Or equivalently u e−x = f (y − x) solving for u we conclude u = ex f (y − x) By direct substitution one can verify that this is a solution to the differential equation

Chapter 2 second order 2.1

constant coefficients

The characteristic method for solving the second order partial differential equations is based on the book by [3]. Let us consider the equation of the form Au ≡ a

∂ 2u ∂ 2u ∂ 2u + 2 b + c = 0, ∂t2 ∂t ∂x ∂x2

x∈R t>0

(2.1.1)

with initial conditions u(x, 0) = ϕ(x) ∂u(x, t) = ψ(x). ∂t t=0

(2.1.2) We further assume that b2 − ac > 0 (so that equation 2.1.1 is hyperbolic). The idea behind the method of characteristics is to change this second order equation into two simpler equation in terms of first order derivatives. We do this by factoring the equation. To simplify notation let us define T =

∂ ∂t

X=

∂ ∂x

Now using this notation into equation 2.1.1 we find aT 2 + 2bT X + cX 2 = 0.

(2.1.3)

This changes equation 2.1.1 into an algebraic equation which is easy to factor. The roots of equation 2.1.3 are given by √ √ −2 bX ± 4b2 X 2 − 4a cX 2 −b ± b2 − a c = X = Xλ1,2 2a a

12

second order

with λ1,2 =

−b ±

√

b2 − a c X. a

since b2 − 4a c > 0 the two roots are real. We then write (T − λ1 X)(T − λ2 X) = 0. We write the factored differential equation as    ∂ ∂ ∂ ∂ Au ≡ − λ1 − λ1 u = 0, ∂t ∂x ∂t ∂x

(2.1.4)

We split this equation into two     ∂ ∂ ∂ ∂ − λ1 , u= u · (−λ1 , 1) = ∇u · (−λ1 , 1) = 0 L(−λ1 ,1) u = ∂t ∂x ∂x ∂t     ∂ ∂ ∂ ∂ − λ2 , u= u · (−λ2 , 1) = ∇u · (−λ2 , 1) = 0 L(−λ2 ,1) u = ∂t ∂x ∂x ∂t (2.1.5) and by using the definition of directional derivative, we conclude that the derivative along the directions (−λ1 , 1) and (−λ2 , 1) is zero. Any line in the (x, t) plane, with direction (−λ1 , 1) has the slope −1/λ1 and so can be written as t − t0 = −

1 (x − x0 ). λ1

which we can write as x + λ1 t = λ1 t0 + x0 = ξ along the same lines we find that the directional derivative of u is zero along any line of the form x + λ2 t = λ2 t0 + x0 = η That is the directional derivatives of u along the lines ξ = const and η = const are zero. This lines are know as the “characteristics” for the differential equation 2.1.1. This provides a method to change variables in a way to simplify the second order differential equation 2.1.1. We define the coordinate transformation u = u(x(ξ, η), t(ξ, η)). Then ∂ξ ∂ξ − λ1 = λ1 − λ1 = 0 ∂t ∂x ∂η ∂η L(−λ2 ,1) η = − λ1 = λ2 − λ2 = 0 ∂t ∂x L(−λ1 ,1) ξ =

2.1 constant coefficients

13

L(−λ1 ,1) is the operator of differentiation along the lines ξ = const, while L(−λ2 ,1) is the operator of differentiation along the lines η = const. So, ∂ ∂ L(−λ1 ,1) = c1 ; L(−λ2 ,1) = c2 (2.1.6) ∂η ξ=const ∂ξ η=const From equations 2.1.4 and 2.1.6 we find ∂ 2u = 0. ∂η ∂ξ

(2.1.7)

This is the canonical form of the second order hyperbolic PDE. To solve this equation, let denote ∂u (ξ, η) = v(ξ, η). ∂η

(2.1.8)

so equation 2.1.7 can be written as ∂v dv = = 0. ∂ξ dξ η=const It then follows that v|η=const does not depend on ξ, that is, v(ξ, η) = c(η), and from 2.1.8 d u = c(η). dη ξ=const Integrating this ordinary differential equation, we obtain Z u = c(η)dη + c1 (ξ). ξ=const

Thus, u(x, t) = f (ξ) + g(η) = f (x + λ1 t) + g(x + λ2 t). where f and g are functions of one variable. To find the form of f and g we use the initial conditions 2.1.2. From the first of these conditions we see that u(x, 0) = f (x) + g(x) = ϕ(x), and from the second condition ∂u(x, t) = λ1 f 0 (x) + λ2 g 0 (x) = ψ(x); ∂t t=0

(2.1.9)

14

second order

we integrate this equation between 0 and x and find Z x ψ(s) ds λ1 f (x) + λ2 g(x) =

(2.1.10)

0

by inserting 2.1.9 into 2.1.10 Z λ1 [ϕ(x) − g(x)] + λ2 g(x) =

x

ψ(s) ds

(2.1.11)

0

so g(x) = f (x) = = =

Z x 1 λ1 ψ(s) ds − ϕ(x) λ2 − λ1 0 λ2 − λ1 ϕ(x) − g(x)   Z x 1 λ1 ϕ(x) − ψ(s) ds 1+ λ2 − λ1 λ2 − λ1 0 Z x λ2 1 ϕ(x) − ψ(s) ds λ2 − λ1 λ2 − λ1 0

so u(x, t) = f (x + λ1 t) + g(x + λ2 t) λ2 λ1 ϕ(x + λ1 t) − ϕ(x + λ1 t) + = λ2 − λ1 λ2 − λ1 Z x+λ2 t Z x+λ1 t 1 1 ψ(s) ds − ψ(s) ds. λ2 − λ1 0 λ2 − λ1 0 Finally we obtain the d’Alembert solution λ2 ϕ(x + λ1 t) − λ1 ϕ(x + λ2 ) 1 u(x, t) = + λ2 − λ1 λ2 − λ1

Z

x+λ2 t

ψ(s) ds. x+λ1

example Let us consider the simple equation 2 ∂ 2 u(x, t) 2 ∂ u(x, t) = a = 0, ∂t2 ∂x2

−∞ < x < ∞,

t > 0.

with 2.1.2 initial conditions. The corresponding characteristic polynomial is given by λ2 − a2 = (λ + a)(λ − a). with roots λ1,2 = ∓a The d’Alembert solution is then given by Z x+at ϕ(x + at) + ϕ(x − at) 1 u(x, t) = + ψ(s) ds. 2 2a x−at

(2.1.12)

2.1 constant coefficients

2.1.1

15

varying coefficients: general hyperbolic equation

We assume the equation ∂ 2u ∂ 2u ∂ 2u ∂u + c(x, t) + 2 b(x, t) + d(x, y) 2 2 ∂t ∂t ∂x ∂x ∂t ∂u + e(x, t) + f (x, y)u = g, x ∈ R t>0 ∂x

Au ≡ a(x, t)

(2.1.13)

where now the coefficients are varying. For the moment we will ignore the initial conditions and focus in the general solution. We want to find a solution around a point (x0 , y0 ) that 2 ∆ = b (x, t) − a(x, t) c(x, t) > 0. x=x0 ,y=y0

We call ∆ the discriminant and the point (x0 , y0 ) a hyperbolic point. At this moment I consider convenient to change notation. New notation From now on, partial derivatives are recognized by a subindex. For example

···

∂u ∂t ∂ 2u ∂t2 ∂u ∂x ∂ξ ∂x and so on

= ut = utt = ux = ξx

observe that each time we use this notation we change 5 characters by 2 characters on a first derivative and from 7 to 3 in a second derivative, saving lots in typing without loosing the clarity in exposition. In the new notation equation 2.1.14 becomes a(x, t) utt + 2 b(x, t) utx + c(x, t) uxx + d(x, t) ut + e(x, t)ux + f (x, t)u = g, x∈R t>0 (2.1.14) Change of coordinates A way to solve equations is by a change of variables that simplify the equation. In the case of equation 2.1.14 we want to defined new coordinates ξ, η such that the equation

16

second order

in the new coordinate system is simpler (or canonical). That is we would like to have equation 2.1.14 converted to the form b∗ (ξ, η)uξ η (ξ, η) + d∗ (ξ, η)uξ + e∗ (ξ, η)uη + f ∗ (ξ, η)u = g ∗ (ξ, η). with some new coefficients b∗ (ξ, η), d∗ (ξ, η), e∗ (ξ, η), f ∗ (ξ, η) and, g ∗ (ξ, η). Note that the idea is to annihilate the coefficients a∗ (ξ, η) and c∗ (ξ, η). In what follows we write a∗ instead of a∗ (ξ, η) to simplify notation. Let us then define the coordinate transformation ξ = ξ(x, t) η = η(x, t) (2.1.15) This coordinate transformation is well defined if ξt ξx = ξt ηx − ηt ξx 6= 0. J = ηt ηx

(2.1.16)

We first find the first partial derivatives of u with respect to t and x and then the second partial derivatives. By using the chain rule ut = uξ ξt + uη ηt ux = uξ ξx + uη ηx Applying the chain rule again utx = (ut )x = (uξ ξt + uη ηt )x = (uξ ξt )x + (uη ηt )x = (uξ )x ξt + uξ ξtx + (uη )x ηt + uη ηt x we use the second of equations 2.1.17 and so utx = (uξξ ξx + uηξ ηx ) ξt + uξ ξt x + (uηξ ξx + uηη ηx )ηt + uη ηtx and rearranging, utx = uξξ ξt ξx + uξη (ξt ηx + ξx ηt ) + uηη ηt ηx + uξ ξtx + uη ηtx .

(2.1.17)

if the derivative is with respect to tt, change x by t in the previous equation, and if the derivative is with respect to xx, change t by x in the previous equation, so utt = uξξ ξt2 + 2uξη ξt ηt + uηη ηt2 + uξ ξtt + uη ηtt .

(2.1.18)

uxx = uξξ ξx2 + 2uξη ξx ηx + uηη ηx2 + uξ ξxx + uη ηxx .

(2.1.19)

and

2.1 constant coefficients

17

we plug equations 2.1.17, 2.1.18, and 2.1.19 into the original equation 2.1.14 and find a∗ uξξ + 2 b∗ uξη + c∗ uηη + d∗ uξ + e∗ uη + f ∗ u = g ∗ ,

(2.1.20)

where the new coefficients a∗ , b∗ and c∗ are the following functions of ξ and η a∗ b∗ c∗ d∗ e∗ f∗ g∗

= = = = = = =

a ξt2 + 2b ξt ξx + c ξx2 a ξt ηt + b (ξt ηx + ξx ηt ) + c ξx ηx a ηt2 + 2b ηt ηx + cηx2 a ξtt + b ξtx + cξtt + dξt + eξx a ηtt + b ηtx + cηxx + dηt + eηx f g. (2.1.21)

We now show that the discriminant ∆∗ in the new coordinate system is still greater than zero. That is, the point (x0 , y0 ) in the new coordinate system (ξ, η) is still hyperbolic. Let us see, after sorting in the order b2 , ac, a2 , c2 , ab, bc ∆∗ = (b∗ )2 − (a∗ )(c∗ ) = b2 [(ξt ηx + ξx ηt )2 − 4 ξt ξx ηt ηx ] − ac[ξt2 ηx2 + ηt2 ξx2 − 2ξt ηt ξx ηx ] + a2 [ξt2 ηt2 − ξt2 ηt2 ] + c2 [ξx2 ηx2 − ξx2 ηx2 ] + ab[2 ξt ηt (ξt ηx + ξx ηt ) − 2 ξt2 ηt ηx − 2 ηt2 ξt ξx ] + bc[2 ξx ηx (ξt ηx + ξx ηx ) − 2 ηx2 ξt ξx − 2 ξx2 ηt ηx = b2 (ξt ηx − ξx ηt )2 − ac(ξt ηx − ξx ηt )2 + a2 0 + a b 0 + b c 0 = (b2 − ac)J 2 where J is given in equation 2.1.16, and since J 6= 0 then ∆∗ > 0 and the change of coordinates preserve the hyperbolic point. The idea is to make a∗ = c∗ = 0. Since J 6= 0, then b 6= 0. That is, we want a ζt2 + 2b ζt ζx + c ζx2 = 0. (2.1.22) where ζ is either ξ or η. We divide by ζx2 and find    2 ζt ζt + 2b +c = 0 a ζx ζx (2.1.23) If we pick ζ(x, t) = const then dζ = ζt dt + ζx dx = 0,

18

second order

which is ζt dx =− , ζx dt and equation 2.1.23 becomes  a

dx dt

2

 − 2b

dx dt

 + c = 0. (2.1.24)

We arrived to a quadratic equation of ordinary derivatives. A mnemotechnic trick to get to here is to make the substitution ∂ 7→ dx ∂t

∂ 7→ −dt ∂x

(2.1.25)

on the original equation 2.1.13 and then divide by dt2 . Equation 2.1.25 turns into the two ordinary differential equations √ b ± b2 − ac dx = λ1,2 = dt a The solutions, which are the characteristics, are integrated and the two constants of integration named ξ and η. We proceed to illustrate the methods by solving the problems in Problem 4.6. section of Komech’s [3]’s book.

2.2

Problems

1. . Find the general solution for the following equations: (a) uxx + 2uxy − 3uyy + 2ux + 6uy = 0. solution: To find the coordinates that transform this form into a canonical form, by using the mapping 2.1.25, on the quadratic terms, then the characteristics are given by the solutions of 

dy dx

2

 −2

dy dx

 − 3 = 0.

Let us identify the following coefficients a = 1,

b = 1,

c = −3,

d = 2,

e = 6,

f = 0,

g = 0.

2.2 Problems

19

This equation can be factored out as    dy dy −3 + 1 = 0, dx dx from which we find two solutions Z Z dy = 3 dx y − 3x = ξ and Z

Z dy = −

dx

(2.2.26)

y+x=η From equations 2.2.26 we find ξx = −3, ξy = 1, ξxx = 0, ξxy = 0, ξyy = 0 ηx = 1, ηy = 1, ηxx = 0, ηxy = 0, ηyy = 0 (2.2.27) With these values we build the coefficients of the new differential equation. As checking we first verify that a∗ = c∗ = 0. a∗ = a ξx2 + 2b ξx ξy + c ξy2 = 9 − 6 − 3 = 0 c∗ = a ηx2 + 2b ηx ηy + cηy2 = 1 + 2 − 3 = 0 so indeed a∗ = b∗ = 0. Let us now compute the rest of coefficients.

b∗ d∗ e∗ f∗ g∗

= = = = =

a ξx ηx + b (ξx ηy + ξy ηx ) + c ξy ηy = −3 + (−3 + 1) − 3 = −8 a ξxx + b ξxy + cξyy + dξx + eξy = 0 + 0 + 0 − 6 + 6 = 0 a ηxx + b ηxy + cηyy + dηx + eηy = 0 + 0 + 0 + 2 + 6 = 8 f =0 g = 0,

with this, equation 2.1.20 turns out to be −16uξη + 8uη = 0. That is 2 uξη − uη = 0.

(2.2.28)

20

second order This equation is in canonical form. Let v = uη . So we have the equation 2vξ − v = 0. This is an ordinary differential equation in ξ that we can integrate. That is Z Z dv dv 2 =v⇒2 + C(η) = dξ, dξ v where C(η) does not depend on ξ. Then 2 ln v + C(η) = ξ and this can be written as v(ξ, η) = ϕ(η) eξ/2 . and because v = uη then du = φ(η) eξ/2 ⇒ du = φ(η) eξ/2 dη dη we integrate this equation with respect to η and find Z Z ξ/2 ξ/2 u = ϕ(η)e dη + ψ(ξ) = e ϕ(η)dη + ψ(ξ). we can rename Z F (η) =

ϕ(η)dη

G(ξ) = ψ(ξ) and rewrite u = eξ/2 F (η) + G(ξ). in terms of x and y this is u = e−(3x−y)/2 F (x + y) + G(y − 3x). (b) 1 x uxx − yuyy + (ux − uy ) = 0, 2

x > 0 y > 0.

2.2 Problems

21

solution: Let us first, normalize the equation. Since x > 0 we can divide by x so we can write the equation as 1 y (ux − uy ) = 0, uxx − uyy + x 2x

x > 0 y > 0.

Here the coefficients are: a = 1,

b = 0,

y c=− , x

d=

1 , 2x

e=−

1 , 2x

f = g = 0. (2.2.29)

The characteristics are given by the solution of  2 dy y − = 0, dx x which can be factored out as r  r   dy y dy y − + = 0. dx x dx x The solutions are given by Z

Z dx dy √ +ξ √ = y x Z Z dx dy −√ + η √ = y x

That is √ √ ξ = 2( y − x) √ √ η = 2( y + x),

(2.2.30)

and from here 1 1 1 √ , ξx = − √ , ξy = √ , ξxx = y x 2x x 1 1 1 ηx = √ , ηy = √ , ηxx = − √ , y x 2x x

ηxy = 0 ηyy

We now verify that a∗ = c∗ = 0. 1 y1 − =0 x xy 1 y1 = a ηx2 + 2b ηx ηy + cηy2 = − = 0. x xy

a∗ = a ξx2 + 2b ξx ξy + c ξy2 = c∗

1 √ 2y y 1 =− √ . 2y y

ξxy = 0 ξyy = −

22

second order The other coefficients are: b∗ = a ξx ηx + b (ξx ηy + ξy ηx ) + c ξy ηy 1 y 1 = −√ − √ xy x xy 1  y = −√ 1+ xy x ∗ d = a ξxx + b ξxy + cξyy + dξx + eξy 1 y 1 1 1 1 1 √ + √ − = √ √ − 2x x x 2y y 2x x 2x y = 0 ∗ e = a ηxx + b ηxy + c ηyy + d ηx + e ηy 1 1 y 1 1 = − √ + √ − √ + √ 2x x x 2y y 2x x 2x y = 0 f∗ = f = 0 g∗ = g = 0 We now build the differential equation 2.1.20 in the new coordinates (ξ, η), y 2  1+ uξη = 0. −√ xy x (2.2.31) which simplifies to uξη = 0. The solutions of this equation are given by √ √ √ √ u = F (ξ) + G(η) = F [2( y − x)] + G[2( y + x)]. (c) x2 uxx − y 2 uyy − 2yuy = 0. We further assume that x 6= 0 and y 6= 0 otherwise the equation would not be hyperbolic. We can divide by x2 and find y  y 2 uxx − uyy − 2 uy = 0. (2.2.32) x x The coefficients are a = 1,

b = 0,

c=−

 y 2 x

,

d = 0,

e = −2

y , x2

f = 0,

g = 0.

2.2 Problems

23

We map the partial derivatives into total derivatives by using mapping 2.1.25 and find 

dy dx

2 −

 y 2 x

= 0.

(2.2.33)

The factorization of this quadratic equation This yields 

y dy − dx x



dy y + dx x

 = 0,

(2.2.34)

from which we find Z

Z dy dx − = ξ y x Z Z dx dy + = η y x (2.2.35)

so the canonical transformation is given by ξ = ln y − ln x η = ln y + ln x and from here 1 1 1 ξx = − , ξy = , ξxx = 2 , x y x 1 1 1 ηx = , ηy = , ηxx = − 2 , x y x

1 y2 1 = − 2. y

ξxy = 0 ξyy = − ηxy = 0 ηyy

and also x = e(η−ξ)/2 ,

y = e(η+ξ)/2

Let us verify that a∗ = b∗ = 0 1 y2 1 − =0 x2 x2 y 2 y2 1 1 2 2 = a ηx + 2b ηx ηy + cηy = 2 − 2 2 = 0 x x y

a∗ = a ξx2 + 2b ξx ξy + c ξy2 = c∗

(2.2.36)

24

second order The other coefficients are: b∗ = a ξx ηx + b (ξx ηy + ξy ηx ) + c ξy ηy y2 1 1 = − 2 2 xy x y   1 1 1 − = x y x d∗ = a ξxx + b ξxy + cξyy + dξx + eξy 1 y2 1 2y 1 = 2+ 2 2− 2 x x y x y = 0 ∗ e = a ηxx + b ηxy + c ηyy + d ηx + e ηy y2 1 2y 1 1 = − 2+ 2 2− 2 + x x y x y 2 = − 2 x f∗ = f = 0 g∗ = g = 0 Equation 2.1.20 turns out to be 2 2 (x − y)uξη − 2 uη = 0, 2 xy x after we multiply by x2 y/2 we find (x − y)uξη − yuη = 0, and from 2.2.36 (e(η−ξ)/2 − e(η+ξ)/2 )uξη − e(η+ξ)/2 uη = 0. We further multiply by e−η/2 eξ/2 and find (1 − eξ ) uξη − eξ uη = 0. Let us call v = uη so (1 − eξ ) vξ − eξ v = 0.

(2.2.37)

we consider this as an ordinary differential equation in ξ. To solve this equation we have two options

2.2 Problems

25

• ξ = 0 (this happens along the line y = x) This implies v = 0, so uη = 0, ⇒ u = F (ξ) = F (ln y − ln x) = F (0), That is the solution is a constant, u(x, y) = constant. • ξ 6= 0 We divide 2.2.37 by 1 − eξ vξ −

2eξ v = 0. 1 − eξ

(2.2.38)

Let us use the method integrating factor. Consider a function M (ξ). We multiply both sides of 2.2.39 by M (ξ) and get M (ξ) vξ − M (ξ)

eξ v = 0. 1 − eξ

(2.2.39)

We want the left hand side to be in the form of the derivative of a product, such that this equation can be written as (M (ξ)v)0 = 0,

(2.2.40)

which after integration becomes vM (ξ) = ϕ(η) ⇒ v =

ϕ(η) . M (ξ)

(2.2.41)

We apply the product rule to 2.2.40 and equating to the left hand side of 2.2.39 M 0 (ξ)v + M (ξ)vξ = M (ξ) vξ − M (ξ)

eξ v 1 − eξ

(note that v 0 = vξ ) so M 0 (ξ) eξ =− M (ξ) 1 − eξ and from here M (ξ) = e−

Rξ o

es /(1−es )

ds

(2.2.42)

26

second order where o is any constant. Now Z ξ es − d s = ln(1 − eξ ) + constant s 1 − e o so equation 2.2.42 becomes M (ξ) = C(1 − eξ )2

(2.2.43)

for a constant C. So, from 2.2.41 v=

ϕ(η) . (1 − eξ )

where now the constant C was absorbed into ϕ(η). Since v = uη we have uη =

ϕ(η) . 1 − eξ

and so, by integration 1 u= 1 − eξ

Z ϕ(η) + G(ξ).

and so u(ξ, η) =

1 F (η) + G(ξ) 1 − eξ

Finally, using 2.2.36 we see that u(x, y) =

1 F (ln y + ln x) + G(ln y + ln x) 1 − xy

(d) ∂ ∂x

  ∂ 2u 2 ∂u x = x2 2 . ∂x ∂y

We start by using the product rule of differentiation x2

2 ∂ 2u ∂u 2∂ u + 2x − x = 0. ∂x2 ∂x ∂y 2

Dividing by x2 and re–arranging we find ∂ 2 u ∂ 2 u 2 ∂u − + = 0. ∂x2 ∂y 2 x ∂x

2.2 Problems

27

This defines the coefficients a = 1,

c = −1,

b = 0,

d=

2 , x

e = f = g = 0.

(2.2.44)

The characteristics are given by solutions of  2 dy − 1 = 0. dx That is, the solutions of 

  dy dy −1 + 1 = 0. dx dx

Integrating, Z

Z dy =

Z dx + ξ

Z dy = −

dx + η

ξ = y−x η = y+x

(2.2.45) (2.2.46)

So, ξx = −1, ξy = 1, ξxx = 0, ξxy = 0 ξyy = 0 ηx = 1, ηy = 1, ηxx = 0, ηxy = 0 ηyy = 0. Let us verify that a∗ = b∗ = 0 a∗ = a ξx2 + 2b ξx ξy + c ξy2 = 1 − 1 = 0 c∗ = a ηx2 + 2b ηx ηy + cηy2 = 1 − 1 = 0 The other coefficients are: b∗ = a ξx ηx + b (ξx ηy + ξy ηx ) + c ξy ηy = −1 − 1 = −2 2 d∗ = a ξxx + b ξxy + cξyy + dξx + eξy = − x 2 e∗ = a ηxx + b ηxy + c ηyy + d ηx + e ηy = x f∗ = f = 0 g∗ = g = 0 Equation 2.1.20 turns out to be uξη −

1 1 uξ + uη = 0. 2x 2x

(2.2.47)

28

second order From equation 2.2.45 we find y=

η+ξ 2

x=

η−ξ 2

(2.2.48)

and equation 2.2.47 becomes uξη −

1 1 uξ + uη = 0. · · · η−ξ η−ξ

(2.2.49)

Problems with initial conditions 1. (a) uxx + 4uxy − 5uyy + ux − uy = 0,

u|y=0 = f (x),

∂u = F (x). ∂y y=0

We identify the coefficients a = 1,

b = 2,

c = −5,

e = −1,

d = 1,

f = 0,

g = 0.

The characteristics are given by the solution of the quadratic equation 

dy dx

2

 −4

dy dx

 − 5 = 0.

We factor this quadratic form into    dy dy −5 +1 dx dx which produces the following change coordinate transformations y − 5x = ξ y+x=η so ξx = −5, ξy = 1, ξxx = 0, ξxy = 0, ξyy = 0 ηx = 1, ηy = 1, ηxx = 0, ηxy = 0, ηyy = 0 first verify that a∗ = c∗ = 0. a∗ = a ξx2 + 2b ξx ξy + c ξy2 = 25 − 20 − 5 = 0 c∗ = a ηx2 + 2b ηx ηy + cηy2 = 1 + 4 − 5 = 0.

2.2 Problems

29

the other coefficients are b∗ d∗ e∗ f∗ g∗

= = = = =

a ξx ηx + b (ξx ηy + ξy ηx ) + c ξy ηy = −5 − 8 − 5 = −18 a ξxx + b ξxy + cξyy + dξx + eξy = −5 − 1 = −6 a ηxx + b ηxy + cηyy + dηx + eηy = 1 − 1 = 0 f =0 g = 0,

Equation 2.1.20 turns out to be −36uξη − 6uξ = 0. which is 6uξη + uξ = 0. We factor this equation as ∂ ∂ξ

  ∂u 6 + u = 0. ∂η

and so   ∂u + u = h(η), 6 ∂η for some h function only on η. We can look this as an ordinary differential equation du u + + h(η) = 0. dη 6 To solve this we multiply it by the integrating factor eη/6 eη/6

du u d(eη/6 u) + eη/6 + eη/6 h(η) = + eη/6 h(η) = 0. dη 6 dη

which can be integrated as u(ξ, η) = e−η/6 g(η) + p(ξ).

(2.2.50)

where Z g(η) = −

eη/6 h(η)dη.

To find the exact form of f and g we use the initial conditions. Let us transform 2.2.50 into the original variables x, y, before applying the initial conditions. That is

30

second order

u(x, y) = e−(x+y)/6 g(x + y) + p(y − 5x). From the first initial condition u(x, 0) = f (x) we find u(x, 0) = e−x/6 g(x) + p(−5x) = f (x),

(2.2.51)

and from the second condition du(x, 0) = F (x), dy y=0 we see 1 F (x) = − e−x/6 g(x) + e−x/6 g 0 (x) + p0 (−5x). 6

(2.2.52)

We should find g(x) and p(x) as functions of f (x) and F (x). From 2.2.51 taking the derivative with respect to x 1 f 0 (x) = e−x/6 g 0 (x) − e−x/6 g(x) − 5p0 (−5x) 6

(2.2.53)

From 2.2.51 and 2.2.53 into 2.2.52 1 1 F (x) = − [f (x) − p(−5x)] + f 0 (x) + [f (x) − p(−5x)] 6 6 0 0 0 +5p (−5x) + p (−5x) = f (x) + 6p0 (−5x). so 6 f (s) − p(−5s) = 5

s

Z

F (t)dt 0

or 6 p(−5s) = 5

  Z s f (s) − F (t)dt 0

and from equation 2.2.51 g(s) = e

s/6

  Z 5 5 s [f (s) − p(−5s)] = e f (s) − f (s) + F (t)dt 6 6 0   Z s es/6 = f (s) − 5 F (t)dt 6 0 s/6

and the solution for u is 1 u(x, y) = 6 6 5

!

−5(x+y)

Z f (x + y) − 5

F (s)ds

+

0

 f

5−y 5



Z −

((5−y)/5

! F (s)ds

0

2.3 Semi–infinite string

2.3 2.3.1

31

Semi–infinite string Mixed problem for the d’Alembert equation

The semi–infinite string can be modeled with the d’Alembert equation utt = a2 uxx ,

x > 0 t > 0,

where one end is located at x = 0 and the other is located at infinite (which for practical purposes is  at) and a represents the wavespeed. We add to this initial and boundary conditions. The initial conditions are given by u(x, 0) = ϕ(x),

∂u (x, 0) = ψ(x), ∂t

x > 0.

The boundary condition, for the string fixed at the left (x = 0) extreme should be a zero condition (since the string is fixed), that is u(t, 0) = 0,

t > 0.

(2.3.54)

The solution for the infinite string is given by equation 2.1.12, that is Z x+at ϕ(x + at) + ϕ(x − at) 1 u(x, t) = + ψ(s) ds. 2 2a x−at

(2.3.55)

This equation is valid as long as x > at. Since t, a, x > 0 then x + at > 0, however x − at could be negative and we have that ϕ is not defined in the negative real line. Therefore we have to modify this solution. The general solution of the equation 2.3.54 is given by u(x, t) = f (x − at) + g(x + at).

(2.3.56)

We know that at x = 0, u(0, t) = f (−at) + g(at) = 0,

t > 0.

(2.3.57)

so we can have the function f in its negative real line, defined as a function of g which is defined in the positive real line. So, we need to find g and that is what we do now. From the initial conditions (t = 0) u(x, 0) = f (x) + g(x) = ϕ(x)

(2.3.58)

∂ u(x, 0) = −a f 0 (x) + a g 0 (x) = ϕ(x) ∂t

(2.3.59)

from the second of these equations, after integration Z x −af (x) + ag(x) = ψ(s)ds + c. 0

32

second order

which when added with the first of equations 2.3.58 (scaled by a) yields Z x ψ(s)ds + c. 2 a g(x) = a ϕ(x) + 0

so ϕ(x) 1 g(x) = + 2 2a

Z

x

ψ(s)ds + 0

c . 2a

(2.3.60)

From 2.3.57, we see that for negative argument, −s < 0 f (−s) = −g(s), so for x − at < 0, Z a t−x ϕ(at − x) 1 c f (x − at) = −g(at − x) = − − ψ(s)ds − 2 2a 0 2a Z 0 ϕ(x − a t) 1 c = − + ψ(s)ds − . 2 2 a a t−x 2a

(2.3.61)

and then for 0 < x < a t, the solution is given by u(x, t) = f (x − a t) + g(x + a t) Z x+a t 1 ϕ(x + at) − ϕ(a t − x) + ψ(s)ds = 2 2 a a t−x

2.3.2

Other boundary conditions

Instead of specifying the displacement data at x = 0 (Dirichlet), we can specify the derivative (Neumann) boundary conditions. Particularly ∂u (0, t) = 0, ∂x

t > 0.

(2.3.62)

Below the principal characteristic curve x > at, equation 2.3.55 is still valid. Let us consider the case above the principal characteristic curves, that is, for x < at. We start with the general equation 2.3.56, which after substituting the boundary condition 2.3.62 is f 0 (−at) + g 0 (at) = 0,

t > 0.

We use the substitution z = −at so f (z) + g 0 (−z) = 0,

z < 0.

Integrating, we get f (z) − g(−z) = c1 = const,

z < 0,

2.3 Semi–infinite string

33

so from 2.3.60 , and for x < at 1 1 f (x − at) = g(at − x) + c1 = ϕ(at − x) + 2 2a

Z

at−x

ψ(s)ds + 0

c + c1 2a (2.3.63)

and using 2.3.60 again we find 1 ϕ(at − x) + ϕ(x + at) + u(x, t) = 2 2a

Z

at−x

ψ(s)ds + c2 . 0

The constant c2 could be found continuity. That is, we assume that u(x, t) is continuous at the characteristic curve x = at. Continuity is necessary for a string or a rod, but it is not for gases where shock waves could be present.

2.3.3

Discontinuities of a solution along a principal characteristic curve. Continuity conditions

In the previous section we solved the initial value and boundary at one side, one– dimensional (string) wave problem. We found that the solution is expressed as two branches. One branch for x > at and another for x < at. We wonder about the continuity along the characteristic x = at. From the general solution u(x, t) = f (x − a t) + g(x + a t) Z x+a t 1 ϕ(x + at) − ϕ(a t − x) + ψ(s)ds = 2 2 a a t−x The wave g(x + at) is continuous when passing through the principal characteristic curve, since its level curves x + at = const intersect the line x = at. The wave f (x − at) below the principal characteristic curve x − at = 0 has a limit equal to f (0−). Thus, u − u = f (0−) − f (0+). x−at=0−

x−at=0+

where f (a±) := lim f (x). x→a±0

A solution u(x, t) is continuous on the principal curve characteristic if f (0−) = f (0+).

(2.3.64)

We now find, under which conditions the semi–infinite string problem has continuous solutions through the principal characteristic curve x − at = 0. Since g is continuous, we focus in the f component. Let us recall the initial conditions u(x, 0) = ϕ(x) ∂u(x, t) = ψ(x). ∂t t=0

34

second order

which imply f (x) + g(x) = ϕ(x) −f 0 (x)a + g 0 (x)a = ψ(x)

(2.3.65) (2.3.66)

integrating the second of equations 2.3.65 1 −f (x) + g(x) = a

Z 0

x

c ϕ(s)dx + . a

Subtracting this, from the first of equations 2.3.65, we find, after scaling by 1/2, Z x 1 1 c f (x) = ϕ(x) − ψ(s)ds − , 2 2a 0 2a so f (0+) =

c ϕ(0) − , 2 2a

(2.3.67)

and from 2.3.61 f (0−) = −

c ϕ(0) − , 2 2a

Therefore, continuity is achieved if −

ϕ(0) ϕ(0) = , 2 2

hence ϕ(0) = 0.

(2.3.68)

Due to the boundary condition u(0, t) = 0, t > 0 we have that making u(0, 0) = 0 is imposing continuity at the point (0, 0). Let us now study the Neumann boundary condition 2.3.62 ∂u (0, t) = 0, ∂x

t > 0.

and how continuity is related to this condition. We use the relation 2.3.63, and find f (0−) =

ϕ(0) c + + c1 . 2 2a

From, here, equation 2.3.67 and the continuity condition 2.3.64 we see that ϕ(0) c ϕ(0) c + + c1 = − , 2 2a 2 2a

hence c1 +

c = c2 = 0. a

In a string or a road, discontinuity c2 6= 0 does not make sense since this implies the braking of the string or the road. In gas dynamics discontinuity is achieved along the principal characteristics and this is known in the fluid dynamics literature as shock waves. This type of discontinuity needs additional conditions not considered on this document.

2.3 Semi–infinite string

35

Example: Find a continuous solution to the problem  utt = 9uxx x > 0, t > 0; −x u(x, 0) = e , ut (x, 0) = cos 5x; ux (0, t) = u(0, t) + t. We identify ϕ(x) = e−x ψ(x) = cos 5x. and so from the general solution for x ≥ 3t 2.1.12 we have Z e−(x−3t) + e−(x+3t) 1 x+3t u(x, t) = + cos 5s ds 2 6 x−3t e−(x−3t) + e−(x+3t) 1 sin[5(x + 3t)] − sin[5(x − 3t)] = + . 2 6 5 (2.3.69) Let us now assume x < 3t. The general solution for u is given by u(x, t) = f (x − 3t) + g(x + 3t) Let us first find g(s). From the initial conditions u(x, 0) = e−x = f (x) + g(x) x > 0 ut (x, 0) = cos 5x = −3f 0 (x) + 3g 0 (x) x > 0

(2.3.70)

From the second equation f (x) − g(x) = −

sin 5x 15

and subtracting this from the first of equations 2.3.70 g(x) =

e−x sin 5x + . 2 30

We can not obtain the function f (x) as we did for g(x) because, the argument of f in u(x, t) is x − 3t which is negative for x < 3t, and we have that f is not defined (see for example equation 2.3.70) for negative argument. Hence we have u(x, t) = f (x − 3 t) + g(x + 3 t) e−(x+3t) sin 5(x + 3t) = f (x − 3 t) + + . 2 30 We use the boundary condition in u, that is f 0 (−3t) −

e−3t cos 15t e−3t sin 15t + = f (−3t) + + + t, 2 6 2 30

t > 0.

36

second order

Before solving this ordinary differential equation we make the change of variables y = −3t.

f 0 (y) − f (y) − ey +

cos 5y sin 5y y + + = 0, 6 30 3

y < 0,

multiplying by the integrating factor e−y e−y cos 5y e−y sin 5y e−y y + + , e f (y) − e f (y) − 1 + 6 30 3 −y 0

−y

y < 0,

so [e−y f (y)]0 − 1 +

e−y cos 5y e−y sin 5y e−y y − + = 0, 6 30 3

y < 0,

which we integrate   Z e−y cos 5y e−y sin 5y e−y y −y e f (y) = dy 1 − + − + c, 6 30 3

y < 0,

for some constant c. That is 1 −y e−y e (cos 5y − 5 sin 5y) − (sin 5y − 5 cos 5y) 156 780 e−y (y + 1) + c, + 3

e−y f (y) = y +

so f (y) = ey y +

1 1 y+1 (cos 5y − 5 sin 5y) − (sin 5y − 5 cos 5y) + + cey 156 780 3

and after simplifying f (y) = ey y + c ey +

y+1 1 1 + cos 5y − sin 5y 3 78 30

We now use the continuity condition f (0−) = f (0+) For f (0+) we use, from 2.3.69, the partition with the x − 3t argument, that is f (y) = and for f (0−) we use 2.3.71.

e−y sin 5y − . 2 30

(2.3.71)

2.3 Semi–infinite string

37

c+

1 1 1 + = , 3 78 2

so c=

1 1 1 2 − − = , 2 3 78 13

and so 1 2 x−3t x − 3t + 1 e + + cos 5(x − 3t) 13 3 78 1 e−(x+3t) sin 5(x + 3t) − sin 5(x − 3t) + + . 30 2 30

u(x, t) = ex−3t (x − 3t) +

38

second order

Chapter 3 Green Functions 3.1 3.1.1

Problems The sifting property of a delta distribution

Show that hδ(x − x0 ), φ(x)i = φ(x0 ). We show this formally using distribution and then present a heuristic view through delta sequences. sln. Proof using distributions Z

∞

hδ(x − x0 ), φ(x)i =

Z

∞

δ(x − x0 )φ(x)dx = −∞

−∞

δ(y)φ(y + x0 ) = φ(y + x0 )

= φ(x0 ). y=0

Note that we use the change of variable (substitution) on integration y = x − x0 Proof using delta sequences Let us assume that fn is a delta sequence around x0 . That is for each n ∈ N: 1. fn (x) ≥ 0,

x ∈ R,

2. Z

∞

fn (x)dx = 1, −∞

(3.1.1)

40

Green Functions 3. Z

b

fn (x)dx = 1,

lim

n→∞

for any a < 0 < b.

(3.1.2)

a

We take the inner product with a test function φ(x). That is Z

∞

fn (x − x0 )φ(x)dx

hfn (x − x0 ), φ(x)i = Z−∞ ∞ =

fn (y)φ(y −∞ Z −1/n

=

+ x0 )dy Z

1/n

fn (y)φ(y + x0 )dy + −∞

fn (y)φ(y + x0 )dy −1/n

∞

Z +

fn (y)φ(y + x0 )dy 1/n

From equations 3.1.1 and 3.1.2 the first and third integrals above go to zero as n → ∞. The second integral, after using the mean theorem for integrals is Z

1/n

Z

1/n

fn (y) = φ(y0 + x0 )

fn (y)φ(y + x0 )dy = φ(y0 + x0 ) −1/n

−1/n

as n → ∞, since −1/n < y0 < 1/n, we see that y0 → 0, and from equation 3.1.2, Z

1/n

fn (y)dy = 1 −1/n

and so lim hfn (x − x0 ), φ(x)i = φ(x0 ).

n→∞

3.1.2

Problem 23.7 c

Solve x2 u00 (x) + 2xu0 (x) = f (x), u0 (1) = 0, u(2) + 5u0 (2) = 0,

sln.

1<x 0, then we are integrating in a domain outside of the support of the Dirac delta which is the source function

68

Higher dimensions

where we rename the constant of integration as c1 = dE/dx(r) and rewrite dE (r) = Θ(r) + c1 dx We proceed in exact the same way and do one more integration to get E(r) = r Θ(r) + c1 r + c2 .

(4.1.31)

There is an infinite number of fundamental solutions. The only way to constrain the solution to be unique is by imposing initial condition and/or radiation conditions (conditions at ∞). Not all solutions have the form shown in equation 4.1.31. Since for each solution E(r), E(−r) is also a solution, then the equation E(−r) = −rΘ(−r) − c1 r + c2 is also a solution. In particular if c1 = c2 = 0 then E+ (r) = r Θ(r) and E− (r)

− rΘ(−r)

are solutions and so d2 E + d2 E − + = 2δ(r), dr2 dr2 Now rΘ(r) − rΘ(−r) |x| E+ (r) + E− (r) = = 2 2 2 and so E(r) =

E+ (r) + E− (r) |x| |r| = = , 2 2 2

is also a fundamental solution of u00 (r) = δ(r), up to any linear shift. The Laplacian equation in n > 1 dimensions We should solve   d n−1 dE δ(r) ∇ E(r) = n−1 r = . r dr dr Sr 2

1

We use the solution of the homogeneous problem 4.1.27, which is valid for all r > 0 and find appropriate c1 and c2 such that in the limit as r → 0 the solution is fundamental. To do this we integrate the Laplace equation 4.1.23 along a ball of radius r around 0. That is Z Z 2 ∇ E(x) = δ(x) = 1 kxk≤r kxk≤r

4.1 Problems

69

then we apply the divergence theorem on ∇2 E(x) = ∇ · ∇E(x) and find Z kxk=1

dE(x) dS = 1. dn

where d/dn is the normal derivative out of the surface. That is, the derivative with respect to r. We then have Z dE(x) dS = 1. dr kxk=1 We know, from equation 4.1.26 (for r > 0) that (n − 2) c1 dE0 = dr rn−1 so we set up the equation Z r=1

c1 drdS n−1 r

=

c1 Sn (r) n−1 r

= 1,

and from 4.1.10 Γ n+2 1 2 c1 = = n/2 Sn (r) nπ




So back to problem 4.1.27, we rewrite   c1 ln r + c2 n = 2 E0 (r) =  c1 + c2 n>2 rn−2 and make c2 zero for convenience 2 . We found  1   2π ln r E(r) =  Γ( n  2)

2(n−2)π n/2 rn−2

n=2 n>2

In particular, for n = 4 we find E(r) = 2

1 . 4π 2 r2

solutions should have finite energy, then c2 = 0, this is a radiation condition to infinity

70

Higher dimensions

4.1.3

The Helmholtz equation

Find the fundamental solution of the Helmholtz equation in three dimensions. That is, find E(x) such that ∇2 E(x) ± k 2 E(x) = δ(x). We transform the equation into polar spherical coordinates. A lot of the work is already done because the Laplacian is most of what we have in the Helmholtz equation and we already did the transformation for the Laplacian and for δ(x) in section 4.1.2. Now, k 2 is constant, so the transformation to any coordinate system will remain as k 2 .   δ(r) 1 d n−1 dE 2 r ± k 2 E(r) = . ∇ E(r) = n−1 r dr dr Sr In particular for 3 dimensions we have   1 d 2 dE δ(r) 2 ∇ E(r) = 2 r ± k 2 E(r) = . r dr dr 4πr

(4.1.32)

We proceed in two steps. First we solve the homogeneous problem, which is valid everywhere for r > 0, in terms of two constants c1 and c2 . Then by integrating 4.1.32 we find the constants c1 and c2 . 3 The homogeneously problem which we rewrite equation 4.1.32 as  2 0 0 r E (r) ± k 2 r2 G(r) = 0

r > 0.

This is a Bessel function and one trick to solve it is to make the change of variables E=

f r

So we find  2 0 0  0 k 2 f r2 r E (r) ± k 2 G(r) = r2 (f /r)0 ±  0  r 00 f f = 2r + r2 ± k2f r r r  0 0 2r 0 2rf f f 2 = f − 2 +r − 2 ± k2f r r r r r  00  0 2f f f f 0 2f 0 2 = 2f − +r − 2 − 2 + 3 ± k2f r r r r r r 2f 2f = 2f 0 − + rf 00 − 2f 0 + ± k2f r = 0 r r 3

And assuming finite energy physical principles

4.1 Problems

71

From which the simpler equation results (in terms of f ) (f 00 ± k 2 )f = 0. The obvious solutions for this homogeneous equations are • Complex exponential functions (+k 2 ). The two solutions here are given by c± e±ikr E(r) = 4πr Constants of integration The constants of integration c± are chosen from the two solutions. The solution with the “minus” sign is picked, since it represents an outgoing wave (causal). By convenience we can pick c− = 1, and so E(r) =

e−ikr 4πr

r > 0.

• Real exponential functions. (−k 2 ) If instead we choose to solve (f 00 − k 2 )f = 0. The solutions are real exponentials. That is E(r) =

c± e−±kr 4πr

r > 0.

To have finite energy we pick the negative exponent, and for convenience c− = 1. So E(r) =

4.1.4

e−kr 4πr

r > 0.

The Cauchy–Riemann equation

Find the fundamental function of the equation ∂ u ∂ z¯ Let us first give meaning to the expression ∂/∂ z¯. We assume z = x + iy, with x, y real √ numbers and i = −1. Then for f (z) = f (x + iy),     ∂f ∂f ∂f 1 ∂f i df = dx + dy = (dz + d¯ z) + − ( dz − d¯ z) ∂x ∂y ∂x 2 ∂y 2     1 ∂f ∂f 1 ∂f ∂f = −i dz + +i d¯ z. 2 ∂x ∂y 2 ∂x ∂y

72

Higher dimensions

from which we see that 1 ∂f = ∂z 2



∂f ∂f −i ∂x ∂y



∂f 1 = ∂ z¯ 2



∂f ∂f +i ∂x ∂y

 .

Then dz dz dz d¯ z d¯ z dz d¯ z d¯ z

= = = =

  1 1 + 0 − 0 − (−1) = 1 2   1 1+0+0−1 =0 2   1 1+0−0−1 =0 2   1 1 + 0 + 0 − (−1) = 1. 2

We observe that ∂ ∂ 1 f= ∂ z¯ ∂z 4



∂2 ∂2 + ∂x2 ∂y 2



1 f = ∇2 f. 4

(4.1.33)

We can take advantage of knowing that in 2D the fundamental solution of the Laplacian is ln r ∇2 E(x, y) = δ(x, y) with E(x, y) = 2π p with r = z z¯ = x2 + y 2 . Then ln r = ln(z z¯)1/2 =

1 ln z z¯ 2

and ∂z z¯ ∂(x2 + y 2 ) 1 = = ∂z ∂z 2



∂(x2 + y 2 ) ∂(x2 + y 2 ) −i ∂x ∂y

Now if f (x, y) = E(x, y) then ∂f 1 1 1 1 = z¯ = ∂z 2π 2 z z¯ 4πz. From 4.1.33 ∂ ∂ 1 δ(x, y) f = ∇2 f = . ∂ z¯ ∂z 4 4 Therefore by setting u(x) = 4

∂f 1 = , ∂z πz

we found the fundamental solution of ∂/∂ z¯.

 = z¯.

Appendix A The generalized polar–spherical coordinates .

A.1

Derivation

The natural generalization of the polar coordinates for 2D and spherical coordinates for 3D can be built upon the concepts of linear algebra. Figure A.1 shows an sketch of the polar–spherical coordinates in 3D. Let us then assume that we have an n–dimensional euclidean space En , and a basis {ei }, i = 1, 2, · · · , n. A vector x can be written in these basis as x=

n X

xi e i .

(A.1.1)

i=1

The idea is to find xi in terms of certain polar angles φi , i = 0, · · · n − 2 and the vector radius r = kxk. There is not a unified convention for the notation of the n parameters used, except for the radius r (...and not even here; some authors use ρ for the radius depending on the dimension). I will use the following convention for the angles. The first parameter introduced is the angle θ that the vector makes with the en axis. Note that in the 2D case, this parameter is taken with respect to the first axis, but after that, θ is used for the polar angle, the angle that the radio vector makes with the en axis. In 3D the parameters are as shown in figure A.1. That is the polar angle θ and the azimuthal angle φ. The starting point is in 3D. The 2D case can be seen embeded in 3D when the polar angle φ0 = θ = 0. What about 4D and on? We introduce the azimuthal angles φi angles one–by–one, in reverse order, and to make the notation more homogeneous define φ0 = θ. That is, first φ0 = θ angle is the angle between the vector x and the en axis. Since

74

The generalized polar–spherical coordinates

Figure A.1: Polar–Spherical Coordinates in 3D e3

θ (x1 , x2 , x3 )

φ

e2

e1

for each xi xi = hx, ei i. then xn = hx, en i = r cos φ0 = r cos θ and we can write A.1.1 as

x = r cos φ0 en +

n−1 X

xi ei ,

(A.1.2)

i=1

0 ≤ φ0 ≤ π. This restriction guarantees the unique value of the xn coordinate as a function of the angle φ0 , since the cosine function is one–to–one in this interval. The angle φ1 will be the angle between the projection of x into the hyper–space span{en }⊥ and the axis en−1 . This is the new “polar” angle in the n − 1 projected dimensional space. If we call the unit projection vector above as u1 then Pn−1 i=1 xi ei u1 = q P n−1 2 i=1 xi

A.1 Derivation

75

The angle between en−1 and u1 is cos φ1 , that is cos φ1 = hu1 , en−1 i xn−1 = qP n−1 2 i=1 xi xn−1 = p r2 − x2n xn−1 = p r2 − r2 cos2 φ0 xn−1 = p r 1 − cos2 φ0 xn−1 = , r sin φ0 so xn−1 = r sin φ0 cos φ1 , where 0 ≤ φ1 ≤ π. We rewrite A.1.2 as x = r cos φ0 + r sin φ0 cos φ1 +

n−2 X

xi e i .

i=1

Following the same inductive idea, we define φ2 as the angle of the projection of x into the orthogonal complement of the hyper–plane spanned by the set {en , en−1 }, and the vector en−2 . A unit vector in this hyperplane is written as Pn−2 i=1 xi ei , u2 = qP n−2 2 i=1 xi and the angle between x and u2 is found by cos φ2 = hu2 , en−2 i xn−2 = qP n−2 2 i=1 xi xn−2 = p r2 − x2n − x2n−1 xn−2 = p r2 − r2 cos2 φ0 − r2 sin2 φ0 cos2 φ1 xn−2 = p 2 r sin φ0 − sin2 φ0 cos2 φ1 xn−2 = , r| sin φ0 | sin φ1

76

The generalized polar–spherical coordinates

Since, sin φ0 ≥ 0, for 0 ≤ φ0 ≤ π, we have then

xn−2 = r sin φ0 sin φ1 cos φ2 ,

When the projection of the vector x falls into a three dimensional space, the expression for x3 inductively would be

x3 = r sin φ0 sin φ1 sin φ2 · · · sin φn−4 cos φn−3

with all angles in the interval [0, π] except for φ1 which satisfies 0 ≤ φ1 ≤ 2π. That is,

n−j−1

xj = r cos φn−j

Y

sin φi

,

j = 3 · · · n.

i=0

The last step is to project x into the plane spanned by axis {e1 , e2 }. That is,the unit vector projection is x1 e1 + x2 e2 un−2 = p 2 x1 + x22 and so x2 sin φn−2 = hun−2 , e2 i = p 2 x1 + x22 (A.1.3) x1

cos φn−2 = hun−2 , e1 i = p 2 x1 + x22

At this point we note that the angle φn−2 can be in the interval [0, 2π] and with the radius r there is a unique mapping between (x1 , x2 ) and ρ, φn−2 , where ρ is the length of the projected x vector along the x1 − x2 plane.

A.1 Derivation

77

Now, x21 + x22 = r2 − x2n − x2n−1 − · · · − x23  2 = r 1 − cos2 φ0 − cos2 φ1 sin2 φ0 − cos2 φ2 sin2 φ0 sin2 φ1 − · · · 2

− cos φn−3

n−4 Y



2

sin φi

i=0



= r2 sin2 φ0 − cos2 φ1 sin2 φ0 − cos2 φ2 sin2 φ0 sin2 φ1 − · · · 2

− cos φn−3

n−4 Y



2

sin φi

i=0 2

2



2

2

2

2

= r sin φ0 1 − cos φ1 − cos φ2 sin φ1 − · · · − cos φn−3

n−4 Y

2



sin φi

i=1 2

2



2

2

2

2

= r sin φ0 sin φ1 − cos φ2 sin φ1 − · · · − cos φn−3

n−4 Y

2



sin φi

i=1 2

2



2

2

2

= r sin φ0 sin φ1 1 − cos φ2 − · · · − cos φn−3

n−4 Y

2



sin φi

i=2

.. . = r

2

n−3 Y

sin2 φi ,

i=0

so q

x21

+

x22

=r

n−3 Y

sin φi ,

i=0

and from A.1.3 x1 = r

n−2 Y

sin φi ,

x2 = r cos φn−2

i=0

n−3 Y

sin φi ,

i=0

In summary

xj =

 Qn−j−1  r cos φn−j i=0 sin φi 

Q

r

Qn−2 i=0

sin φi

, j = 2···n j=1

The product is assumed to be 1 if the upper index is smaller than the lower index. The angle θ is the polar angle between 0 and π, the last azimuthal angle φn−2 between

78

The generalized polar–spherical coordinates

0 and 2π and rest of the angles φi are polar angles in their own domain, between 0 and π. We can rewrite equation A.1.4 in a more compact form by calling φn−1 = 0. n−j−1

xj = r cos φn−j

Y

sin φi

,

j = 1···n

(A.1.4)

i=0

I have seen notations that are in the reverse order with respect to the notation here. See Wikipedia, for example.

A.2

The Jacobian

The Jacobian is, by definition, the determinant ∂x1 ∂x1 dr dφ · · · ∂x ∂x0 2 2 dr dφ0 · · · J = . .. .. .. . ∂x ∂x. n n ··· dr dφ0

∂x1 dφn−2 ∂x2 dφn−2

.. .

∂xn dφn−2



In order to compute the Jacobian J we need to estimate the partial derivatives. This is done in the next section.

A.2.1

Partial Derivatives

We split the new variables in two groups. • The variable r. From equation A.1.4 ∂xj ∂r • The angle φk ,

n−j−1

= cos φn−j

Y

sin φi

j = 1···n

,

i=0

0≤k ≤n−2            

n−j−1

−r sin φn−j

Y

sin φi

if

k =n−j

i=0

∂xj n−j−1 = Y ∂φk   r cos φ cot φ sin φi if n−j k    i=0       0 if

k n−j

A.2 The Jacobian

79

For the differential matrix we will understand a matrix where its j column is given by the partial derivatives (xi ),k , for xi = 1 · · · n, k = 1, n . Here the value k = 1 corresponds to the radius r, the value k = 2 corresponds to the theta angle φ0 , and for any other value k > 2, the variable of differentiation is φk−2 . Here is the differential matrix, (recall that cos φn−1 = 1) 

n−2 Y  sin φi   i=0    n−3  Y  cos φ sin φi n−2   i=0     n−4 Y  sin φi D=  cos φn−3  i=0    ..   .     sin φ0 cos φ1     cos φ0

r cos φ0

n−2 Y

···

sin φi

i=1

r cos φ0 cos φn−2

n−2 Y

r cos φn−3

n−3 Y

···

sin φi

n−4 Y

n−3 Y

r cos φn−3 cos φn−2

sin φi

i6=n−3

···

sin φi

−r sin φn−3

n−4 Y

sin φi

i=0

i=1

.. .

..

.

.. .

r cos φ0 cos φ1

···

0

−r sin φ0

···

0

n−2 Y

The next step is to find det D. First, let us see that r is in each column after the second, and so by properties of the determinant det D = rn−1 det D0 where D0 is the same matrix D except that it does not have any r factor. Qn−3 Qn−2 Next we pull out Q0 i=0 sin φi from the first raw, i=0 sin φi from the second raw, and so forth until i=0 sin φi from the n − 1–th row. Then we get det D = r

n−1

n−2 Y

sin φi

i=0

n−3 Y

sin φi · · ·

i=0

0 Y

sin φi det G

(A.2.5)

i=0

where the matrix G is represented as follows:            G=         

1

cot φ0

···

cos φn−2 cot φ0 cos φn−2 · · · cos φn−3 cot φ0 cos φn−3 · · ·

cot φn−3

cot φn−2

cot φn−3 cos φn−2 − sin φn−2 − sin φn−3

0

.. .

.. .

..

.

.. .

.. .

cos φ1

cot φ0 cos φ1

···

0

0

cos φ0

− sin φ0

···

0

0



sin φi    i6=n−2    n−3  Y sin φi  −r sin φn−2   i=0       0      ..   .     0     0

r cos φn−2

i6=n−3

i=1

r cos φ0 cos φn−3

sin φi

                    

80

The generalized polar–spherical coordinates

We, working on G, now pull out cot φ0 from the second column, cot φ1 from the third column and so up to cot φn−2 from the last (n–the) column to obtain

det G =

n−2 Y

cot φi det G1

(A.2.6)

i=0

with



1 1  cos φn−2 cos φn−2    cos φn−3 cos φn−3 G1 =  .. ..   . .   cos φ1 cos φ1 2φ 0 cos φ0 − sin cos φ0

··· ···

1 cos φn−2 2

φn−3 · · · − sin cos φn−2 .. ... . ··· 0 ··· 0

1



2φ n−2  − sin cos φn−2    0  ..   .   0 0

Now, working on G1 we pull out cos φn−2 from the second raw, cos φn−3 from the third raw, and so forth until cos φ0 from the last raw to find

det G1 =

n−2 Y

cos φi det H1 ,

(A.2.7)

i=0

with



1 1 1 .. .

1 1 1 .. .

··· 1 1 ··· 1 − tan2 φn−2 · · · − tan2 φn−3 0 .. .. .. . . .

    H1 =     1 1 ··· 2 1 − tan φ0 · · ·

0 0

0 0

        

A.2 The Jacobian

81

To evaluate det H we use Gaussian elimination. That is   tan2 φn−2 tan2 φn−2 · · · tan2 φn−2 tan2 φn−2  1 1 ··· 1 − tan2 φn−2    2   1 1 · · · − tan φ 0 n−3   2 det H1 = cot φn−2 det   .. .. .. .. . .   . . . . .     1 1 ··· 0 0 2 1 − tan φ0 · · · 0 0   tan2 φn−2 tan2 φn−2 · · · tan2 φn−2 tan2 φn−2  sec2 φn−2 sec2 φn−2 · · ·  sec2 φn−2 0   2   1 1 · · · − tan φ 0 n−3   = cot2 φn−2 det   .. .. . . . . . .   . . . . .     1 1 ··· 0 0 2 1 − tan φ0 · · · 0 0   1 1 ··· 1 1  1 1 · · · 1 0    2   1 1 · · · − tan φ 0 n−3   = sec2 φn−2 det  .. .. .. ..  ..  . . . . .     1 1 ··· 0 0  1 − tan2 φ0 · · · 0 0 = (−1)n sec2 φn−2 det H2 with H2 the n − 1 × n − 1 matrix 

1 1 .. .

1 1 .. .

··· 1 · · · − tan2 φn−3 .. .. . .

   H2 =    1 1 ··· 2 1 − tan φ0 · · ·

0 0

    .  

We now that H2 has the same structure as H1 , but instead of being an n × n matrix it is an n − 1 × n − 1, so we have the following recursion det H1 = (−1)n sec2 φn−2 det H2 = (−1)n (−1)n−1 sec2 φn−2 sec2 φn−3 det H3 = (−1)1 sec2 φn−2 sec2 φn−3 det H3 .. . n−2 Y = (−1)n−1 sec2 φi det Hn−1 , i=1

with  Hn−1 =

1 1 1 − tan2 φ0

 ,

82

The generalized polar–spherical coordinates

and det H1 = (−1)

n−2

n−2 Y

sec2 φi .

(A.2.8)

i=0

We will ignore the sign in what follows. The sign is related to, if the system is right or left hand oriented. I assume it is right hand oriented and ignore the “minus” sign. We are now ready to put all the pieces together, by invoking equations A.2.5, A.2.6, A.2.7, and A.2.8. That is det D = rn−1 = r = r = r

n−2 Y

n−1

i=0 n−2 Y

n−1

i=0 n−2 Y

n−1

i=0 n−2 Y i=0

sin φi sin φi sin φi sin φi

n−3 Y i=0 n−3 Y i=0 n−3 Y i=0 n−3 Y

sin φi · · · sin φi · · · sin φi · · · sin φi · · ·

i=0

0 Y i=0 0 Y i=0 0 Y i=0 0 Y i=0

sin φi det G sin φi sin φi sin φi

n−2 Y i=0 n−2 Y i=0 n−2 Y

cot φi det G1 , cot φi cot φi

i=0

n−2 Y i=0 n−2 Y i=0

cos φi det H1 cos φi

n−2 Y

sec2 φi

i=0

Q0 Qn−2  Qn−3  n−1  i=0 sin φi i=0 sin φi · · · i=0 sin φi = r Qn−2   i=0 sin φi  n−3 0 Y Y n−1 = r sin φi · · · sin φi = rn−1

i=0 n−2 Y

i=0

sinn−i−1 φi−1 .

i=1

A.3

The Volume and Surface Area

With the help of the Jacobian we can now compute the volume and the surface of a sphere of radius a. The formula for the volume is given by Z

Z

dV = rn−1 sinn−2 φ0 sinn−3 φ1 · · · sin φn−3 sin φ1 dr dφ0 dφ1 · · · dφn−2 SZ S = 2π rn−1 sinn−2 φ0 sinn−3 φ1 · · · sin φn−3 sin φ1 dr dφ0 dφ1 · · · dφn−3 ,

Vn =

S

Qn−3

with S = [0, a] i=0 [0, π/2] Given that both r and sin φi are bounded functions in the integration domains, we can apply Fubini’s rule and separate the integrals in decoupled variables as follows:

A.3 The Volume and Surface Area

Z

a n−1

Z

π

83

Z

n−2

π n−3

π

Z

sin φn−3 dφn−3 sin φ0 dφ0 sin φ1 dφ1 · · · r dr 0 0 0 Z π Z π Z π 2π r = sinn−2 φ0 dφ0 sinn−3 φ1 dφ1 · · · sin φn−3 dφn−3 n 0 0 0

Vn = 2π

0 n

Let us test this equation with simple cases. In 2D, n = 2 and V2 =

2πa2 = πa2 . 2

In 3D 2πr3 V3 = 3

Z 0

π

2πa3 4πa3 sin φ0 dφ0 = 2 = . 3 3

For higher dimensions the problem is more interesting. We have the recursion Z Z π Z π 2π rn+1 π n−1 n−2 Vn+1 = sin φ0 dφ0 sin φ1 dφ1 · · · sin φn−2 dφn−2 n+1 0 0 0 Z π nr = sinn−1 φ0 dφ0 n+1 0   Z Z π Z π 2πrn π n−2 n−3 sin φ1 dφ1 sin φ1 dφ1 · · · sin φn−2 dφn−2 n 0 0 0 Z π nr = Vn sinn−1 φ0 dφ0 , n+1 0 and after using equation 4.1.4 where Z π/2 Z 1 π n−1 n−1 sin φ0 dφ0 = sin φ0 dφ0 , an−1 = 2 0 0 we find the recursion: Vn+1

nr(2 an−1 ) nr Γ n2 = Vn = n+1 n+1 Γ




nr Γ n2 = n+1 Γ




Γ

n+1 2

Γ

n+1 2

1
2




1
2




Vn (n − 1)r Γ n

n−1 2

Γ




Γ
n

1 2


Vn−1

2

where instead of writing a we use the symbol r for the radius. To easy on this recursion we will separate member factors as follows. • We will stop at V2 = πr2 . • The r factors will go up to rn−1

84

The generalized polar–spherical coordinates • The fractions with n–like factors will be arranged as follows n − 2 n − 1 2  ···  n  n+1  n −1 3
n 

=

2 . n+1

• Finally, the Gamma function factors go like
1





n n−1 n−1 1 Γ 1 22Γ 12   Γ Γ Γ Γ Γ  2 2 2 2 2  



··· = n+1 3 n  Γ n+1 Γ Γ Γ 2 2  2 2 Putting all pieces togheter we find Vn+1 = πr

n+1



2πrn+1 Γn−1 21 2 Γn−1 21
=
n + 1 Γ n+1 (n + 1)Γ n+1 2 2

or Vn

πrn Γn−2 12
= n Γ n2 2 =




π n/2 rn
Γ n+2 2

which agrees with equation 4.1.5 for the volume of an n–dimensional ball. The surface area is easily found by taking the derivative of the volume formula. So the n–dimensional surface area formula is Sn =

nrn−1 π n/2 dVn (r)
= dr Γ n+2 2

which is a repeat of what was done in equation 4.1.10.

Appendix B Derivations of wave equations Which waves? Waves are everywhere. Waves are acoustic or elastic or electromagenetic etc. Waves are studied in 1D, 2D, 3D and n–dimensional spaces. Waves are studied in homogeneous/hetereogeneous, isotropic/anisotropic media. Waves carry energy and the manifestation of that energy is motion of particles which is registred as for example, particle displacement, particle velocity along some direction, pressure, acceleration along certain direction, voltage etc. Among the characteristics of wave propagation are frequency, wave speed, wavelength, phase etc. These characteristics depend on the medium where they travel, for example if it is a string, they depend on the string length, linear density and tension of the string. If it is a system of springs, it depends on the Hook’s constant k and the mass. If it is an acoustic wave in water or the air, it depends on density and compressibility parameter K. Because of acoustic waves we can hear, because of electromagnetic waves we can see. Elastic waves let us see through the earth’s interior with the use of seismology. Sometimes waves are characterized according to the direction of particle vibration. For example, if the particle vibration is aligned with the energy displacement direction, they are called compressional, while if they are orthogonal to that direction (or quasi–orthogonal in the case of anisotropic media) they are called transverse or shear. In 1D, the spring is a simple example of a generator of compressional waves, while a guitar string is a simple example of transverse or shear waves. This appendix show some simple derivations of these two types of waves. All wave equation derivations in elastic solids are the result of two basic laws from continuum mechancis: • Hook’s Law • Newton’s second law. In the case of a string only Newton’s second law is needed. In fluids the continuity equation and Newton’s second law are needed. So, in general, all elastic or acoustic wave equations are consequences of Newton second’s law. In the case of electromagnetic waves, they can be derived directly from Maxwell equations.

86

Derivations of wave equations

Figure B.1: System of springs and masses m1

B.1

m2

m3

m4

u(x, t)

u(x + h, t)

u(x + 2h, t)

One dimensional compressional wave equation

I present two derivations. The first derivation illustrates the case of a one–dimensional compressional wave in a system of springs. It is interesting because it shows the transicion between the discrete to the continuum. The second derivation is a more classical derivation based on continuum mechanic concepts.

B.1.1

From a system of springs to a bar

Figure B.1 shows a system of springs with masses. Let us for the moment focus on the mass m3 and study the interaction of this mass with its two neighbors m2 and m4 . For simplicity we assume that all springs have the same Hook’s constant k and mass m. The quantity u(x, t) represents the displacement of the mass at the location x from its center of equilibrium (where the force is zero). Then by Newton second law mass m3 is moving with some force FH = ma where a=

∂ 2 u(x + h, t) ∂t2

The force FH is the resultant of two forces by its neighbor masses, that is FH = k[u(x + 2h, t) − u(x + h, t)] − k[u(x + h, t) − u(x, t)] = k [u(x + 2h, t) − 2u(x + h, t) + u(x, t)] so m

∂ 2 u(x + h, t) = k [u(x + 2h, t) − 2u(x + h, t) + u(x, t)] ∂t2

or ∂ 2 u(x + h, t) k = [u(x + 2h, t) − 2u(x + h, t) + u(x, t)] 2 ∂t m

(B.1.1)

We are ready to pass to the limit, that sends us from a discrete system of springs into the contiuum. Assume we have an infinite number of masses between springs and that h → 0. That is, we are converting an infinite system of springs into a bar with density ρ and cross–section area A0 .

B.1 One dimensional compressional wave equation

87

We now use the relationship between the Hook’s constant and the Young’s modulus E, that is k=

EA0 h

with A0 the cross section of the string. We see that m = ρV = ρA0 h and so EA0 E k = 2 = 2 m ρA0 h ρh Now taking the limit in B.1.1 we find ∂ 2 u(x, h) E ∂u(x, t) = . ∂t2 ρ ∂x2

(B.1.2)

For this particular wave equation, the wave speed is given by s E v= . ρ

B.1.2

From continuum mechanics tools

Assume a homogenous bar with stress σ distributions along the x direction as shown in figure B.2. Here σ=

P A

is the axial stress and A is the cross–section area of the bar. As shown in the figure, ∂σ the stress in two close points is approximately given by σ at the point x and σ + ∆x ∂x at x + ∆x. Now, from Newton’s second law the total the total force resulting from the stresses on left and right is −σA + (σ +

∂σ ∂ 2u ∆x)A = mρ 2 . ∂x ∂t

Now since m = A∆x we find −σ + (σ +

∂σ ∂ 2u ∆x) = ρ 2 ∆x. ∂x ∂t

88

Derivations of wave equations

Figure B.2: A solid bar

P ∆x

σ

σ+

∂σ ∆x ∂x

or ∂σ ∂ 2u =ρ 2. ∂x ∂t

(B.1.3)

We now apply Hook’s law, that is σ = E, where E is the Young’s modulus and  =

∂u , so equation B.1.3 becomes ∂x

ρ ∂ 2u ∂ 2u = , ∂x2 E ∂t2 which agrees with equation B.1.2.

B.2

One dimensional transverse wave equation

This derivation is taken from Armstead and Karls article. Figure B.3 shows a piece of string with two tension vectors at the ends. The angles and bending of the string are exagerated to make the picture easy to visualize. We assume that the motion is pure transverse. That is, that the energy propagates along the horizontal while the vibration is taken place along the vertical. So, since there is not horizontal vibration the horizontal sum of forces is zero. For this to happen the horizontal tension should be constant. That is

T (x) cos θ = T (x + ∆x) cos(θ + ∆θ) = T.

(B.2.4)

B.2 One dimensional transverse wave equation

89

Figure B.3: A piece of string (brown). The piece of string represents the displacement function u(x, t). T (x + ∆x)

θ + ∆θ

θ T (x)

For the vertical direction we apply Newton’s second law. That is, −T (x) sin θ + T (x + ∆x) sin(θ + ∆θ) = m

∂ 2u (x, t) ∂t2

(B.2.5)

where u(x, t) measures the displacement in the vertical direction. Now from B.2.4 T (x) =

T cos θ

and T (x + ∆x) =

T cos(θ + ∆θ)

and plugging this into B.2.5 we find −T tan θ + T tan(θ + ∆θ) = µ∆x

∂ 2u (x, t). ∂t

(B.2.6)

Note that we use m = µ∆x, where µ is the linear density of the string. Since tan θ and tan(θ + ∆θ) are the tangents at each end of the string piece, they measure the slope of the displacement function u(x, t). That is tan θ =

∂u (x, t) and ∂x

tan(θ + ∆θ) =

∂u (x + ∆x, t); ∂x

We now sustitue this last equation in B.2.6, and find −T

∂u ∂u ∂ 2u (x, t) + T (x + ∆x, t) = ρ ∆x 2 (x, t), ∂x ∂x ∂t

from which, after dividing by T, ∆x is ∂u ∂u (x + ∆x, t) − (x, t) ρ ∂ 2u ∂x ∂x = (x, t). ∆x T ∂t2

90

Derivations of wave equations The final step is to take the limit ad ∆x → 0, so ∂ 2u ρ ∂ 2u . = ∂x2 T ∂ 2t

The wave speed for this particular wave is s v=

B.3

T . ρ

Elastic Waves in 3D Anisotropic Media

This derivation is taken from Aki and Richards [1] The idea of the derivation of the wave equation is again based on • Hook’s Law, and • Newton’s second law. The general point in the three–dimensional space is noted by x. The function ui (x, t) in the three–dimensional space and time is the particle displacement function along the i–direction, (i = 1, 2, 3, for xi coordinate axis) In general strain is a tensor (a matrix) defined as 1 eij = (ui,j + uj,i ) 2 where the shorthand notation ui,j means ui,j =

∂ui (x, t) . ∂xi

The total forces acting on a volume consist of • Body Forces: for example gravity and magnetic forces • Surface Forces: these are contact forces, for example shear and normal forces along the surface of a body. Figure B.3 illustrates a vector of a surface force acting on a piece of volume associated with a normal vector n. A force f acting on a point ξ and in an given specific time τ and along a given direction i can be described with the help of the Dirac and Kronercker delta symbols as follows fi (x, t) = Aδ(x − ξ)δ(t − τ )δin .

B.3 Elastic Waves in 3D Anisotropic Media

91

Figure B.4: A body and a surface force acting on it

T (n) n

Newton’s second law applied to the body, indicates that the total force on the body is equal to ma, where a is the acceleration of the body’s center of mass. If we consider the body as a continuum of particles and on each particle we assume that Newton’s second law is applied, the total forces are integrated by Z Z Z ∂ 2 ui f dV + T (s) dS . ρ dV = (B.3.7) ∂t2 | V {z } | S {z } | V {z } Mass times acceleration

Total body forces

Total surface forces

We want to study the Cauchy stress tensor, which from definition is ∆F ∆S→0 ∆S

T (n) = lim

where n denotes the normal to the surface element. An infinite number of traction vectors can act on a point with arbitrary directions. We want to prove the Cauchy’s law that shows how any arbitrary stress tensor can be written in term of some basis functions which are stress tensors along the main coordinate directions, but before we can prove the Cauchy’s law, we first prove the Cauchy’s lemma.

B.3.1

Cauchy’s Lemma

Cauchy’s lemma states that traction vectors acting on opposite sides of a surface are equal and opposite. This is T (n) = −T (−n).

(B.3.8)

This is equivalent to the third Newton’s law of action and reaction. Initially we want to take a volume slab around a particle of interest as shown in figure B.5, with an infinitesimal thickness δ. We treat each term on equation B.3.7 individually, and use the mean value theorem for integrals. The first integral, Z ∂ 2u ρ 2 dV = a∆m ∂t ∆V

92

Derivations of wave equations

Figure B.5: A small disc within a stressed medium. T (n)

n -n

T (−n)

where the last equality comes from the mean value theorem for integrals. Here a is the acceleration vector at some point inside the volume. The body force acting in the volume is, similarly Z f dV = f ∆V. ∆V

where f is a force vector inside the body. The total mass is Z ∆m = ρdV = ρ∆V ∆V

where ρ is the density at some point interior to the volume. The resultant force from the two tractions (last integral ) is T (n) ∆S + T (−n) ∆S, where ∆S is the area of the slab. So we find a∆m = f ∆V + T (n) ∆S + T (−n) ∆S. We divide by ∆S, a ρ ∆V ∆V =f + T (n) + T (−n), ∆S ∆S now, clearly ∆V = 0. δ→0 ∆S lim

since ∆V = ∆Sδ, so we find the Cauchy lemma T (n) = −T (−n).

B.3 Elastic Waves in 3D Anisotropic Media

93

x3 T (n) n τ11

τ12 τ22 τ32

τ21

τ13

τ23

τ31

x2

τ33

x1

B.3.2

Cauchy’s Law

The Cauchy’s law states that there exists a Cauchy stress tensor τ which maps the normal to a surface to the traction vector acting on that surface according to T = τ n, or in short hand index notation Ti = τij nj , where we assume the Einstein’s notation of sum over repeated indices. This is a very useful formula because each traction can be computed along arbitrary normal directions with the help of the Cauchy stress tensor τ . It is similar to what we do in calculus when we know the gradient and use it to compute directional derivatives, only that here we have one more dimension. The gradient is a tensor of order 1 and the Cauchy stress tensor is a tensor of order 2.

94

Derivations of wave equations

Bibliography [1] K. Aki and R. P. Quantitative Seismology: Theory and Methods, volume I, pages 156–158. W. H. Freeman and Company, 1980. [2] N. Bleistein. Mathematical methods for wave phenomena. Academic Press, Inc, 1984. [3] A. Komech and K. A. Principles of Partial Differential Equations. Springer-Verlag, 2012.