Solutions Partial Differential Equations 1a4s4e

Contents 1 Classification of Differential Equations 2 Models in one dimension 2.1 Heat Flow in a bar; Fourier’s law . . . . 2.2 The hanging bar . . . . . . . . . . . . . 2.3 The wave equation for a vibrating string 2.4 Advection; kinematic waves . . . . . . .

1

. . . .

. . . .

. . . .

. . . .

. . . .

. . . .

. . . .

. . . .

. . . .

. . . .

. . . .

. . . .

. . . .

. . . .

. . . .

. . . .

. . . .

. . . .

. . . .

. . . .

. . . .

. . . .

. . . .

. . . .

. . . .

3 3 4 5 5

3 Essential Linear Algebra 3.1 Linear systems as linear operator equations . . . . . 3.2 Existence and uniqueness of solutions to Ax = b . . 3.3 Basis and dimension . . . . . . . . . . . . . . . . . . 3.4 Orthogonal bases and projections . . . . . . . . . . . 3.5 Eigenvalues and eigenvectors of a symmetric matrix

. . . . .

. . . . .

. . . . .

. . . . .

. . . . .

. . . . .

. . . . .

. . . . .

. . . . .

. . . . .

. . . . .

. . . . .

. . . . .

. . . . .

. . . . .

. . . . .

. . . . .

. . . . .

. . . . .

. . . . .

. . . . .

. . . . .

. . . . .

. . . . .

7 7 9 10 11 12

4 Essential Ordinary Differential Equations 4.1 Background . . . . . . . . . . . . . . . . . . . 4.2 Solutions to some simple ODEs . . . . . . . . 4.3 Linear systems with constant coefficients . . . 4.4 Numerical methods for initial value problems 4.5 Stiff systems of ODEs . . . . . . . . . . . . .

. . . . .

. . . . .

. . . . .

. . . . .

. . . . .

. . . . .

. . . . .

. . . . .

. . . . .

. . . . .

. . . . .

. . . . .

. . . . .

. . . . .

. . . . .

. . . . .

. . . . .

. . . . .

. . . . .

. . . . .

. . . . .

. . . . .

. . . . .

. . . . .

15 15 17 19 21 23

5 Boundary Value Problems in Statics 5.1 The analogy between BVPs and linear algebraic systems 5.2 Introduction to the spectral method; eigenfunctions . . . 5.3 Solving the BVP using Fourier series . . . . . . . . . . . 5.4 Finite element methods for BVPs . . . . . . . . . . . . . 5.5 The Galerkin method . . . . . . . . . . . . . . . . . . . 5.6 Piecewise polynomials and the finite element method . .

. . . . . .

. . . . . .

. . . . . .

. . . . . .

. . . . . .

. . . . . .

. . . . . .

. . . . . .

. . . . . .

. . . . . .

. . . . . .

. . . . . .

. . . . . .

. . . . . .

. . . . . .

. . . . . .

. . . . . .

. . . . . .

. . . . . .

. . . . . .

. . . . . .

. . . . . .

25 25 27 29 32 33 35

6 Heat Flow and Diffusion 6.1 Fourier series methods for the heat equation . . . . . . . 6.2 Pure Neumann conditions and the Fourier cosine series . 6.3 Periodic boundary conditions and the full Fourier series 6.4 Finite element methods for the heat equation . . . . . . 6.5 Finite elements and Neumann conditions . . . . . . . . .

. . . . .

. . . . .

. . . . .

. . . . .

. . . . .

. . . . .

. . . . .

. . . . .

. . . . .

. . . . .

. . . . .

. . . . .

. . . . .

. . . . .

. . . . .

. . . . .

. . . . .

. . . . .

. . . . .

. . . . .

. . . . .

. . . . .

39 39 41 44 47 49

7 Waves 7.1 The homogeneous wave equation without boundaries 7.2 Fourier series methods for the wave equation . . . . 7.3 Finite element methods for the wave equation . . . . 7.4 Resonance . . . . . . . . . . . . . . . . . . . . . . . .

. . . .

. . . .

. . . .

. . . .

. . . .

. . . .

. . . .

. . . .

. . . .

. . . .

. . . .

. . . .

. . . .

. . . .

. . . .

. . . .

. . . .

. . . .

. . . .

. . . .

. . . .

. . . .

53 53 54 57 58

. . . .

. . . .

. . . .

. . . .

. . . . .

i

. . . .

. . . . .

. . . .

. . . . .

. . . . .

. . . .

. . . .

ii

CONTENTS 7.5

Finite difference methods for the wave equation . . . . . . . . . . . . . . . . . . . . . . . . . . .

61

8 First-order PDEs and the Method of Characteristics 8.1 The simplest PDE and the method of characteristics . . . . . . . . . . . . . . . . . . . . . . . . 8.2 First-order quasilinear PDEs . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . 8.3 Burgers’s equation . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . .

63 63 65 67

9 Green’s Functions 9.1 Green’s functions for BVPs in ODEs: Special cases . . . . 9.2 Green’s functions for BVPs in ODEs: the symmetric case 9.3 Green’s functions for BVPs in ODEs: the general case . . 9.4 Introduction to Green’s functions for IVPs . . . . . . . . . 9.5 Green’s functions for the heat equation . . . . . . . . . . . 9.6 Green’s functions for the wave equation . . . . . . . . . .

. . . . . .

. . . . . .

. . . . . .

. . . . . .

. . . . . .

. . . . . .

. . . . . .

. . . . . .

. . . . . .

. . . . . .

. . . . . .

. . . . . .

. . . . . .

. . . . . .

. . . . . .

. . . . . .

. . . . . .

. . . . . .

. . . . . .

. . . . . .

. . . . . .

71 71 72 75 76 79 80

10 Sturm-Liouville Eigenvalue Problems 10.2 Properties of the Sturm-Liouville operator . . . . . . . 10.3 Numerical Methods for Sturm-Liouville problems . . . 10.4 Examples of Sturm-Liouville problems . . . . . . . . . 10.5 Robin boundary conditions . . . . . . . . . . . . . . . 10.6 Finite element methods for Robin boundary conditions 10.7 The theory of Sturm-Liouville problems: an outline . .

. . . . . .

. . . . . .

. . . . . .

. . . . . .

. . . . . .

. . . . . .

. . . . . .

. . . . . .

. . . . . .

. . . . . .

. . . . . .

. . . . . .

. . . . . .

. . . . . .

. . . . . .

. . . . . .

. . . . . .

. . . . . .

. . . . . .

. . . . . .

. . . . . .

83 83 84 85 85 86 88

11 Problems in Multiple Spatial Dimensions 11.1 Physical models in two or three spatial dimensions . . . . . . 11.2 Fourier series on a rectangular domain . . . . . . . . . . . . . 11.3 Fourier series on a disk . . . . . . . . . . . . . . . . . . . . . . 11.4 Finite elements in two dimensions . . . . . . . . . . . . . . . . 11.5 The free-space Green’s function for the Laplacian . . . . . . . 11.6 The Green’s function for the Laplacian on a bounded domain 11.7 Green’s function for the wave equation . . . . . . . . . . . . . 11.8 Green’s functions for the heat equation . . . . . . . . . . . . .

. . . . . . . .

. . . . . . . .

. . . . . . . .

. . . . . . . .

. . . . . . . .

. . . . . . . .

. . . . . . . .

. . . . . . . .

. . . . . . . .

. . . . . . . .

. . . . . . . .

. . . . . . . .

. . . . . . . .

. . . . . . . .

. . . . . . . .

. . . . . . . .

. . . . . . . .

. . . . . . . .

. . . . . . . .

89 89 90 94 96 100 102 107 109

12 More about Fourier Series 12.1 The complex Fourier series . . . . . . . . . . . 12.2 Fourier series and the FFT . . . . . . . . . . . 12.3 Relationship of sine and cosine series to the full 12.4 Pointwise convergence of Fourier series . . . . . 12.5 Uniform convergence of Fourier series . . . . . 12.6 Mean-square convergence of Fourier series . . . 12.7 A note about general eigenvalue problems . . .

. . . . . . .

. . . . . . .

. . . . . . .

. . . . . . .

. . . . . . .

. . . . . . .

. . . . . . .

. . . . . . .

. . . . . . .

. . . . . . .

. . . . . . .

. . . . . . .

. . . . . . .

. . . . . . .

. . . . . . .

. . . . . . .

. . . . . . .

. . . . . . .

. . . . . . .

111 111 113 116 116 117 117 119

. . . . . . . . . . methods . . . . .

. . . .

. . . .

. . . .

. . . .

. . . .

. . . .

. . . .

. . . .

. . . .

. . . .

. . . .

. . . .

. . . .

. . . .

. . . .

. . . .

. . . .

. . . .

121 121 122 124 124

. . . . . .

. . . . . .

. . . . . . . . . . . . . . . . Fourier series . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . .

13 More about Finite Element Methods 13.1 Implementation of finite element methods . . . . . . . 13.2 Solving sparse linear systems . . . . . . . . . . . . . . 13.3 An outline of the convergence theory for finite element 13.4 Finite element methods for eigenvalue problems . . . .

Chapter 1

Classification of Differential Equations 1.

(a) First-order, homogeneous, linear, nonconstant-coefficient, scalar ODE. (b) Second-order, inhomogeneous, linear, constant-coefficient, scalar PDE. (c) First-order, nonlinear, scalar PDE.

3.

(a) Second-order, nonlinear, scalar ODE. (b) First-order, inhomogeneous, linear, constant-coefficient, scalar PDE. (c) Second-order, inhomogeneous, linear, nonconstant-coefficient, scalar PDE.

5.

(a) No. (b) Yes. (c) No.

7. There is only one such f : f (t) = t cos (t). 9. For any constant C 6= 1, the function w(t) = Cu(t) is a different solution of the differential equation.

1

Chapter 2

Models in one dimension 2.1

Heat Flow in a bar; Fourier’s law

1. The units of κ are energy per length per time per temperature, for example, J/(cm s K) = W/(cm K). 3. The units of Aρcu∆x are g J Kcm = J. cm2 cm3 g K r 5. ∂u (`, t) = Aκ , t > t0 . ∂x 7. Measuring x in centimeters, the steady-state temperature distribution is u(x) = 0.1x + 20, the solution of

−

d2 u = 0, u(0) = 20, u(100) = 30. dx2

The heat flux is

du (x) = −0.802 · 0.1 = −0.0802 W/cm2 . dx . Since the area of the bar’s cross-section is π cm2 , the rate at which energy is flowing through the bar is 0.0802π = 0.252 W (in the negative direction). 9. Consider the right endpoint (x = `). Let the temperature of the bath be u` , so that the difference in temperature between the bath and the end of the bar is u(`, t) − u` . The heat flux at x = ` is, according to Fourier’s law, −κ

∂u (`, t), ∂x so the statement that the heat flux is proportional to the temperature flow is written −κ

∂u (`, t) = α (u(`, t) − u` ) , ∂x where α > 0 is a constant of proportionality. This can be simplified to −κ

αu(`, t) + κ

∂u (`, t) = αu` . ∂x

The condition at the right end is similar: ∂u (0, t) = αu0 . ∂x (The sign change appears because, at the left end, a positive heat flux means that heat flows into the bar, while at the right end the opposite is true.) 11. The IBVP is ∂2u ∂u − D 2 = f (x, t), 0 < x < `, t > t0 , ∂t ∂x u(x, t0 ) = ψ(x), 0 < x < `, ∂u (0, t) = 0, t > t0 , ∂x ∂u (`, t) = 0, t > t0 . ∂x αu(0, t) − κ

3

4

CHAPTER 2. MODELS IN ONE DIMENSION 13.

(a) The rate is u − u 0 ` . ` (b) The rate of mass transfer varies inversely with ` and directly with the square of r. −πr2 D

2.2

The hanging bar

1. The sum of the forces on the cross-section originally at x = ` must be zero. As derived in the text, the internal (elastic) force acting on this cross-section is du −Ak(`) (`), dx while the external traction results in a total force of pA (force per unit area times area). Therefore, we obtain −Ak(`)

du (`) + pA = 0, dx

or k(`) 3.

du (`) = p. dx

The other boundary condition, of course, is u(0) = 0. (a) The stiffness is the ratio of the internal restoring force to the relative change in length. Therefore, if a bar of length ` and cross-sectional area A is stretched to a length of 1 + p (0 < p < 1) times `, then the bar will pull back with a force of 195pA gigaNewtons. Equivalently, this is the force that must be applied to the end of the bar to stretch the bar to a length of (1 + p)`. (b) 196/195 m, or approximately 1.005 m. (c) The BVP is

d2 u = 0, 0 < x < 1, dx2 u(0) = 0, 9 du (1) = 109 . 195 · 10 dx It is easy to show by direct integration that the solution is u(x) = x/195, and therefore u(1) = 1/195. Since u is the displacement (in meters), the final position of the end of the bar is 1 + 1/195 m, that is, the stretched bar is 196/195 m. 5. The wave equation is 2 ∂2u 11 ∂ u 7.9 · 103 − 1.95 · 10 = f (x, t), 0 < x < 1. ∂t2 ∂x2 3 The units of f are N/m (force per unit volume). The units of the first term on the left are −195 · 109

kgm/s2 kg m N = = 3, 3 2 m s m3 m and the units of the second term on the left are N 1 N = 3. m2 m m 7.

(a) We have ∂2u (x, t) = −c2 θ2 u(x, t), ∂t2 while

∂2u (x, t) = −θ2 u(x, t). ∂x2

Therefore, 2 ∂2u 2∂ u − c = −c2 θ2 u + c2 θ2 u = 0. ∂t2 ∂x2 (b) Regardless of the value of θ, u(0, t) = 0 holds for all t. The only way that u(`, t) = 0 can hold for all t is if

sin (θ`) = 0, so θ must be one of the values θ=

nπ , n = 0, ±1, ±2, . . . . `

2.3. THE WAVE EQUATION FOR A VIBRATING STRING

2.3

5

The wave equation for a vibrating string

1. Units of acceleration (length per time squared). 3. Units of velocity (length per time). 5. The internal force acting on the end of the string at, say, x = `, is T

∂u (`, t). ∂x

(2.1)

If this end can move freely in the vertical direction, force balance implies that (2.1) must be zero. 7. If u(x, t) = f (x − ct), then ∂v (x, t)=f 0 (x − ct), ∂x ∂v (x, t)= − cf 0 (x − ct), ∂t

∂2v (x, t)=f 00 (x − ct), ∂x2 ∂2v (x, t)=c2 f 00 (x − ct), ∂t2

Therefore, ∂2u ∂2u (x, t) − c2 2 (x, t) = c2 f 00 (x − ct) − c2 f 00 (x − ct) = 0, 2 ∂t ∂x as desired. Similarly, if v(x, t) = f (x + ct), then ∂v (x, t)=f 0 (x + ct), ∂x ∂v (x, t)=cf 0 (x + ct), ∂t and hence

2.4

∂2v (x, t)=f 00 (x + ct), ∂x2 ∂2v (x, t)=c2 f 00 (x + ct), ∂t2

2 ∂2v 2∂ v (x, t) − c (x, t) = c2 f 00 (x + ct) − c2 f 00 (x + ct) = 0. ∂t2 ∂x2

Advection; kinematic waves

1. The solution u of ∂u ∂u +c = 0, 0 < x < `, t > 0, ∂t ∂x u(x, 0) = u0 (x), 0 < x < `, u(0, t) = φ(t), t > 0. can be found by reasoning similar to that used to solve the pure IVP (see (2.23) and following in the text). The cross-section at point x at time t was ct units to the left at time t = 0. If x − ct > 0, then the value of u(x, t) is given by the initial condition: u(x, t) = u0 (x − ct). If x − ct < 0, then u(x, t) is determined by the boundary condition u(0, t) = φ(t), and we must ask: At what time was the cross-section, now at point x at time t, found at x = 0? The answer is t0 , where c(t − t0 ) = x, that is, t0 = t − (1/c)x. Thus u(x, t) = φ(t − (1/c)x) if x − ct < 0. It is easy to by a direct calculuation that u0 (x − ct), x − ct > 0, u(x, t) = φ t − 1c x , x − ct < 0 solves the given IBVP. 3. Suppose u solves ∂u ∂u +c = 0, −∞ < x < ∞, t > 0, ∂t ∂x u(x, 0) = φ(x), −∞ < x < ∞, where c > 0, and suppose that φ(x) = 0 for all x ≤ a. The solution is u(x, t) = φ(x − ct), so u(x, t) = 0 if x − ct ≤ a, that is, if t > (x − a)/c.

Chapter 3

Essential Linear Algebra 3.1

Linear systems as linear operator equations

1. A function of the form f (x) = ax + b is linear if and only if b = 0. Indeed, if f : R → R is linear, then f (x) = ax, where a = f (1). 3. Let f : R2 → R2 be defined by f (x) =

x21 + x22 x2 − x21



If we take x = (1, 1), then f (x) = (2, 0), while 2x = (2, 2) and hence f (2x) = (8, −2) 6= 2f (x). Thus f is not linear. 5.

(a) Vector space (subspace of C[0, 1]). (b) Not a vector space; does not contain the zero function. (Subset of C[0, 1], but not a subspace.) (c) Vector space (subspace of C[0, 1]). (d) Vector space. (e) Not a vector space; does not contain the zero polynomial. (Subset of Pn , but not a subspace.)

7. Suppose u ∈ C 1 [a, b] is any nonzero function. Then L(2u) = 2

du du + (2u)3 = 2 + 8u3 , dx dx

while 2Lu = 2

du + 2u3 . dx

Since L(2u) 6= 2Lu, L is not linear. 9. Define K : C 2 [a, b] → C[a, b] by Ku = x2

du d2 u − 2x + 3u. dx2 dx

Then K is linear: d2 d (αu) − 2x (αu) + 3(αu) dx2 dx du d2 u = αx2 2 − 2αx + 3αu dx dx d2 u du = α x2 2 − 2x + 3u = αKu, dx dx

K(αu) = x2

d2 d (u + v) − 2x (u + v) + 3(u + v) dx2 dx 2 d u d2 v du dv = x2 + − 2x + + 3(u + v) dx2 dx2 dx dx d2 u du d2 v dv = x2 2 − 2x + 3u + x2 2 − 2x + 3v = Ku + Kv. dx dx dx dx

K(u + v) = x2

7

8

CHAPTER 3. ESSENTIAL LINEAR ALGEBRA Notice how we used the linearity of the derivative operator. 11. Let ρ, c, and κ be constants, and define L : C 2 (R2 ) → C(R2 ) by Lu = ρc

∂u ∂2u −κ 2. ∂t ∂x

We prove that L is linear as follows: L(αu) = ρc

L(u + v) = ρc

13.

∂2 ∂u ∂2u ∂u ∂2u ∂ (αu) − κ 2 (αu) = αρc − ακ 2 = α ρc − κ 2 = αLu, ∂t ∂x ∂t ∂x ∂t ∂x

2 ∂2 ∂u ∂v ∂ u ∂ ∂2v (u + v) − κ 2 (u + v) = ρc + −κ + ∂t ∂x ∂t ∂t ∂x2 ∂x2 2 ∂u ∂ u ∂2v ∂v = ρc − κ 2 + ρc − κ 2 = Lu + Lv. ∂t ∂x ∂t ∂x

(a) The proof is a direct calculation:

A(αx + βz)

= = = = = =

x1 z1 α +β x2 z2 a11 a12 αx1 + βz1 a21 a22 αx2 + βz2 a11 (αx1 + βz1 ) + a12 (αx2 + βz2 ) a21 (αx1 + βz1 ) + a22 (αx2 + βz2 ) α(a11 x1 + a12 x2 ) + β(a11 z1 + a12 z2 ) α(a21 x1 + a22 x2 ) + β(a21 z1 + a22 z2 ) a11 x1 + a12 x2 a11 z1 + a12 z2 α +β a21 x1 + a22 x2 a21 z1 + a22 z2 a11 a21

a12 a22

αAx + βAz.

This shows that the operator defined by A is linear. (b) We have (A(αx + βy))i =

n X

aij (αxj + βyj )

j=1

and (αAx + βAy)i = α

X

aij xj + β

j=1

n X

aij yj .

j=1

Now, (A(αx + βy))i

=

n X

aij (αxj + βyj )

j=1

=

n X

(αaij xj + βaij yj )

j=1

=

n X

αaij xj +

j=1

=

α

n X j=1

=

n X

βaij yj

j=1

aij xj + β

n X

aij yj

j=1

(αAx + βAy)i ,

as desired. The third equality follows from the commutative and associative properties of addition, while the fourth equality follows from the distributive property of multiplication over addition.

3.2. EXISTENCE AND UNIQUENESS OF SOLUTIONS TO AX = B

3.2

9

Existence and uniqueness of solutions to Ax = b

1. The range of A is {α(1, −1) : α ∈ R} , which is a line in the plane. See Figure 3.1. 3

2

x2

1

0

−1

−2

−3 −3

−2

−1

0 x

1

2

3

1

Figure 3.1: The range of A (Exercise 3.2.1). 3.

(a) A is nonsingular. (b) A is nonsingular. (c) A is singular. Every vector in the range is of the form 1 1 x1 Ax = 1 1 x2 x1 + x2 = . x1 + x2 That is, every y ∈ R2 whose first and second components are equal lies in the range of A. Thus, for example, Ax = b is solvable for 1 b= , 1 but not for 1 b= . 2

5. The solution set is not a subspace, since it cannot contain the zero vector (A0 = 0 6= b). 7. The null space of L is the set of all first-degree polynomials: N (L) = {u : [a, b] → R : u(x) = mx + c for some m, c ∈ R} . 9.

(a) The only functions in C 2 [a, b] that satisfy d2 u =0 dx2 are functions of the form u(x) = C1 x + C2 . The Neumann boundary condition at the left endpoint implies that C1 = 0, and the Dirichlet boundary condition at the right endpoint implies that C2 = 0. Therefore, only the zero function is a solution to Lm ˜ u = 0, and so the null space of Lm ˜ is trivial. 2 (b) Suppose f ∈ C[a, b] is given and u ∈ Cm [a, b] satisfies ˜ −

d2 u (x) = f (x), a < x < b. dx2 Integrating once yields, by the fundamental theorem of calculus, Z x 2 Z x d u − (s) ds = f (s) ds 2 a dx a Z x du du ⇒ − (x) + (a) = f (s) ds dx dx a Z x du ⇒ (x) = − f (s) ds. dx a −

10

CHAPTER 3. ESSENTIAL LINEAR ALGEBRA The last step follows from the Neumann condition at x = a. We now integrate again: Z x du f (s) ds (x) = − dx a Z b Z bZ z du ⇒ (z) dz = − f (s) ds dz x dx x a Z bZ z ⇒ u(b) − u(x) = − f (s) ds dz Z aZ z x a f (s) ds dz. ⇒ u(x) = x

a

This shows that Lm ˜ u = f has a unique solution for each f ∈ C[a, b]. 11. By the fundamental theorem of calculus, x

Z u(x) =

f (s) ds a

satisfies Du = f . However, this solution is not unique; for any constant C, Z x u(x) = f (s) ds + C a

is another solution.

3.3 1.

Basis and dimension (a) Both equal (14, 4, −5). P (b) Both equal n j=1 (vj )i xj .

3. The given set is not a basis. If A is the matrix whose columns are the given vectors, then N (A) is not trivial. 5. There is a typo in the problem statement; the given polynomials are linearly independent: c1 (1 − x + 2x2 ) + c2 (1 − 2x2 ) + c3 (1 − 3x + 7x2 ) = 0 ⇔(c1 + c2 + c3 ) + (−c1 − 3c3 )x + (2c1 − 2c2 + 7c3 )x2 = 0 ⇔

c1 + c2 + c3 −c1 − 3c3 2c1 − 2c2 + 7c3

= 0, = 0, = 0.

A direct calculation (using Gaussian elimination) shows that this last system of equations has only the trivial solution, and hence the three polynomials form a linearly independent set. 7. We know that P2 has dimension 3, and therefore it suffices to show either that the given set of three vectors spans P2 or that it is linearly independent. If p ∈ P2 , then we write c1 = p(x1 ), c3 = p(x3 ), c3 = p(x3 ) and define the polynomial q(x) = c1 L1 (x) + c2 L2 (x) + c3 L3 (x). Then q(x1 ) = c1 L1 (x1 ) + c2 L2 (x1 ) + c3 L3 (x1 ) = c1 · 1 + c2 · 0 + c3 · 0 = c1 = p(x1 ), and, similarly, q(x2 ) = p(x2 ), q(x3 ) = p(x3 ). But then p and q are second-degree polynomials that agree at three distinct points, and three points determine a quadratic (just as two points determine a line). Therefore, p(x) = c1 L1 (x) + c2 L2 (x) + c3 L3 (x), and we have shown that every p ∈ P2 can be written as a linear combination of L1 , L2 , L3 . This completes the proof.

3.4. ORTHOGONAL BASES AND PROJECTIONS

11

2 9. Let L : CN [a, b] → C[a, b] be the second derivative operator. We wish to find a basis for the null space of L. A 2 function u ∈ CN [a, b] belongs to the null space of L if and only if it satisfies

du d2 u du (0) = (`) = 0. (x) = 0 for all x ∈ [a, b], dx2 dx dx The only functions satisfying the differential equation are first-degree polynomials, u(x) = mx + c. Moreover, the boundary conditions are satisfied if and only if the slope m is zero. Thus u belongs to the null space of L if and only if u(x) = c for some c ∈ R, that is, if and only if u is a constant function. A basis is therefore the set {u1 }, where u1 (x) = 1 for all x ∈ [a, b].

3.4 1.

Orthogonal bases and projections (a) The verification is a direct calculation of vi · vj for the 6 combinations of i, j. For example, v1 · v1

= =

v1 · v2

= =

1 1 1 1 1 1 √ √ +√ √ +√ √ 3 3 3 3 3 3 1 1 1 + + = 1, 3 3 3 1 1 1 1 √ √ + √ 0+ √ 3 2 3 3 1 1 √ − √ = 0, 6 6



1 −√ 2



and so forth. (b) x

= =

(v1 · x)v1 + (v2 · x)v2 + (v3 · x)v3 4 2 √ v1 + 0v2 + − √ v3 . 3 6

3. We have kx + yk2

=

(x + y, x + y)

=

(x, x + y) + (y, x + y)

=

(x, x) + (x, y) + (y, x) + (y, y)

=

(x, x) + 2(x, y) + (y, y)

=

kxk2 + 2(x, y) + kyk2 .

Therefore, kx + yk2 = kxk2 + kyk2 if and only if (x, y) = 0. 5. First of all, if (y, z) = 0 for all z ∈ W, then since wi ∈ W for each i, we have, in particular, (y, wi ) = 0, i = 1, 2, . . . , n. Suppose, on the other hand, that (y, wi ) = 0, i = 1, 2, . . . , n. If z is any vector in W , then, since {w1 , w2 , . . . , wn } is a basis for W , there exist scalars α1 , α2 , . . . , αn such that z=

n X i=1

αi wi .

12

CHAPTER 3. ESSENTIAL LINEAR ALGEBRA We then have (y, z)

=

y,

n X

! αi wi

i=1 n X

=

i=1 n X

=

αi (y, wi ) αi · 0

i=1

=

0.

This completes the proof. 7. The best approximation to the given data is y = 2.0109x − 0.0015151. See Figure 3.2. 6 5.5 5 4.5 y

4 3.5 3 2.5 2 1

1.5

2 x

2.5

3

Figure 3.2: The data from Exercise 3.4.7 and the best linear approximation. 9. The projection of g onto P2 is

or

√ 2 5(π 2 − 12) 2 q1 (x) + 0q2 (x) + q3 (x) π π3 10(π 2 − 12) 2 π 2 − 12 π 2 − 12 2 + − 60 x + 60 x . 3 3 π π π π3

See Figure 3.3. 1

0.5

0

y=sin(π x) quadratic approx. −0.5 0

0.2

0.4

0.6

0.8

1

Figure 3.3: The function g(x) = sin (πx) and its best quadratic approximation over the interval [0, 1] (see Exercise 3.4.9).

3.5

Eigenvalues and eigenvectors of a symmetric matrix

1. The eigenvalues of A are λ1 = 200 and λ2 = 100, and the corresponding (normalized) eigenvectors are −0.8 0.6 u1 = , u2 = , 0.6 0.8

3.5. EIGENVALUES AND EIGENVECTORS OF A SYMMETRIC MATRIX

13

respectively. The solution is x=

b · u2 b · u1 u1 + u2 = λ1 λ2



1 1



3. The eigenvalues and eigenvectors of A are λ1 = 2, λ2 = 1, λ3 = −1 and    1   1 √

 u1 =  

3 √1 3 √1 3

√

   , u2 =   

2

   0   , u3 = 

− √12

.

√1 6 − √26 √1 6

  . 

The solution of the Ax = b is   −3/2 b · u1 b · u2 b · u3 x= u1 + u2 + u3 =  1/2  . λ1 λ2 λ3 5/2 5. The computation of ui · b/λi costs 2n operations, so computing all n of these ratios costs 2n2 operations. Computing the linear combination is then equivalent to another n dot products, so the total cost is 2n2 + n ∗ (2n − 1) = 4n2 − n = O(4n2 ). 7.

(a) For k = 2, 3, . . . , n − 1, we have 1 (− sin ((k − 1)jπh) + 2 sin (kjπh) − sin ((k + 1)jπh)) Ls(j) = h2 k 1 (− sin (kjπh) cos (jπh) + cos (kjπh) sin (jπh) + 2 sin (kjπh) = h2 − sin (kjπh) cos (jπh) − cos (kjπh) sin (jπh)) 1 = (2 − 2 cos (jπh)) sin (kjπh) h2 1 (j) = (2 − 2 cos (jπh)) sk . h2 Thus

1 (j) Ls(j) = 2 (2 − 2 cos (jπh)) sk , k = 2, 3, . . . , n − 1. h k

Similar calculations show that this formula holds for k = 1 and k = n also. Therefore, s(j) is an eigenvector of L with eigenvalue 2 − 2 cos (jπh) . λj = h2 (b) The eigenvalues λj are all positive and are increasing with the frequency j. (c) The right-hand side b can be expressed as b = s(1) · b s(1) + s(2) · b s(2) + · · · + s(n) · b s(n) , while the solution x of Lx = b is x=

s(1) · x (1) s(2) · x (2) s(n) · x (n) s + s + ··· + s . λ1 λ2 λn

Since λj increases with j, this shows that, in producing x, the higher frequency components of b are dampened more than are the lower frequency components of b. Thus x is smoother than b.

Chapter 4

Essential Ordinary Differential Equations 4.1

Background

1. Define x1 = u, x2 =

du . dt

Then dx1 dt dx2 dt 3. Define x1 = u, x2 =

=

x2 ,

=

1 − x1 . 2

du d2 u d3 u , x3 = 2 , x4 = 3 . dt dt dt

Then dx1 dt dx2 dt dx3 dt dx4 dt

=

x2 ,

=

x3 ,

=

x4 ,

=

2x3 − x1 + sin t.

5. Define

dp dq , x3 = q, x4 = . dt dt

x1 = p, x2 = Then the first-order system is dx1 dt dx2 m1 dt dx3 dt dx4 m2 dt

=

x2 ,

=

f1 (x1 , x3 , x2 , x4 ),

=

x4 ,

=

f2 (x1 , x3 , x2 , x4 ).

7. Let u1 (t) = e−t , u2 (t) = e−2t .

15

16

CHAPTER 4. ESSENTIAL ORDINARY DIFFERENTIAL EQUATIONS (a) We have du1 d 2 u1 (t) = −e−t , (t) = e−t , dt dt2 and therefore du1 d 2 u1 (t) + 3 (t) + 2u(t) = e−t − 3e−t + 2e−t = 0. dt2 dt Similarly, du2 d 2 u2 (t) = −2e−2t , (t) = 4e−2t , dt dt2 and thus d2 u2 du2 (t) + 3 (t) + 2u(t) = 4e−2t − 6e−2t + 2e−2t = 0. dt2 dt (b) The Wronskian of u1 , u2 is u2 (t) e−t = du2 (t) −e−t dt

u1 (t) W (t) = du1 (t) dt

e−2t = −e−3t . −2e−2t

Since W (t) 6= 0 for all t, we see that {u1 , u2 } is linearly independent. 9.

(a) Let u1 , u2 be two functions defined on an interval, let t0 be a point in that interval, and suppose u1 (t0 ) du 1 (t0 ) dt

u2 (t0 ) 6 0. = du2 (t0 ) dt

We wish to prove that {u1 , u2 } is linearly independent. We argue by contradiction, and suppose {u1 , u2 } is linearly dependent. Then there exist c1 , c2 ∈ R such that c1 u1 +c2 u2 = 0. This means that c1 u1 (t)+c2 u2 (t) = 1 2 0 for all t in the given interval, which in turn implies that c1 du (t) + c2 du (t) = 0 for all t. But then, taking dt dt t = t0 , we obtain c1 u1 (t0 ) + c2 u2 (t0 ) = 0, du2 du1 (t0 ) + c2 (t0 ) = 0. c1 dt dt In matrix-vector form, this system is

u1 (t0 ) du1 (t0 ) dt

u2 (t0 ) du2 (t0 ) dt



c1 c2



=

0 0

.

Since (c1 , c2 ) 6= (0, 0), this implies that the matrix

u1 (t0 ) du1 (t0 ) dt

u2 (t0 ) du2 (t0 ) dt



is singular, contradicting the assumption that its determinant is nonzero. This contradiction completes the proof. (b) Let u1 (t) = cos (t), u2 (t) = cos (2t). Then u1 (t) du 1 (t) dt

u2 (t) du2 (t) dt

cos (t) = − sin (t)

cos (2t) = cos (2t) sin (t) − 2 cos (t) sin (2t). −2 sin (2t)

We see that W (π/2) = −1 6= 0, and hence, by part (a), {u1 , u2 } is linearly independent. Nevertheless, W (0) = 0, which shows that the fact that {u1 , u2 } is linearly independent does not imply that W (t) 6= 0 for all t, unless u1 and u2 are both solutions to the same second-order linear differential equation.

4.2. SOLUTIONS TO SOME SIMPLE ODES

4.2 1.

17

Solutions to some simple ODEs (a) We first note that the zero function is a solution of (4.9), so S is nonempty. If u, v ∈ S and α, β ∈ R, then w = αu + βv satisfies a

d2 w dw + cw +b dt2 dt

= = = =

d2 d [αu + βv] + b [αu + βv] + c(αu + βv) dt2 dt 2 dv d u d2 v du +β a α 2 +β 2 +b α + c(αu + βv) dt dt dt dt 2 2 d u du d v dv α a 2 +b + cu + β a 2 + b + cv dt dt dt dt α · 0 + β · 0 = 0. a

Thus w is also a solution of (4.9), which shows that S is a subspace. (b) As explained in the text, no matter what the values of a, b, c (a 6= 0), the solution space is spanned by two functions. Thus S is finite-dimensional and it has dimension it is easy to see that each at most 2. Moreover, of the sets {er1 t , er2 t } (r1 6= r2 ), eµt cos (λt), eµt sin (λt) , and ert , tert is linearly independent, since a set of two functions is linearly dependent if and only if one of the functions is a multiple of the other. Thus, in every case, S has a basis with two functions, and so S is two-dimensional. 3. Suppose that the characteristic polynomial of (4.9) has a single, repeated root r = −b/(2a) (so b2 − 4ac = 0). Then, as shown in the text, u(t) = c1 ert + c2 tert is a solution of (4.9) for each choice of c1 , c2 . We wish to show that, given k1 , k2 ∈ R, there is a unique choice of c1 , c2 such that u(0) = k1 , du/dt(0) = k2 . We have du (t) = rc1 ert + c2 ert + rc2 tert . dt Therefore, the equations u(0) = k1 , du/dt(0) = k2 simplify to 1 0 c = k. r 1 Since the coefficient matrix is obviously nonsingular (regardless of the value of r), there is a unique solution c1 , c2 for each k1 , k2 . 5.

(a) The only solution is u(x) = 0. (b) The only solution is u(x) = 0. (c) The function u(x) = sin (πx) is a nonzero solution. It is not unique, since any multiple of it is another solution.

7. By the product rule and the fundamental theorem of calculus, Z t Z t d d a(t−s) at −as e f (s) ds = e e f (s) ds = dt t0 dt t0

aeat

t

e−as f (s) ds + eat e−at f (t)

t0

Z =

Z

t

ea(t−s) f (s) ds + f (t).

a t0

Therefore, with u(t) = u0 ea(t−t0 ) +

Z

t

ea(t−s) f (s) ds, t0

we have du (t) dt

=

au0 ea(t−t0 ) + a

Z

t

ea(t−s) f (s) ds + f (t) t0

=

au(t) + f (t).

Thus u satisfies the differential equation. We also have Z t0 u(t0 ) = u0 ea(t0 −t0 ) + ea(t−s) f (s) ds = u0 + 0 = u0 , t0

so the initial condition holds as well.

18

CHAPTER 4. ESSENTIAL ORDINARY DIFFERENTIAL EQUATIONS 9. u(t) =

1 2

11. u(t) =

1 2

13.

sin2 (t) e−t + cos (t) + sin (t)



2

(a) The solution is x(t) = et

/2

.

(b) The solution is x(t) = (1 + cos (1) − cos (t))/t3 . 2

(c) The solution is y(t) = 1/2 + (1/2)e−t . 15. Suppose we substitute u(t) = c1 (t)u1 (t) + c2 (t)u2 (t) into a

du d2 u +b + cu = f (t), dt2 dt

where u1 , u2 are solutions of the homogeneous version of this ODE, and assume that dc1 dc2 u1 + u2 = 0. dt dt We then have du1 du2 du = c1 + c2 , dt dt dt d2 u d 2 u1 d 2 u2 dc2 du2 dc1 du1 + . = c1 2 + c2 2 + 2 dt dt dt dt dt dt dt Substituting into the ODE, we obtain d 2 u2 du1 du2 d 2 u1 dc1 du1 dc2 du2 + + b c1 + c2 + c(c1 u1 + c2 u2 ) = f (t) a c1 2 + c2 2 + dt dt dt dt dt dt dt dt 2 2 dc1 du1 d u1 d u2 dc2 du2 du1 du2 ⇒c1 a 2 + b + cu1 + c2 a 2 + b + cu2 + a + = f (t) dt dt dt dt dt dt dt dt dc2 du2 dc1 du1 + = f (t) ⇒a dt dt dt dt dc2 du2 dc1 du1 + = a−1 f (t), ⇒ dt dt dt dt as desired. Notice that we have used the fact that u1 , u2 solve the homogenous ODE: a

d 2 u1 du1 d 2 u2 du2 +b + cu1 = 0, a 2 + b + cu2 = 0. 2 dt dt dt dt

17. Let us define v(s) = u(es ); then u(t) = v(ln (t)). It follows that 1 d2 v 1 dv du 1 dv d2 u (t) = (ln (t)), (t) = 2 2 (ln (t)) − 2 (ln (t)). 2 dt t ds dt t ds t ds It follows that du d2 u (t) + at (t) + bu(t) = 0 dt2 dt d2 v dv dv ⇔ 2 (ln (t)) − (ln (t)) + a (ln (t)) + bv(ln (t)) = 0 ds ds ds d2 v dv ⇔ 2 (ln (t)) + (a − 1) (ln (t)) + bv(ln (t)) = 0 ds ds d2 v dv ⇔ 2 (s) + (a − 1) (s) + bv(s) = 0. ds ds t2

Thus the change of variables t = es transforms the Euler equation into a constant coefficient ODE.

4.3. LINEAR SYSTEMS WITH CONSTANT COEFFICIENTS

4.3

19

Linear systems with constant coefficients

1. The solution is x(t)

1 1 √ u1 − √ e−t u2 + 3 2  1 1 −t − e + 16 e−2t 3 2  1 − 13 e−2t  3

=

=

1 3

1 √ e−2t u3 6   .

+ 12 e−t + 16 e−2t

3. The general solution is x(t) =

c1 et + c2 e−t c1 et − c2 e−t

.

5. x0 must be a multiple of the vector (1, −1). 7.

9.

  55 − 3e−2t − 45e−t + 20t − 7 cos (t) − 11 sin (t) 1   10 + 6e−2t + 20t − 16 cos (t) − 8 sin (t) x(t) = 60 −2t −t −5 − 3e + 45e + 20t − 37 cos (t) + 19 sin (t) (a) The solution is 3 1 3 −t 1 3t e + e , y(t) = e−t − e3t . 2 2 2 2 The population of the first species (x(t)) increases exponentially, while the population of the second species (y(t)) goes to zero in finite time (y(t) = 0 at t = ln (3)/4). Thus the second species becomes extinct, while the first species increases without bound. x(t) =

(b) If the initial populations are x(0) = r, y(0) = s, then the solution to the IVP is x(t) =

1 1 (r + s)e−t + (r − s)e3t , y(t) = (r + s)e−t + (s − r)e3t . 2 2

Therefore, if r = s, both populations decay to zero exponentially; that is, both species die out. 11. Assume v(t; s) satisfies d2 v dv (0; s) = 0 + θ2 v = 0, v(0; s) = 0, dt2 dt for all s, and define t

Z

v(t − s; s) ds.

u(t) = 0

Then du (t) = v(0; t) + dt

dv d2 u (t) = (0; t) + dt2 dt (using the fact that v(0; t) = 0 and

t

Z 0

dv (t − s; s) ds = dt t

Z 0

t

Z

dv (t − s; s) ds, dt Z t 2 d2 v d v (t − s; s) ds = f (t) + (t − s; s) ds 2 dt2 0 dt 0

dv (0; t) dt

= 0 for all t). We then have Z t 2 Z t d2 u d v (t) + θ2 u(t) = f (t) + (t − s; s) ds + θ2 v(t − s; s) ds 2 2 dt 0 dt 0 Z t 2 d v 2 = f (t) + (t − s; s) + θ v(t − s; s) ds dt2 0 = f (t)

since

d2 v (t − s; s) + θ2 v(t − s; s) = 0 for all s. dt2

Also, Z 0 du dv (t) = (t − s; s) ds = 0, dt 0 0 dt and thus u satisfies both the differential equations and the initial conditions. Z

u(t) =

0

v(t − s; s) ds = 0,

20

CHAPTER 4. ESSENTIAL ORDINARY DIFFERENTIAL EQUATIONS

13. Suppose that, for all s, v(t; s) satisfies the ODE dk u dk−1 u du + ak−1 (t) k−1 + · · · + a1 (t) + a0 (t)u = 0 k dt dt dt and also the initial conditions dk−2 u dk−1 u du (s) = 0, . . . , (s) = 0, (s) = f (s). dt dtk−2 dtk−1

u(s) = 0, Define

t

Z

v(t; s) ds.

u(t) = 0

Then

Z t Z t dv dv du (t) = v(t; t) + (t; s) ds = (t; s) ds dt dt 0 0 dt (since v(t; t) = 0 for all t). Similarly, Z t 2 Z t 2 d v d2 u dv d v (t; t) + (t) = (t; s) ds = (t; s) ds, 2 2 dt2 dt dt 0 0 dt .. . Z t k−1 Z t k−1 dk−1 u dk−2 v d v d v (t) = (t; t) + (t; s) ds = (t; s) ds k−1 k−1 dtk−1 dtk−2 dt dt 0 0 and finally dk u dk−1 v (t) = (t; t) + dtk dtk−1

t

Z 0

dk v (t; s) ds = f (t) + dtk

Z 0

t

dk v (t; s) ds. dtk

We then see that dk u dk−1 u du + a0 (t)u + ak−1 (t) k−1 + · · · + a1 (t) dtk dt dt Z t k d v dk−1 v dv =f (t) + (t; s) + a (t) (t; s) + · · · + a (t) (t; s) + a (t)v(t; s) ds 1 0 k−1 dtk dtk−1 dt 0 =f (t), since

dk−1 v dv dk v (t; s) + ak−1 (t) k−1 (t; s) + · · · + a1 (t) (t; s) + a0 (t)v(t; s) = 0 k dt dt dt

by assumption. Also, dj u (0) = dtj

0

Z 0

dj v (t; s) ds = 0, j = 0, 1, . . . , k − 1. dtj

15. The solution to d2 u du +3 + 2u = 0, dt2 dt u(0) = 0, du (0) = f (s) dt is given by v(t; s) = (e−t − e−2t )f (s). Therefore, by Duhamel’s principle, the solution of d2 u du +3 + 2u = f (t), dt2 dt u(0) = 0, du (0) = 0 dt is

t

Z

Z v(t − s; s) ds =

u(t) = 0

0

t

e−(t−s) − e−2(t−s) f (s) ds.

4.4. NUMERICAL METHODS FOR INITIAL VALUE PROBLEMS

4.4 1.

21

Numerical methods for initial value problems (a) Four steps of Euler’s method yield an estimate of 0.71969. (b) Two steps of the improved Euler method yield an estimate of 0.80687. (c) One step of the RK4 method yields an estimate of 0.82380. Euler’s method gives no correct digits, the improved Euler method gives one correct digit, and RK4 gives three correct digits. Each of the methods evaluated f (t, u) four times.

3.

(a) Let u1 = x, u2 = dx/dt, u3 = y, u4 = dy/dt. Then the system is du1 dt du2 dt du3 dt du4 dt

=

u2 ,

=

u1 + 2u4 −

=

u4 ,

=

−2u2 + u3 −

µ2 (u1 + µ1 ) µ1 (u1 − µ2 ) − , r1 (u1 , u3 )3 r2 (u1 , u3 )3

µ2 u 3 µ1 u 3 − . r1 (u1 , u3 )3 r2 (u1 , u3 )3

(b) The routine ode45 from MATLAB (version 5.3) required 421 steps to produce a graph with the ending point apparently coinciding with the initial value. The graph of y(t) versus x(t) is given in Figure 4.1. The graphs of x(t) and y(t), together with the times steps, are given in Figure 4.2. They show, not surprisingly, that the time step is forced to be small precisely when the coordinates are changing rapidly. 1 0.8 0.6 0.4 0.2

y

0 −0.2 −0.4 −0.6 −0.8 −1 −1.5

−1

−0.5

0 x

0.5

1

1.5

Figure 4.1: The orbit of the satellite in Exercise 4.4.3.

x(t)

2 0

y(t)

−2 0 1

2

3

4

5

6

7

4

5

6

7

3 4 Step index

5

6

7

t

0

−1 0 0.15

Time step

1

1

2

3 t

0.1 0.05 0 0

1

2

Figure 4.2: The coordinates of the satellite in Exercise 4.4.3 (top two graphs) with the step lengths taken (bottom graph). (c) The minimum step size taken by ode45 was 3.28 · 10−4 , and using this step length over the entire interval of [0, T ] would require almost 19000 steps. This is to be compared to the 421 steps taken by the adaptive algorithm. (Note: The exact results for minimum step size, etc., will differ according to the algorithm and tolerances used.)

22

CHAPTER 4. ESSENTIAL ORDINARY DIFFERENTIAL EQUATIONS 5.

(a) We first note that a + (t1 − t0 ) ≤ t ≤ b + (t1 − t0 ) if and only if a ≤ t − (t1 − t0 ) ≤ b, so the function v is well-defined. In fact, {v(t) : t ∈ [a + (t1 − t0 ), b + (t1 − t0 )]} =

{u(t − (t1 − t0 )) : t ∈ [a + (t1 − t0 ), b + (t1 − t0 )]}

=

{u(t) : t ∈ [a, b]}

(just replace t − (t1 − t0 ) by t in the last step). We have dv (t) dt

=

d [u(t − (t1 − t0 ))] dt du (t − (t1 − t0 )) dt f (u(t − (t1 − t0 )))

=

f (v(t)),

= =

so v satisfies the ODE. Finally, v(t1 ) = u(t1 − (t1 − t0 )) = u(t0 ) = u0 . Therefore v is a solution to the given IVP. (b) Consider the IVP du dt u(0)

=

t,

=

1,

=

t,

=

1

which has solution u(t) = t2 /2 + 1. The IVP dv dt v(1) has solution v(t) = t2 /2 + 1/2, which is not equal to u(t − (1 − 0)) = u(t − 1) =

1 (t − 1)2 + 1. 2

7. The exact solution is x(t) =

2 − 2e−t , 1 t e, 2

0 < t < ln 2, t > ln 2.

(a) The following errors were obtained at t = 0.5: ∆t 1/4 1/8 1/16 1/32 1/64 1/128

Error 2.4332e-05 1.3697e-06 8.1251e-08 4.9475e-09 3.0522e-10 1.8952e-11

By inspection, we see that as ∆t is cut in half, the error is reduced by a factor of approximately 16, as expected for O(∆t4 ) convergence. (b) The following errors were obtained at t = 2.0:

4.5. STIFF SYSTEMS OF ODES

23 ∆t 1/4 1/8 1/16 1/32 1/64 1/128

Error 5.0774e-03 3.2345e-03 3.1896e-04 1.0063e-04 8.9026e-06 3.4808e-06

The error definitely does not exhibit O(∆t4 ) convergence to zero. The reason is the lack of smoothness of the function 1 + |x − 1|; the rate of convergence given in the text only applies to ODEs defined by smooth functions. When integrating from t = 0 to t = 0.5, the nonsmoothness of the right-hand side is not encountered since x(t) < 1 on this interval. This explains why we observed good convergence in the first part of this problem.

4.5 1.

Stiff systems of ODEs (a) The exact solution is x(t) =

1 2



e−t + et e−t − et

.

(b) The norm of the error is approximately 0.089103. (c) The norm of the error is approximately 0.10658. 3. The largest value is ∆t = 0.04. 5. (a) Applying the methods of Section 4.3, we find the solution 1 3e−t − e−100t x(t) = . 2 3e−t + e−100t (b) By trial and error, we find that ∆t ≤ 0.02 is necessary for stability. Specifically, integrating from t = 0 to t = 1, n = 49 (i.e. ∆t = 1/49) yields kx49 k > kx0 k , while n ≥ 50 yields kxn k ≤ kx0 k . (i)

(c) With xi+1 = xi + ∆tAxi and y1 = u1 · xi , we have u1 · xi+1 = u1 · xi + ∆tu1 · Axi (i+1)

= y1 + ∆t(Au1 ) · xi

(i+1)

= y1 + ∆tλ1 u1 · xi

(i+1)

= y1 + ∆tλ1 y1

⇒

y1

⇒

y1

⇒

y1

⇒

(i+1) y1

(i)

(i)

(i)

= (1 +

(i)

(i) ∆tλ1 )y1 .

A similar calculation shows that (i+1)

y2

(i)

= (1 + ∆tλ2 )y2 .

For stability, we need |1 + ∆tλ1 |

≤

1,

|1 + ∆tλ2 |

≤

1.

Since the eigenvalues of A are −1, −100, it is not hard to see that the upper bound for ∆t is ∆t ≤ 0.02, just as was determined by experiment. (d) For the backward Euler method, a similar calculation shows that y1

(i+1)

=

(1 − ∆tλ1 )−1 y1 ,

(i+1) y2

=

(1 − ∆tλ2 )−1 y2 .

(i)

Stability is guaranteed for any ∆t, since |1 − ∆tλi |−1 < 1 for every ∆t > 0.

(i)

Chapter 5

Boundary Value Problems in Statics 5.1 1.

The analogy between BVPs and linear algebraic systems (a) Since MD is linear, it suffices to show that the null space of MD is trivial. If MD u = 0, then, du/dx = 0, that is, u is a constant function: u(x) = c. But then the condition u(0) = 0 implies that c = 0, and so u is the zero function. Thus N (MD ) = {0}. 1 [0, `] and MD u = f , then we have (b) If u ∈ CD

du (x) = f (x), 0 < x < ` and u(0) = 0, dx which imply that x

Z u(x) =

f (s) ds. 0

But we also must have u(`) = 0, which implies that `

Z

f (s) ds = 0. 0

If f ∈ C[0, `] does not satisfy this condition, then it is impossible for MD u = f to have a solution. 3.

(a) If v, w ∈ S are both solutions of Lu = f , then Lv = Lw, or (by linearity) L(v − w) = 0. Therefore v − w ∈ N (L). But, since S is a subspace, v − w is also in S. If the only function in both N (L) and S is the zero function, then v − w must be the zero function, that is, v and w must be the same function. Therefore, if N (L) ∩ S = {0}, then Lu = f can have at most one solution for any f . (b) We have already seen that Lu = f has a solution for any f ∈ C[0, `] (see the discussion immediately preceding Example 5.1). We will use the result from the first part of this exercise to show that the solution is unique. The null space of L is the space of all first degree polynomials: N (L) = {u : [0, `] → R : u(x) = ax + b for some a, b ∈ R} . i. Suppose u ∈ N (L) ∩ S. Then u(x) = ax + b for some a, b ∈ R and du dx

` = 0 ⇒ a = 0. 2

Then u(x) = b, and Z

`

u(x) dx = 0 ⇒ b` = 0 ⇒ b = 0. 0

Therefore u is the zero function, and so N (L) ∩ S = {0}. The uniqueness property then follows from the first part of this exercise.

25

26

CHAPTER 5. BOUNDARY VALUE PROBLEMS IN STATICS ii. Suppose u ∈ N (L) ∩ S. Then u(x) = ax + b for some a, b ∈ R and u(0) = 0 ⇒ b = 0. Then u(x) = ax, and du (`) = 0 ⇒ a = 0. dx Therefore u is the zero function, and so N (L) ∩ S = {0}. The uniqueness property then follows from the first part of this exercise. 5. The condition

`

Z 0



du (x) dx

2 dx = 0

(5.1)

implies that du/dx is zero, and hence that u is constant. The boundary conditions on u would then imply that u is the zero function. But, by assumption, u is nonzero (we assumed that (u, u) = 1). Therefore, (5.1) cannot hold. 7.

(a) If Lm ˜ u = 0, then u has the form u(x) = ax + b. The first boundary condition, du/dx(0) = 0, implies that a = 0, and then the second boundary condition, u(`) = 0, yields b = 0. Therefore, u is the zero function and N (Lm ˜ ) is trivial. (b) For any f ∈ C[0, `], the function Z

`

z

Z

u(x) =

f (s) ds dz x

0

2 belongs to Cm ˜ [0, `] and satisfies

d2 u (x) = f (x), 0 < x < `. dx2 This shows that Lm ˜ u = f has a solution for any f ∈ C[0, `], and therefore R(Lm ˜ ) = C[0, `]. −

2 (c) Suppose u, v ∈ Cm ˜ [0, `]. Then

= = = = = =

`

d2 u (x)v(x) dx 2 0 dx ` Z ` du dv du − (x)v(x) + (x) (x) dx dx dx dx 0 0 Z ` du dv (x) (x) dx (since v(`) = 0, du/dx(0) = 0) dx 0 dx ` Z ` d2 v dv u(x) 2 (x) dx u(x) (x) − dx dx 0 0 Z ` 2 d v − u(x) 2 (x) dx (since u(`) = 0, dv/dx(0) = 0) dx 0 (u, Lm ˜ v). Z

(Lm ˜ u, v)

−

Thus Lm ˜ is symmetric. (d) Suppose λ is an eigenvalue of Lm ˜ with corresponding eigenfunction u, and assume that u has been normalized so that (u, u) = 1. Then λ = λ(u, u) = (λu, u)

= = = = >

(Lm ˜ u, u) Z ` 2 d u − (x)u(x) dx 2 0 dx ` Z ` 2 du du − (x)u(x) + (x) dx dx dx 0 0 Z ` 2 du (x) dx (since u(`) = 0, du/dx(0) = 0) dx 0 0.

2 The last step follows because every nonzero function in Cm ˜ [0, `] has a nonzero derivative.

5.2. INTRODUCTION TO THE SPECTRAL METHOD; EIGENFUNCTIONS

27

9. Define u(x) = x(1 − x), v(x) = x2 (1 − x). Then a direct calculation shows that (M u, v) =

7 , 30

(u, M v) =

4 . 15

but

11.

2 (a) Suppose u, v ∈ CR [0, `]. Then `

d2 u (x)v(x) dx 2 0 dx ` Z ` du dv du (x) (x) dx − κ (x)v(x) + κ dx dx 0 dx 0 ` ` Z ` du dv d2 v − κ (x)v(x) + κu(x) (x) − κ u(x) 2 (x) dx dx dx dx 0 0 0 du du dv −κ (`)v(`) + κ (0)v(0) + u(`)κ (`) dx dx dx dv −u(0)κ (0) + (u, LR v). dx Z

(LR u, v)

−κ

= = = =

Applying the boundary conditions on u, v, we obtain du du dv dv (`)v(`) + κ (0)v(0) + u(`)κ (`) − u(0)κ (0) dx dx dx dx αu(`)v(`) + αu(0)v(0) − αu(`)v(`) − αu(0)v(0) = 0, −κ

=

so (LR u, v) = (u, LR v), as desired. (b) If LR u = 0, then u must be a first degree polynomial: u(x) = ax + b. The boundary conditions imply that a and b must satisfy the following equations: −κa + αb

=

0,

(α` + κ)a + αb

=

0.

The determinant is computed as follows: −κ α` + κ

α = −α(2κ + α`) < 0. α

The only solution is a = b = 0, that is, u = 0, so N (LR ) is trivial.

5.2

Introduction to the spectral method; eigenfunctions

1. If n 6= m, then `

Z (un , um )

=

nπx

sin

mπx

dx ` ` (n + m)πx (n − m)πx − cos dx cos ` ` 0 ` (n − m)πx (n + m)πx ` ` sin − sin . (n − m)π ` (n + m)π ` 0 sin

0

=

1 2

=

1 2

Z

`

This last expression simplifies to four , each including the sine function evaluated at an integer multiple of π, and hence each equal to zero. Thus (un , um ) = 0 for n 6= m. 3. The eigenpairs are (2n − 1)2 π 2 , cos 4`2



(2n − 1)πx , n = 1, 2, 3, . . . . 2`

28

CHAPTER 5. BOUNDARY VALUE PROBLEMS IN STATICS 5. The characteristic polynomial is r2 − r + (λ − 5), so the characteristic roots are √ √ 1 + 21 − 4λ 1 − 21 − 4λ , r2 = . r1 = 2 2

(5.2)

Case 1: λ = 21/4. In this case, the characteristic roots are r = 1/2, 1/2 and the general solution of the ODE is u(x) = c1 ex/2 + c2 xex/2 . The boundary condition u(0) = 0 yields c1 = 0, and then the boundary condition u(1) = 1 implies that c2 = 0. Thus there is no nonzero solution to the BVP, and λ = 21/4 is not an eigenvalue. Case 2: λ < 21/4. In this case, the characteristic roots, given by (5.2), are real and distinct. The general solution of the ODE is u(x) = c1 er1 x + c2 er2 x . The boundary conditions lead to the system c1 + c2

=

0,

c1 r1 er1 + c2 r2 er2

=

0.

This system has the unique solution c1 = c2 = 0, so there is no nonzero solution for λ < 21/4. Hence no λ < 21/4 is an eigenvalue. Case 3: λ > 21/4. In this case, the roots are complex conjugate: √ 1 4λ − 21 1 r1 = − θi, r2 = + θi, θ = . 2 2 2 The general solution of the ODE is u(x) = (c1 cos (θx) + c2 sin (θx)) ex/2 . The boundary condition u(0) = 0 yields c1 = 0, so any eigenfunction is of the form u(x) = sin (θx)ex/2 . The Neumann condition at x = 1 is equivalent to the equation tan (θ) = −2θ, θ > 0. Although this equation cannot be solved explicitly, a simple graph shows that there are infinitely many solutions 0 < θ1 < θ2 < · · · , with 2k − 1 2k + 1 θk ∈ π, π , k = 1, 2, . . . . 2 2 Define λk by √ 4λk − 21 θk = , 2 that is, 21 λk = θk2 + , k = 1, 2, . . . . 4 Then λk , k = 1, 2, . . ., are eigenvalues, and the corresponding eigenfunctions are vk (x) = sin (θk x)ex/2 . Using Newton’s method, we find that which yields

. . θ1 = 1.8366, θ2 = 4.8158, . . λ1 = 8.6231, λ2 = 28.4423.

The eigenfunctions are v1 , v2 , as given by the above formula. A direct calculation shows that . (v1 , v2 ) = −0.2092, so v1 and v2 are not orthogonal.

5.3. SOLVING THE BVP USING FOURIER SERIES 7.

(a)

P∞

2(−1)n+1 nπ

sin (nπx)

(b)

P∞

4 sin nπ 2 n2 π 2

sin (nπx)

(c)

P∞

n=1

(

n=1

)

n+1

n=1

12(−1) n3 π 3

29

sin (nπx)

720(−1)n+1 n=1 n5 π 5

P∞

(d) sin (nπx) The errors in approximating the original functions using 10 of the Fourier sine series are graphed in Figure 5.1. 1

0.03 0.02

0.5 0.01 0 0 −0.5 0

0.5 x

1

−0.01 0

−3

1.5

0.5 x

1

0.5 x

1

−5

x 10

5

x 10

1 0.5 0 0 −0.5 −1 0

0.5 x

−5 0

1

Figure 5.1: Errors in approximating four functions by 10 of their Fourier sine series (see Exercise 5.2.7). Top left: g(x) = x; Top right: h(x) = 21 − x − 21 ; Bottom left: m(x) = x − x3 ; Bottom right: k(x) = 7x − 10x3 + 3x5 . 9. sin (3πx) (That is, all of the Fourier sine coefficients are zero, except the third, which is one.) 11. The series have the form

∞ X

an cos

n=1

(2n − 1)πx , 2

where n

(2n−1)π) (a) an = − 4(2+(−1) ; π 2 (2n−1)2

(b) an =

32 cos ((2n−1)π/4) sin2 ((2n−1)π/8) ; π 2 (2n−1)2

(c) an = −

8(24+12(−1)n (2n−1)π+(2n−1)2 π 2 ) π 4 (2n−1)4

;

n

(d) an = −

8(5760+2880(−1) (2n−1)π+240(2n−1)2 π 2 +7(2n−1)4 π 4 ) π 6 (2n−1)6

.

The errors in approximating the original functions using 10 of the Fourier quarter-wave cosine series are graphed in Figure 5.2.

5.3 1.

Solving the BVP using Fourier series (a)

P∞

n=1

2(1−(−1)n ) n3 π 3

sin (nπx) n

P∞

3.

) (b) x + n=1 2(1−(−1) sin (nπx) n3 π 3 P∞ 32(−1)n+1 (2n−1)πx (a) n=1 (2n−1)4 π 4 sin 2

(b) 1 + (c)

P∞

16((−1)n π+2(−1)n+1 nπ−2)

n=1

2(1−(−1)n ) n=1 nπ (2+n2 π 2 )

P∞

(2n−1)4 π 4

sin (nπx)

cos



(2n−1)πx 2



30

CHAPTER 5. BOUNDARY VALUE PROBLEMS IN STATICS 1

0.04 0.02

0.5 0 0 −0.02 −0.5 0

0.5 x

1

0.01

−0.04 0

0.5 x

1

0.5 x

1

0.05

0

0

−0.01

−0.05

−0.02

−0.1

−0.03 0

0.5 x

1

−0.15 0

Figure 5.2: Errors in approximating four functions by 10 of their Fourier quarter-wave cosine series (see Exercise 5.2.11). Top left: g(x) = x; Top right: h(x) = 12 − x − 21 ; Bottom left: m(x) = x − x3 ; Bottom right: k(x) = 7x − 10x3 + 3x5 .

(d) x +

2(−1)n n=1 nπ (1+n2 π 2 )

P∞

sin (nπx)

5. We have d2 u (x) = −x, 0 < x < 1, dx2 so du (x) dx

= = =

Z x 2 du d u (0) + (s) dx 2 dx 0 dx Z x du (0) − s dx dx 0 x2 du (0) − . dx 2

Write C for du/dx(0); then another integration yields x

Z

u(x)

=

du (s) ds dx Z x s2 0+ C− ds 2 0 u(0) +

0

= =

Cx −

x3 . 6

Finally, we choose C so that u(1) = 0, which yields C = 1/6. This yields u(x) =

1 x − x3 . 6

7. The Fourier quarter-wave sine coefficients of −T d2 u/dx2 are bn = −

2T `

Z

2 `

Z

`

0

d2 u (x) sin dx2



(2n − 1)πx 2`

dx,

while those of u are an =

`

u(x) sin

0

(2n − 1)πx 2`

dx.

5.3. SOLVING THE BVP USING FOURIER SERIES

31

Using integration by parts twice, almost exactly as in (5.19), we can express bn in of an :

=

=

=

=

=

`

(2n − 1)πx d2 u (x) sin dx 2 2` 0 dx ( ` (2n − 1)πx du 2T (x) sin − ` dx 2` x=0 Z ` (2n − 1)πx (2n − 1)π du (x) cos − dx 2` 2` 0 dx Z ` 2T (2n − 1)π (2n − 1)πx du (x) cos dx 2`2 dx 2` 0 (since sin (0) = du/dx(`) = 0) ( ` (2n − 1)πx 2T (2n − 1)π u(x) cos 2`2 2` x=0 Z ` (2n − 1)π (2n − 1)πx + u(x) sin dx 2` 2` 0 Z T (2n − 1)2 π 2 2 ` (2n − 1)πx dx u(x) sin 4`2 ` 0 2` (since u(0) = cos ((2n − 1)π/2) = 0)

2T − `

Z

T (2n − 1)2 π 2 an . 4`2

This gives the desired result. (Actually, since it is known that the negative second derivative operator is symmetric under the mixed boundary conditions, we can just appeal to (5.25), which is the above calculation written abstractly.) 9. Let a1 , a2 , a3 , . . . be the Fourier quarter-wave sine coefficients of u; then −κ

∞ X κ(2n − 1)2 π 2 (2n − 1)πx d2 u (x) = a sin , n dx2 2002 200 n=1

where κ = 3/2. We also have 0.001 =

∞ X n=1

0.004 sin (2n − 1)π



(2n − 1)πx . 200

Setting the two series equal and solving for an , we find u(x) =

∞ X n=1

3.2 · 102 sin 3(2n − 1)3 π 3



(2n − 1)πx . 200

The temperature distribution is graphed in Figure 5.3. 3.5 3

temperature

2.5 2 1.5 1 0.5 0 0

20

40

60

80

100

x

Figure 5.3: The temperature distribution (in degrees Celsius) in Exercise 5.3.9.

32

CHAPTER 5. BOUNDARY VALUE PROBLEMS IN STATICS

5.4

Finite element methods for BVPs

2 1. Suppose f, g ∈ CD [0, `], so that

f (0) = f (`) = g(0) = g(`) = 0. Then −

= = = = =

`

k(x)

df dx



,g

df k(x) (x) g(x) dx dx 0 ` Z ` df df dg − k(x) (x)g(x) + k(x) (x) (x) dx (integration by parts) dx dx dx 0 0 Z ` df dg k(x) (x) (x) dx (using g(0) = g(`) = 0) dx dx 0 ` Z ` d dg dg f (x) k(x) (x) dx (integration by parts) k(x)f (x) (x) − dx dx dx 0 0 Z ` d dg − f (x) k(x) (x) dx (using f (0) = f (`) = 0) dx dx 0 d dg f, − k(x) . dx dx Z

=

d dx

−

d dx

3. Take p(x) = (x − c)3 (d − x)3 and v[c,d] as suggested in the hint. 2 [0, `] be a test function. We have 5. Let v ∈ CD d du − k(x) (x) + p(x)u(x) = f (x), 0 < x < ` dx dx d du ⇒ − k(x) (x) v(x) + p(x)u(x)v(x) = f (x)v(x), 0 < x < ` dx dx Z ` Z ` Z ` du d ⇒ − k(x) (x) v(x) dx + p(x)u(x)v(x) dx = f (x)v(x) dx dx 0 dx 0 0 ` Z ` Z ` Z ` du dv du ⇒ −k(x) (x)v(x) + k(x) (x) (x) dx + p(x)u(x)v(x) dx = f (x)v(x) dx dx dx dx 0 0 0 0 Z ` Z ` Z ` dv du p(x)u(x)v(x) dx = f (x)v(x) dx ⇒ k(x) (x) (x) dx + dx dx 0 0 0 Z ` Z ` du dv ⇒ k(x) (x) (x) + p(x)u(x)v(x) dx = f (x)v(x) dx. dx dx 0 0

The last step follows from the fact that v(0) = v(`) = 0. Thus the weak form of the given BVP is Z ` Z ` du dv 2 2 find u ∈ CD [a, b] such that k(x) (x) (x) + p(x)u(x)v(x) dx = f (x)v(x) dx for all v ∈ CD [a, b]. dx dx 0 0 7. The calculation is the same as in Exercise 5; we integrate by parts only in the first integral, and obtain the 2 following weak form: Find u ∈ CD [a, b] such that Z ` Z ` du dv du 2 k(x) (x) (x) + c(x) (x)v(x) + p(x)u(x)v(x) dx = f (x)v(x) dx for all v ∈ CD [a, b]. dx dx dx 0 0 Notice that now the left-hand side is not symmetric in u and v. 9. Repeating the calculation beginning on page 167, we obtain " Z # 2 2 Z ∂ 1 ` ∂u 1 ` ∂u Aρ(x) dx + Ak(x) dx ∂t 2 0 ∂t 2 0 ∂x Z ` 2 ∂u = − c dx, t > 0. ∂t 0

5.5. THE GALERKIN METHOD

33

This shows that derivative with respect to time of the total energy is always negative, and hence that the total energy is always decreasing.

5.5 1.

The Galerkin method (a) Yes. (b) No, a(f, f ) = 0 for f (x) = 1, but f 6= 0. (c) No, a(f, f ) = 0 for f (x) = 1, but f 6= 0.

3. Let FN = span {sin (πx), sin (2πx), . . . , sin (N πx)} . We will apply the Galerkin method to the BVP d2 u = f (x), 0 < x < `, dx2 u(0) = 0,

−k

u(`) = 0, using FN as the approximating subspace. The weak form of the BVP is Z ` Z ` dv du 2 2 find u ∈ CD [a, b] such that f (x)v(x) dx for all v ∈ CD [a, b], k (x) (x) dx = dx dx 0 0 and the Galerkin method takes the form find vN =

N X

Z

`

uj sin (jπx) such that

k 0

j=1

dvn dv (x) (x) dx = dx dx

Z

`

f (x)v(x) dx for all v ∈ FN . 0

We find the coefficients uj , j = 1, 2, . . . , N , by solving Ku = f , where Z ` Z ` dφi dφj (x) (x) dx, Fi = f (x)φi (x) dx Kij = k dx dx 0 0 and φj (x) = sin (jπx). We have `

Z 0

dφi dφj k (x) (x) dx = dx dx

(

`

Z

2

kijπ cos (jπx) cos (iπx) dx = 0

kj 2 π 2 , 2

0,

i = j, i 6= j.

This last step follows from the fact that the functions cos (πx), cos (2πx), . . . , cos (N πx) are mutually orthogonal with respect to the L2 (0, 1) inner product. Also, notice that Z ` Z ` cj f (x)φj (x) dx = f (x) sin (jπx) dx = , 2 0 0 where cj is the Fourier sine coefficient of f . Since K is diagonal, it is easy to solve Ku = f : uj =

cj /2 cj = 2 2. kj 2 π 2 /2 kj π

Thus vN (x) =

N X j=1

cj sin (jπx). kj 2 π 2

This is exactly the solution obtained by the method of Fourier series in Section 5.3.2. 5.

(a) The set S is the span of {x, x2 } and therefore is a subspace. (b) As discussed in the text, a(·, ·) automatically satisfies two of the properties of an inner product. Every function p in S satisfies p(0) = 0, so a(·, ·) also satisfies the third property (a(u, u) = 0 implies that u = 0) for exactly the same reason as given in the text for the subspace V (see page 174). (c) The best approximation is p(x) = (9 − 3e)x2 + (4e − 10)x.

34

CHAPTER 5. BOUNDARY VALUE PROBLEMS IN STATICS (d) The bilinear form is not an inner product on P2 , since a(1, 1) = 0 but 1 6= 0 (a constant is “invisible” to the energy inner product and norm). Therefore, every polynomial of the form (9 − 3e)x2 + (4e − 10)x + C ∈ S is equidistant from f (x) = ex in the energy norm. 7. Write p(x) = x(1 − x) and q(x) = x(1/2 − x)(1 − x). The approximation will be v(x) = u1 p(x) + u2 q(x), where Ku = f . The stiffness matrix K is " R1 K

= " =

dp dp (1 + x) dx (x) dx (x) dx

0 R1 (1 0 1 2 1 − 30

dq dp (x) dx (x) dx + x) dx # 1 − 30 , 3

R1

dq dp (1 + x) dx (x) dx (x) dx

0 R1 (1 0

#

dq dq + x) dx (x) dx (x) dx

40

and the load vector f is R1 " 1 # xp(x) dx 12 0 f = R1 = . 1 xq(x) dx − 120 0 Therefore, . u=





0.16412 −0.038168

.

The exact solution is 1 1 1 u(x) = − x2 + x − ln (1 + x). 4 2 4 ln 2 The exact and approximate solutions are graphed in Figure 5.4. 0.045 0.04 0.035 0.03 0.025 0.02 0.015 exact approximate

0.01 0.005 0 0

0.2

0.4

0.6

0.8

1

Figure 5.4: The exact and approximate solutions from Exercise 5.5.7.

9.

(a) It will take about 103 times as long, or 1000 seconds (almost 17 minutes), to solve a 1000 × 1000 system, and 1003 times as long, or 106 seconds (about 11.5 days) to solve a 10000 × 10000 system. (b) Gaussian elimination consists of a forward phase, in which the diagonal entries are used to eliminate nonzero entries below and in the same column, and a backward phase, in which the diagonal entries are used to eliminate nonzero entries above and in the same column. During the forward phase, at a typical step, there is only 1 nonzero entry below the diagonal, and only 5 arithmetic operations are required to eliminate it (1 division to compute the multiplier, and 2 multiplications and 2 additions to add a multiple of the current row to the next). Thus the forward phase requires O(5n) operations. A typical step of the backward phase requires 3 arithmetic operations (a multiplication and an addition to adjust the right-hand side and a division to solve for the unknown). Thus the backward phase requires O(3n) operations. The grand total is O(8n) operations. (c) It will take about 10 times as long, or 0.1 seconds, to solve a 1000 × 1000 tridiagonal system, and 100 times as long, or 1 second, to solve a 10000 × 10000 tridiagonal system.

5.6. PIECEWISE POLYNOMIALS AND THE FINITE ELEMENT METHOD

5.6

35

Piecewise polynomials and the finite element method

1. We wish to compute the piecewise linear finite element approximation to the solution of d2 u = x, 0 < x < 1, dx2 u(0) = 0,

−

u(`) = 0 using a uniform mesh with four elements. Such a mesh corresponds to h = 1/4, and there are three basis functions. The stiffness matrix is   8 −4 0 8 −4  K =  −4 0 −4 8 (as computed in Example 5.18 in the text). The load vector F ∈ R3 is defined by 1

Z Fi =

f (x)φi (x) dx. 0

A direct calculation shows that F1 = 1/16, F2 = 1/8, F3 = 3/16. Solving Ku = f yields  5   u=

128 1 16 7 128

 .

The exact solution to the BVP is u(x) = (x − x3 )/6. Figure 5.5 shows the exact solution and the piecewise linear approximation. 0.07 0.06 0.05 0.04 0.03 0.02 0.01 0 0

0.2

0.4

0.6

0.8

1

Figure 5.5: The exact (dashed curve) and approximate (solid curve) solutions from Exercise 5.6.1. 3.

(b) Here are the errors for n = 10, 20, 40, 80: n 10 20 40 80

maximum error 1.6586 · 10−3 4.3226 · 10−4 1.1036 · 10−4 2.7881 · 10−5

We see that, when n is doubled, the error decreases by a factor of approximately four. Thus 1 error = O . n2 5. The exact solution is u(x) = x(1 − x3 )/12. Here are the errors for n = 10, 20, 40, 80:

36

CHAPTER 5. BOUNDARY VALUE PROBLEMS IN STATICS n 10 20 40 80

maximum error 1.1286 · 10−3 2.9710 · 10−4 7.6186 · 10−5 1.9288 · 10−5

We see that, when n is doubled, the error descreases by a factor of approximately four. Thus 1 error = O . n2 7. The weak form is `

Z find u ∈ V such that 0

dv du k(x) (x) (x) + p(x)u(x)v(x) dx dx dx

`

Z

f (x)v(x) dx for all v ∈ V.

= 0

The bilinear form is Z

`

a(u, v) = 0

dv du k(x) (x) (x) + p(x)u(x)v(x) dx. dx dx

9. We wish to apply the piecewise linear finite element method to the BVP −

d2 u 1 + 2u = − x, 0 < x < 1, dx2 2 u(0) = 0, u(1) = 0.

We use a sequence of increasingly fine uniform meshes. The weak form of the BVP is Z 1 Z 1 Z 1 dv du 1 2 2 u ∈ CD [0, 1], (x) (x) dx + − x v(x) dx for all v ∈ CD [a, b], 2u(x)v(x) dx = dx 2 0 dx 0 0 and the Galerkin method requires the solution of Z 1 Z 1 Z 1 dvn 1 dv vn ∈ Vn , (x) (x) dx + 2 − x v(x) dx for all v ∈ Vn . vn (x)v(x) dx = dx dx 2 0 0 0 The second integral on the left side of the variational equation leads to the matrix M ∈ R(n−1)×(n−1) , where Z 1 Mij = φj (x)φi (x) dx 0

and φj is the jth standard basis function for the space of continuous piecewise linear functions. The matrix M can be computed explicitly (similarly to how K was computed in Example 5.18 in the text), and the result is  2h , j = i,  3 Mij = −f rach6, j = i − 1 or j = i + 1,  0, otherwise. We then must solve (K + 2M)u = f , where K is the usual stiffness matrix (as in Example 5.18) and the load vector f is defined by Z 1 1 h fi = − x φi (x) dx = (1 − 2ih), i = 1, 2, . . . , n − 1. 2 2 0 After solving (K + 2M)u = f to get the piecewise linear approximation for each n, we estimate the maximum error by evaluating both the exact function and the approximation on a finer mesh. The results are shown in the following table. n 10 20 40 80 160 320

error 5.4953 · 10−4 1.4660 · 10−4 3.7842 · 10−5 9.6121 · 10−6 2.4222 · 10−6 6.0794 · 10−7

5.6. PIECEWISE POLYNOMIALS AND THE FINITE ELEMENT METHOD

37

These results show that when n is doubled, the error is reduced by a factor of approximately 4. Thus it appears that the error is O(1/n2 ).

Chapter 6

Heat Flow and Diffusion 6.1

Fourier series methods for the heat equation

1. The steady-state solution of Example 6.2 satisfies the BVP d2 u dx2 u(0)

=

10−7 , 0 < x < 100,

=

0,

u(100)

=

0.

−D

Either by solving this BVP or by taking the limit (as t → ∞) of the solution of Example 6.2, we find that the steady-state solution is ∞ nπx X 2 · 10−3 (1 − (−1)n ) us (x) = sin . Dn3 π 3 100 n=1 3. The solution is u(x, t) =

∞ X 2(−1)n+1 −n2 π2 t e sin (nπx). nπ n=1

A graph of u(·, 0.1) is given in Figure 6.1. 1 u(⋅,0) u(⋅,0.1)

0.9 0.8 0.7 0.6 0.5 0.4 0.3 0.2 0.1 0 0

0.2

0.4

0.6

0.8

1

x

Figure 6.1: The snapshot u(·, 0.1), together with the initial temperature distribution, for Exercise 6.1.3. 5. With p(x, t) = g(t) +

x (h(t) − g(t)) `

and v(x, t) = u(x, t) − p(x, t), we have ρc

∂2v ∂v −κ 2 ∂t ∂x

= =

∂2u ∂2p ∂u ∂p − κ 2 − ρc −κ 2 ∂t ∂x ∂t ∂x x dh dg dg f (x, t) − ρc (t) + (t) − (t) . dt ` dt dt ρc

39

40

CHAPTER 6. HEAT FLOW AND DIFFUSION We also have v(x, t0 ) = u(x, t0 ) − p(x, t0 ) = ψ(x) − g(t0 ) −

x (h(t0 ) − g(t0 )) `

and v(0, t)

=

u(0, t) − p(0, t) = g(t) − g(t) = 0,

v(`, t)

=

u(`, t) − p(`, t) = h(t) − h(t) = 0.

We define

g(x, t) = f (x, t) − ρc



dg x (t) + dt `

and φ(x) = ψ(x) − g(t0 ) −

dh dg (t) − (t) dt dt

x (h(t0 ) − g(t0 )). `

Then v satisfies ρc

∂2v ∂v −κ 2 ∂t ∂x v(x, t0 )

=

φ(x), 0 < x < `,

v(0, t)

=

0, t > t0 ,

v(`, t)

=

0, t > t0 .

=

g(x, t), 0 < x < `, t > 0,

7. Define v(x, t) = u(x, t) − x cos (t). Then v satisfies the following IBVP (with homogeneous boundary conditions, but with a nonzero source term): ∂v ∂2v − ∂t ∂x2 v(x, 0)

=

x sin (t), 0 < x < 1, t > 0,

=

0, 0 < x < 1,

v(0, t)

=

0, t > 0,

v(1, t)

=

0, t > 0.

This IBVP has solution

∞ X

v(x, t) =

an (t) sin (nπx),

n=1

where

2 2 2(−1)n cos (t) − e−n π t − n2 π 2 sin (t) . 4 4 nπ (n π + 1) The solution to the original IBVP is then u(x, t) = v(x, t) + x cos (t). A graph of u(·, 1.0) is given in Figure 6.2. an (t) =

1 u(⋅,0) u(⋅,0.1)

0.9 0.8 0.7 0.6 0.5 0.4 0.3 0.2 0.1 0 0

0.2

0.4

0.6

0.8

1

x

Figure 6.2: The snapshot u(·, 1.0), together with the initial temperature distribution, for Exercise 6.1.7. 9. The temperature u(x, t) satisfies the IBVP ρc

∂u ∂2u −κ 2 ∂t ∂x u(x, 0)

=

0, 0 < x < 100, t > 0,

=

5, 0 < x < 100,

u(0, t)

=

0, t > 0,

u(100, t)

=

0, t > 0.

6.2. PURE NEUMANN CONDITIONS AND THE FOURIER COSINE SERIES The solution is u(x, t) =

41

∞ nπx X κn2 π 2 t 10(1 − (−1)n ) − 100 2 ρc sin e . nπ 100 n=1

The temperature at the midpoint after 20 minutes is . u(50, 1200) = 1.58 degrees Celsius. 11.

(a) The steady-state temperature is us (x) = x/20. (b) The temperature u(x, t) satisfies the IBVP ρc

∂2u ∂u −κ 2 ∂t ∂x u(x, 0)

=

0, 0 < x < 100, t > 0,

=

5, 0 < x < 100,

u(0, t)

=

0, t > 0,

u(100, t)

=

5, t > 0.

The solution is

∞

u(x, t) =

nπx X 10 − κn2 π2 t x + e 1002 ρc sin . 20 n=1 nπ 100

Considering the results of Exercise 6.1.9, it is obvious that at least several thousand seconds will elapse before the temperature is within 1% of steady state, so we can accurately estimate u(x, t) using a single term of the Fourier series: πx κπ 2 t 10 − 100 . x 2 ρc u(x, t) = + e . sin 20 π 100 We want to find t large enough that |u(x, t) − us (x)| ≤ 0.01 |us (x)| for all x ∈ [0, 100]. Using our approximation for u, this is equivalent to 200 sin πx κπ 2 100 − 1002 ρc t ≤ 0.01. e πx A graph shows that 200 sin πx

πx 100

≤2

for all x, so we need 2 − κπ2

e

100 ρc

t

≤ 0.005.

This yields −1002 ρc ln 0.005 . = 4520. κπ 2 About 75 minutes and 20 seconds are required. t≥

6.2

Pure Neumann conditions and the Fourier cosine series

1. The solution is u(x, t) = d0 +

∞ X

2

dn e−n

π2 t

cos (nπx),

n=1

where Z d0 = 0

1

1 x(1 − x) dx = , dn = 2 6

The graphs are given in Figure 6.3.

1

Z 0

2 (−1)n+1 − 1 x(1 − x) cos (nπx) dx = . n2 π 2

42

CHAPTER 6. HEAT FLOW AND DIFFUSION 0.25

0.2

0.15

0.1 u(x,0) u(x,0.02) u(x,0.04 u(x,0.06) u(x,∞)

0.05

0 0

0.2

0.4

0.6

0.8

1

x

Figure 6.3: The solution u(x, t) from Exercise 6.2.4 at times 0, 0.02, 0.04, and 0.06, along with the steady-state solution. These solutions were estimated using 10 in the Fourier series. 3.

2 [0, `], then (a) If u, v ∈ CN

= = = = = =

`

d2 u (x)v(x) dx 2 0 dx ` Z ` dv du du (x)v(x) + (x) (x) dx − dx dx 0 dx 0 Z ` du dv (x) (x) dx (since du (0) = du (`) = 0) dx dx dx dx 0 ` Z ` d2 v dv u(x) 2 (x) dx u(x) (x) − dx dx 0 0 Z ` 2 d v dv dv − (0) = dx (`) = 0) u(x) 2 (x) dx (since dx dx 0 (u, LN v). Z

(LN u, v)

−

This shows that LN is symmetric. (b) Suppose λ is an eigenvalue of LN and u is a corresponding eigenvector, normalized so that (u, u) = 1. Then λ = λ(u, u)

=

(λu, u)

=

(LN u, u) Z ` 2 d u − (x)u(x) dx 2 0 dx ` Z ` 2 du du − (x)u(x) + (x) dx dx dx 0 0 2 Z ` du (x) dx (since du (0) = du (`) = 0) dx dx dx 0 0.

= = = ≥

Thus LN cannot have any negative eigenvalues. 5.

(a) Suppose u is a solution to the BVP. Then Z ` Z ` 2 d u du du f (x) dx = −κ (x) dx = κ (0) − (`) = κ(a − b). 2 dx dx 0 0 dx This is the compatibility condition: Z

`

f (x) dx = κ(a − b). 0 2 (b) The operator K : CN [0, `] → C[0, `] defined by

Ku = −κ

d2 u +u dx2

6.2. PURE NEUMANN CONDITIONS AND THE FOURIER COSINE SERIES

43

has eigenpairs nπx κn2 π 2 , n = 1, 2, 3, . . . . , γ (x) = cos n `2 `

λ0 = 1, γ0 (x) = 1, and λn = 1 +

The method of Fourier series can be applied to show that a unique solution exists for each f ∈ C[0, `]. (The key is that 0 is not an eigenvalue of K, as it is of LN .) 7. As shown in this section, the solution to the IBVP is u(x, t) = d0 +

∞ X

dn e−n

2

π 2 t/`2

cos

nπx `

n=1

,

where d0 , d1 , d2 , . . . are the Fourier cosine coefficients of ψ. We have ∞ ∞ X nπx nπx X 2 2 2 −n2 π 2 t/`2 cos dn e |dn | e−n π t/` cos ≤ ` ` n=1

≤

n=1 ∞ X

|dn | e−n

2

π 2 t/`2

n=1

≤

e−π

2

t/`2

∞ X

|dn | e−(n

2

−1)π 2 t/`2

.

n=1

This last series certainly converges (as can be proved, for example, using the comparison test), and e−π This shows that d0 +

∞ X

dn e−n

2

2

t/`2

→ 0 as t → ∞.

π 2 t/`2

cos

n=1

The limit d0 is 1 ` the average of the initial temperature distribution.

Z

nπx `

→ d0 as t → ∞.

`

ψ(x) dx, 0

9. The steady-state temperature is about 0.992 degrees Celsius. 11.

(a) Suppose that the Fourier sine series of u(x, t) on (0, `) is u(x, t) =

∞ X n=1

an (t) sin

Z nπx nπx 2 ` , an (t) = u(x, t) sin dx ` ` 0 `

and u satisfies

∂u ∂u (0, t) = (`, t) = 0 for all t > 0. ∂x ∂x The nth Fourier sine coefficient of −∂ 2 u/∂x2 is computed as follows: Z nπx 2 ` ∂2u − (x, t) sin dx ` 0 ∂x2 ` " # Z nπx ` nπx 2 ∂u nπ ` ∂u = − (x, t) sin − (x, t) cos dx ` ∂x ` 0 ` 0 ∂x ` Z nπx 2nπ ` ∂u = (x, t) cos dx `2 0 ∂x ` Z nπx ` nπx 2nπ nπ ` = u(x, t) cos + u(x, t) sin dx `2 ` ` 0 ` 0 =

2nπ n2 π 2 ((−1)n u(`, t) − u(0, t)) + 2 an (t). 2 ` `

Since the values u(0, t) and u(`, t) are unknown, we see that it is not possible to express the Fourier sine coefficients of −∂ 2 u/∂x2 in of a1 (t), a2 (t), . . ..

44

CHAPTER 6. HEAT FLOW AND DIFFUSION (b) The Fourier sine series of u(x, t) = t is ∞ X 2 (1 − (−1)n ) t sin (nπx), nπ n=1

and the formal calculation of the sine series of −∂ 2 u/∂x2 yields ∞ X

2 (1 − (−1)n ) nπt sin (nπx).

n=1

However, ∂2u (x, t) = 0, ∂x2 and so all of the Fourier sine coefficients of −∂ 2 u/∂x2 should be zero. Thus the formal calculation is wrong. −

6.3

Periodic boundary conditions and the full Fourier series

1. The formula for the solution u is exactly the same as in Example 6.6, with a different value for κ (4.29 instead of 3.17). This implies that the amplitude of the solution is reduced by about 26%. Therefore, there is less variation in the temperature distribution in the silver ring as opposed to the gold ring. 3.

(a) The IBVP is ρc

∂2u ∂u −κ 2 ∂t ∂x u(x, 0) u(−5π, t) ∂u (−5π, t) ∂x

=

0, −5π < x < 5π, t > 0,

=

ψ(x), −5π < x < 5π,

=

u(5π, t), t > 0, ∂u (5π, t), t > 0. ∂x

=

(b) The solution is u(x, t) = c0 +

∞ X

2

cn e

κn t − 25ρc

cos

n=1

nx 5

,

where c0

=

cn

=

π4 , 9 n+1 10(−1) , n = 1, 2, 3, . . . . n4

25 +

(c) The steady-state temperature is the constant us = 25 + π 4 /9 degrees Celsius. (d) We must choose t so that |us − u(x, t)| ≤ 0.01 |us | holds for every x ∈ [−5π, 5π]. By trial and error, we find that about 360 seconds (6 minutes) are required. 5.

(a) To show that Lp is symmetric, we perform the now familiar calculation: we form the integral (Lp u, v) and integrate by parts twice to obtain (u, Lp v). The boundary term from the first integration by parts is −

du (x)v(x) dx

` = −`

du du (−`)v(−`) − (`)v(`). dx dx

Since both u and v satisfy periodic boundary conditions, we have v(−`) = v(`),

du du (−`) = (`), dx dx

so the boundary term vanishes. The boundary term from the second integration by parts vanishes for exactly the same reason.

6.3. PERIODIC BOUNDARY CONDITIONS AND THE FULL FOURIER SERIES

45

(b) Suppose Lp u = λu, where u has been normalized: (u, u) = 1. Then Z ` 2 d u (x)u(x) dx λ = λ(u, u) = (λu, u) = (Lp u, u) = − 2 −` dx ` 2 Z ` du du + = − (x)u(x) (x) dx dx dx −` −` 2 Z ` du dx (x) = dx −` ≥ 0. The boundary vanishes because of the periodic boundary conditions: du du (−`)u(−`) = (`)u(`). dx dx 7.

(a) Since the ring is completely insulated, a steady-state temperature distribution cannot exist unless the net amount of heat being added to the ring is zero. This is exactly the same situation as a straight bar with the ends, as well as the sides, insulated. (b) Suppose u is a solution to (6.21). Then ` Z ` Z ` 2 ∂u ∂u ∂u ∂ u (x, t) dx = −κ (x, t) =κ (−`, t) − (`, t) = 0. f (x) dx = −κ 2 ∂x ∂x ∂x −` −` ∂x −` The last step follows from the periodic boundary conditions. (c) The negative second derivative operator, subject to boundary conditions, has a nontrivial null space, namely, the space of all constant functions on (−`, `). In analogy to the Fredholm alternative for symmetric matrices, we would expect a solution to the boundary value problem to exist if and only if the right-hand-side function is orthogonal to this null space. This condition is Z ` cf (x) dx = 0 for all c ∈ R, −`

or simply Z

`

f (x) dx = 0. −`

9. Consider the IBVP ρc

∂2u ∂u − κ 2 + pu = f (x, t), −` < x < `, t > t0 , ∂t ∂x u(x, t0 ) = ψ(x), −` < x < `, u(−`, t) = u(`, t), t > t0 , ∂u ∂u (−`, t) = (`, t), t > t0 . ∂x ∂x

Assume that p is a constant. (a) We wish to derive the solution to the IBVP using the Fourier series method. The calculation is very similar to that carried out in Section 6.3.3 in the text, and we will use the same notation as in that section. We represent the solution u as the series ∞ n nπx nπx o X u(x) = a0 + an cos + bn sin . ` ` n=1 Substituting this series into the left side of the PDE yields ∂u ∂2u (x, t) − κ 2 (x, t) + pu(x, t) ∂t ∂x 2 2 ∞ nπx X da0 dan κn π = ρc + pa0 (t) + ρc (t) + + p a (t) cos + n dt dt `2 ` n=1 2 2 nπx dbn κn π ρc (t) + + p b (t) sin n dt `2 ` ρc

46

CHAPTER 6. HEAT FLOW AND DIFFUSION Following the derivation in Section 6.3.3, we find a0 , a1 , . . ., b1 , b2 , . . . by solving the IVPs da0 + pa0 = c0 (t), a0 (t0 ) = p0 , dt 2 2 κn π + + p an = cn (t), an (t0 ) = pn , n = 1, 2, . . . , `2 2 2 κn π + + p bn = dn (t), bn (t0 ) = qn , n = 1, 2, . . . . `2 ρc

ρc

dan dt

ρc

dbn dt

(b) The eigenvalues of the spatial operator are λ0 = p and κn2 π 2 + p, n = 1, 2, . . . . `2

λn =

Therefore, λn ≥ λ0 for all n, and all the eigevalues are positive if and only if p > 0. 11. Let f : (−`, `) → R be an even function. We wish to prove that the full Fourier series of f reduces to the cosine series of f , regarded as a function on (0, `). The full Fourier series of f is f (x) = a0 +

∞ n X

an cos

nπx

n=1

+ bn sin

`

nπx o `

,

where Z ` 1 f (x) dx, 2` −` Z ` nπx 1 dx, n = 1, 2, . . . , f (x) cos an = ` −` ` Z nπx 1 ` bn = dx, n = 1, 2, . . . . f (x) sin ` −` ` a0 =

For any function g : (−`, `) → R, if g is even, then Z ` Z g(x) dx = 2 −`

`

g(x) dx, 0

while if g is odd, then Z

`

g(x) dx = 0. −`

Since f is even, a0 =

1 2`

Z

`

f (x) dx = −`

1 `

Z

`

f (x) dx. 0

Also, f (x) cos (nπx/`) is an even function of x (the product of even functions is even), and therefore Z Z nπx nπx 2 ` 1 ` an = f (x) cos dx = f (x) cos dx, n = 1, 2, . . . . ` −` ` ` 0 ` Finally, f (x) sin (nπx/`) is an odd function of x (the product of an even and an odd function is odd), and thus Z nπx 1 ` dx = 0, n = 1, 2, . . . . bn = f (x) sin ` −` ` Thus the full Fourier series reduces to f (x) = a0 +

∞ X n=1

where

an cos

nπx `

,

Z Z nπx 1 ` 2 ` f (x) dx, an = f (x) cos dx, n = 1, 2, . . . . ` 0 ` 0 ` This is precisely the Fourier cosine series of f on the interval (0, `) (cf. Section 6.2.2). a0 =

6.4. FINITE ELEMENT METHODS FOR THE HEAT EQUATION

6.4

47

Finite element methods for the heat equation

1. We give the proof for the general case of a Gram matrix G. Suppose Gx = 0, where x ∈ Rn . Then (x, Gx) = 0 must hold, and ! n n X n n X n n n X X X X X (x, Gx) = (Gx)i xi = Gij xj xi = (uj , ui )xj xi = xj u j , u i u i i=1

i=1 j=1

i=1 j=1

i=1 n X

=

j=1

xj uj ,

j=1

=

n X

! xi u i

i=1

n

2

X

xj uj .

j=1

Pn

Thus Gx = 0 implies that the vector j=1 xj uj is the zero vector. But, since {u1 , u2 , . . . , un } is linearly independent, this in turn implies that x1 = x2 = · · · = xn = 0, that is, that x = 0. Since the only vector x ∈ Rn satisfying Gx = 0 is the zero vector, G is nonsingular. 3.

(a) M=

2/9 1/18

1/18 2/9





6 −3

, K=

−3 6

, f (t) =

1 162



11 cos (t) 11 cos (t)

(b) The system of ODEs is M dα = −Kα + f (t), that is, dt 2 dα1 1 dα2 + 9 dt 18 dt 2 dα2 1 dα1 + 18 dt 9 dt 5.

= =

11 cos (t) , 162 11 cos (t) 3α1 − 6α2 + . 162 −6α1 + 3α2 +

(a) Write ρ1 = 8.97, ρ2 = 7.88, ρ(x) =

ρ1 , ρ2 ,

0 < x < 50, 50 < x < 100,

and similarly for c(x) and κ(x). Then the IBVP is ∂u ∂ ∂u ρ(x)c(x) − κ(x) ∂t ∂x ∂x u(x, 0)

=

0, 0 < x < 100, t > 0,

=

5, 0 < x < 100,

u(0, t)

=

0, t > 0,

u(100, t)

=

0, t > 0.

(b) The mass matrix M is tridiagonal and symmetric, and its nonzero entries are ( ρ c h 1 1 , i = 1, 2, . . . , n2 − 1, 6 Mi,i+1 = ρ2 c 2 h , i = n2 , n2 + 1, . . . , n − 1, 6 and Mii =

      

2ρ1 c1 h , 3 (ρ1 c1 +ρ2 c2 )h , 3 2ρ2 c2 h , 3

i = 1, 2, . . . , i= i=

n , 2 n + 2

1,

n 2

n 2

− 1,

+ 2, . . . , n − 1.

The stiffness matrix K is tridiagonal and symmetric, and its nonzero entries are ( κ − h1 , i = 1, 2, . . . , n2 − 1, Ki,i+1 = − κh2 , i = n2 , n2 + 1, . . . , n − 1, and Kii =

    

2κ1 , h κ1 +κ2 , h 2κ2 , h

i = 1, 2, . . . , i= i=

n , 2 n + 2

1,

n 2

n 2

− 1,

+ 2, . . . , n − 1.

.

48

CHAPTER 6. HEAT FLOW AND DIFFUSION 7. The solution is u(x, t) =

∞ X

an (t) sin

n=1

nπx , 100

where an (t) =

400(2 + (−1)n ) κ2 π 7 n7

2 2 − κn π t 104 ρc e 104 ρc − 1 + κn2 π 2 t .

The errors in Examples 6.8 and 6.9 are graphed in Figure 6.4. 0.05

0

−0.05

−0.1

−0.15

−0.2 Euler Backward Euler −0.25 0

20

40

60

80

100

x

Figure 6.4: The errors in Examples 6.8 and 6.9 (see Exercise 6.4.7).

(a) The IBVP is

ρ(x)c(x)

∂ ∂u − ∂t ∂x

∂u κ(x) ∂x u(x, 0)

=

0, 0 < x < `, t > 0,

=

5, 0 < x < `,

u(0, t)

=

0, t > 0,

u(`, t)

=

0, t > 0.

(b) The weak form is to find u satisfying `

Z 0

∂u ∂u dv ρ(x)c(x) (x, t)v(x) + κ(x) (x, t) (x) dx = 0, t > 0, ∀v ∈ V. ∂t ∂x dx

(c) The temperature distribution after 120 seconds is shown in Figure 6.5. t = 120 (seconds) 5 4.5 4 3.5

temperature

9.

3 2.5 2 1.5 1 0.5 0 0

20

40

60

80

100

x

Figure 6.5: The temperature distribution after 120 seconds in Exercise 6.4.9.

6.5. FINITE ELEMENTS AND NEUMANN CONDITIONS

6.5

49

Finite elements and Neumann conditions

1. The BVP to be solved is d2 u dx2 du (0) dx du (100) dx −k

=

f (x), 0 < x < 100,

=

0,

=

0,

with k = 1.5 and f (x) = 10−7 x(25 − x)(100 − x) + 1/240. The steady-state temperature is not unique; the solution with u(100) = 0 is shown in Figure 6.6. 6 5

temperature

4 3 2 1 0 −1 0

20

40

60

80

100

x

Figure 6.6: The computed steady-state temperature distribution from Exercise 6.5.1. (a) The total amount of heat energy being added to the bar is 0.51A W, where A is the cross-sectional area (0.01A W through the left end and 0.5A W in the interior). Therefore, 0.51A W must be removed through the right end; that is, heat energy must be removed at a rate of 0.51 W/cm2 through the right end. (b) The BVP is d2 u dx2 du k (0) dx du k (100) dx −κ

=

0.005, 0 < x < 100,

=

−0.01,

=

−0.51.

The steady-state temperature is not unique; the temperature with u(100) = 0 is graphed in Figure 6.7. 7 6 5

temperature

3.

4 3 2 1 0 0

20

40

60

80

100

x

Figure 6.7: The computed steady-state temperature distribution from Exercise 6.5.3.

50

CHAPTER 6. HEAT FLOW AND DIFFUSION 5.

(a) We have

˜ = u · Ku

n X n X

n X n X

˜ ij uj ui = K

i=0 j=0

a (φj , φi ) uj ui

=

i=0 j=0

n X

a

i=0

=

a

n X

ui

n X

!

j=0 n X

uj φj ,

j=0

=

! uj φj , φ i

ui φi

i=0

a(v, v).

Therefore, ˜ = 0 ⇒ u · Ku ˜ = 0 ⇒ a(v, v) = 0. Ku (b) Suppose v ∈ V˜ and a(v, v) = 0. By definition of a(·, ·), this is equivalent to `

Z

κ(x)

0

2 dv (x) dx = 0. dx

Since the integrand is nonnegative, this implies that the integrand is in fact zero, and, since κ(x) is positive, we conclude that dv (x) dx

7.

is the zero function. Therefore, v is a constant function. ˜ = 0, it follows that v(x) = Pn ui φi (x) is a constant function, where u0 , u2 , . . . , un are the (c) Thus, if Ku i=0 components of u. But then there is a constant C such that v(xi ) = C, i = 0, 1, 2, . . . , n. These nodal values of v are precisely the numbers u0 , u1 , . . . , un , so we see that u = Cuc , which is what we wanted to prove. (a) Define Vˆ = v ∈ C 2 [0, `] : v(0) = 0 . Then the weak form of the BVP is find u ∈ Vˆ such that a(u, v) = (f, v) for all v ∈ Vˆ .

(6.1)

(b) The fact that a solution of the strong form is also a solution of the weak form is proved by the usual argument: multiply the differential equation by an arbitrary test function v ∈ Vˆ , and then integrate by parts. The proof that a solution of the weak form is also a solution of the strong form is similar to the argument given in Section 6.5.2. Assuming that u satisfies the weak form (6.1), an integration by parts and some simplification shows that

k(`)

du (`)v(`) − dx

Z 0

`



du d k(x) (x) + f (x) v(x) dx = 0 for all v ∈ Vˆ . dx dx

Since V ⊂ Vˆ , this implies that the differential equation d du k(x) (x) + f (x) = 0, 0 < x < ` dx dx holds, and we then have k(`)

du (`)v(`) = 0 for all v ∈ Vˆ . dx

Choosing any v ∈ Vˆ with v(`) 6= 0 shows that the Neumann condition holds at x = `. 9. The temperature distribution after 300 seconds is shown in Figure 6.8.

6.5. FINITE ELEMENTS AND NEUMANN CONDITIONS

51

t=300 (seconds) 8 7

temperature

6 5 4 3 2 1 0 0

20

40

60

80

100

x

Figure 6.8: The temperature distribution from Exercise 6.5.9 (after 300 seconds). ˆ = 0, then 11. Here is a sketch of the proof: If Ku ˆ = 0 ⇒ a(v, v) = 0, u · Ku P where v = n−1 i=0 ui φi . But then, by the usual reasoning, v must be a constant function, and v(xn ) = 0 (since φi (xn ) = 0 for i = 0, 1, 2, . . . , n − 1). Thus v is the zero function, which implies that the nodal values of v are all ˆ is nonsingular. zero. Therefore u = 0, and so K

Chapter 7

Waves 7.1

The homogeneous wave equation without boundaries

1. The solution, by d’Alembert’s formula, is u(x, t) = (ψ(x − ct) + ψ(x + ct)). Figure 7.1 shows three snapshots of u. −4

20

x 10

15

10

5

0

−5 −20

−15

−10

−5

0

5

10

15

20

Figure 7.1: Three snapshots of the solution of Exercise 7.1.1: t = 0 (dashed curve), t = 0.01 (solid curve), and t = 0.025 (dash-dotted curve).

3. The two waves and add constructively, and then separate again. See Figure 7.2. −3

1.5

x 10

t=0 t=0.02 t=0.04

1

0.5

0 −40

−30

−20

−10

0 x

10

20

30

40

Figure 7.2: Constructive interference in Exercise 7.1.3. The “blip” in the center is the result of two waves combining temporarily. (The velocity in this example is c = 1.)

53

54

CHAPTER 7. WAVES 5. Consider the IVP 2 ∂2u 2∂ u − c = 0, −∞ < x < ∞, t > 0, ∂t2 ∂x2 u(x, 0) = ψ(x), −∞ < x < ∞,

∂u = 0, −∞ < x < ∞, ∂t where the of ψ is [a, b]. The solution is u(x, t) = (ψ(x − ct) + ψ(x + ct))/2. Now, since ψ(x − ct) is a right-moving wave, we see that ψ(x1 − ct) = 0 if any of the following conditions hold: • x1 ≤ a; • a < x1 < b and t ≥

x1 −a c

(since ψ(x1 − ct) = 0 after the trailing edge of the wave es x1 );

x1 −b c

• x ≥ b and (t ≤ or t ≥ x1c−a ) (since ψ(x1 − ct) = 0 before the leading edge of the wave reaches x1 and after the trailing edge of the wave es x1 ). Similarly, ψ(x + ct) is a left-moving wave, and therefore ψ(x1 + ct) = 0 if any of the following is true: • x1 ≥ b; • a < x1 < b and t ≥ • x1 ≤ a and (t ≤

b−x1 ; c

a−x1 c

or t ≥

b−x1 ). c

Since u(x1 , t) is guaranteed to be zero if both ψ(x1 − ct) and ψ(x1 + ct) are zero, we see that u(x1 , t) = 0 if any of the following is true: • x1 ≤ a and (t ≤

a−x1 c

• a < x1 < b and t ≥ • x1 ≥ b and (t ≤

or t ≥

x1 −a c

x1 −b c

b−x1 ); c

and t ≥

or t ≥

b−x1 ; c

x1 −a ). c

7. Consider the IBVP ∂2u ∂2u − c2 2 = 0, 0 < x < `, t > 0, 2 ∂t ∂x u(x, 0) = ψ(x), 0 < x < `, ∂u = 0, 0 < x < `, ∂t u(0, t) = 0, t > 0, u(`, t) = 0, t > 0, where the of ψ is [a, b] and 0 < a < b < `. If we define u(x, t) =

1 (ψ(x − ct) + ψ(x + ct)) , 2

then we know that u satisfies the same PDE and initial condition as does the solution of the above IBVP. Also, since the of ψ is [a, b], u(0, t) = 0 for all t ≤ a/c (since u(0, t) cannot be nonzero until the leading edge of the left-moving wave reaches x = 0), and similarly u(`, t) = 0 for all t ≤ (` − b)/c (since u(`, t) cannot be nonzero until the leading edge of the right-moving wave reaches x = `). Therefore, if t < tf = min{a/b, (` − b)/c}, then u satisfies both boundary conditions, in addition to the initial conditions and the PDE, and hence is a solution of the above IBVP.

7.2

Fourier series methods for the wave equation

1. The solution is found by setting cn = 0 in Example 7.4: u(x, t) =

∞ X n=1

bn cos

cnπt 25

sin

nπx . 25

A graph analogous to Figure 7.6 from the text is given in Figure 7.3.

displacement

7.2. FOURIER SERIES METHODS FOR THE WAVE EQUATION

55

0

0

5

10

15

20

25

x

Figure 7.3: Fifty snapshots of the vibrating string from Example 7.4, with no external force.

3. The solution is u(x, t) =

∞ X

an (t) sin

n=1

(2n − 1)πx , 50

where 502 cn 2 c (2n − 1)2 π 2

an (t) = bn − and bn =

8



cos

c(2n − 1)πt 50

+

502 cn − 1)2 π 2

c2 (2n

√ √ 2 sin (nπ/2) − 2 cos (nπ/2) + (−1)n 4000 , cn = . 5(2n − 1)2 π 2 (2n − 1)π

The fundamental frequency is now c/100 instead of c/50. Fifty snapshots are shown in Figure 7.4. 2.5 2 1.5

displacement

1 0.5 0 −0.5 −1 −1.5 −2 −2.5 0

5

10

15

20

25

x

Figure 7.4: Fifty snapshots of the vibrating string from Example 7.4, with the right end free.

5. The solution is u(x, t) =

1 + 50t2 20

+

∞ X

bn cos

n=1

cnπt 25

cos

nπx 25

,

where bn =

2 (2 cos (nπ/2) − 1 − (−1)n ) , n = 1, 2, 3, . . . . 5n2 π 2

Fifty snapshots of the solution are shown in Figure 7.5. The solution gradually moves up; this is possible because both ends are free to move vertically. 7. The solution is u(x, t) =

∞ X n=1

an (t) sin

(2n − 1)πx , 2

56

CHAPTER 7. WAVES 0.6

displacement

0.4

0.2

0 0

5

10

15

20

25

x

Figure 7.5: Fifty snapshots of the vibrating string from Example 7.4, with both ends free.

where an (t)

=

bn

=

cn

=

bn −

4cn c2 (2n − 1)2 π 2



cos

c(2n − 1)πt 2

+

4cn , c2 (2n − 1)2 π 2

2(−1)n+1 , 125(2n − 1)2 π 2 4g . (2n − 1)π

Snapshots of the solution are graphed in Figure 7.6, which should be compared to Figure 7.9 from the text. −3

displacement

2

x 10

1 t=0 t=5e−05 t=0.0001 t=0.00015

0 0

0.2

0.4

0.6

0.8

1

x

Figure 7.6: Snapshots of the displacement of the bar in Example 7.6, taking into the effect of gravity. 9. Consider a string whose (linear) density (in the unstretched state) is 0.25 g/cm and whose longitudinal stiffness is 6500 N. Suppose the unstretched (zero tension) length of the string is 40 cm. We wish to determine the length to which the string must be stretched in order that its fundamental frequency be f0 = 261 Hz (middle C). We will write `0 = 40, k = 6.5 · 108 (the stiffness of the string in g cm/s2 ), and T for the unknown tension. Note that the mass of the string is m = 10 (g). We know that f0 = c/(2`), where c is the wave speed and ` = `0 + ∆`. Also, c2 = T /ρ, where ρ = m/`. Finally, ∆` is related to T by Hooke’s law: k∆`/`0 = T . Putting all these equations together, we obtain q c 1 q f0 = = 2(`0 + ∆`) 2(`0 + ∆`)

k∆` `0

m `0 +∆`

√ k∆` k∆` = p ⇒ = f02 . 4m` (` 0 0 + ∆`) 2 m`0 (`0 + ∆`)

Some straightforward algebra yields ∆` =

4m`20 f02 . = 8.0586. k − 4m`0 f02

Thus the string must be stretched to 48.0586 cm.

7.3. FINITE ELEMENT METHODS FOR THE WAVE EQUATION

7.3

57

Finite element methods for the wave equation

1. The fundamental period is 2`/c = 0.2. Using h = 100/20 = 5, Dt = T /60, and the RK4 method, we obtain the solution shown in Figure 7.7.

displacement

0.2

0

−0.2 0

0.2

0.4

0.6

0.8

1

x

Figure 7.7: The computed solution of the wave equation in Exercise 7.3.1. Shown are 30 snapshots, plus the initial displacement. Twenty subintervals in space and 60 time steps were used. 3.

(a) Since the initial disturbance is 24 cm from the boundary and the wave speed is 400 cm/s, it will take 24/400 = 0.06 s for the wave to reach the boundary. (b) The IBVP is ∂2u ∂2u − c2 2 2 ∂t ∂x u(x, 0) ∂u (x, 0) ∂t u(0, t) u(50, t)

=

0, 0 < x < 50, t > 0,

=

0, 0 < x < 50,

=

γ(x), 0 < x < 50,

=

0, t > 0,

=

0, t > 0,

with c = 400 and γ given in the statement of the exercise. (c) Using piecewise linear finite elements and the RK4 method (with h = 50/80 and ∆t = 6·10−4 ), we computed the solution over the interval 0 ≤ t ≤ 0.06. Four snapshots are shown in Figure 7.8, which shows that it does take 0.06 seconds for the wave to reach the boundary. (The reader should notice the spurious “wiggles” in the computed solution; these are due to the fact that the true solution is not smooth.)

displacement

0.1

0

−0.1 0

10

20

30

40

50

x

Figure 7.8: The computed solution of the wave equation in Exercise 7.3.3. 5. It is not possible to obtain a reasonable numerical solution using finite elements. In Figure 7.9, we display the result obtained using h = 2.5 · 10−3 and ∆t = 3.75 · 10−4 . Three snapshots are shown. 7.

(a) The weak form is Z 0

`

Z ` ∂u dv ∂2u ρ(x) 2 (x, t)v(x) + k(x) (x, t) (x) dx = f (x, t)v(x) dx ∂t ∂x dx 0

58

CHAPTER 7. WAVES 1.2 1

displacement

0.8 0.6 0.4 0.2 0 −0.2 0

0.2

0.4

0.6

0.8

1

x

Figure 7.9: The computed solution of the wave equation in Exercise 7.3.5. for t ≥ t0 and v ∈ V˜ , where V˜ = v ∈ C 2 [0, `] : v(0) = 0 . (b) The system of ODEs has the same form as in the case of a homogeneous bar: 2

˜ d u + Ku ˜ =˜ M f (t), dt2 except that the entries in the mass matrix are now Z ` ˜ ij = M ρ(x)φj (x)φi (x) dx, i, j = 1, 2, . . . , n, 0

and the entries in the stiffness matrix are now Z ` dφi dφj ˜ ij = (x) (x) dx, i, j = 1, 2, . . . , n. K k(x) dx dx 0 The components of the vector ˜ f are unchanged (see Section 7.3.2).

7.4

Resonance

1. We represent the solution as u(x, t) =

∞ X n=1

an (t) sin

nπx . `

Following the usual procedure, an must solve the IVP d2 an cnπ 2 dan + an = cn (t), an (0) = 0, (0) = 0, dt2 ` dt where cn (t), n = 1, 2, . . ., are the Fourier sine coefficients of the right-hand-side function f (t) = sin (2πωt): Z nπx 2 ` 2 sin (2πωt) sin dx = cn (t) = (1 − (−1)n ) sin (2πωt) ` 0 ` nπ 0, if n is odd, = 4 sin (2πωt), if n is even. nπ Therefore, an (t) = 0 if n is even, while, if n is odd, Z t cnπ ` an (t) = sin (t − s) cn (s) ds. cnπ 0 ` If n is odd and ω 6= cn/(2`), we obtain an (t) =

4`2 (cn sin (2πωt) − 2`ω sin cn2 π 2 (c2 n2 − 4`2 ω 2 )

cnπt `

)

.

7.4. RESONANCE

59

If n is odd and ω = cn/(2`), we obtain an (t) =

cnπt `

2`(` sin



− cnπt cos c2 n3 π 3

cnπt `

)

.

Therefore, if ω does not equal cn/(2`) for any odd n, then the Fourier coefficients of u are uniformly bounded with respect to t and go to zero like 1/n3 . Notice, though, that if ω is very close to cn/(2`) for some odd n, then the corresponding an (t) is larger (due to the factor of c2 n2 − 4`2 ω 2 in the denominator), and thus that frequency has greater weight in the Fourier series. On the other hand, if ω = cn/(2`) for some odd n, then the corresponding Fourier coefficient an (t) grows without bound as t → ∞, and resonance is observed. 3. The experiment described in this problem is modeled by the IBVP ∂2u ∂2u − c2 2 = 0, 0 < x < 1, t > 0, ∂t2 ∂x u(x, 0) = 0, 0 < x < 1, ∂u (x, 0) = 0, 0 < x < 1, ∂t ∂u k (0, t) = B sin (2πωt), t > 0, ∂x u(1, t) = 0, t > 0. . Here c2 = k/ρ = 3535.53 m/s. We transform this problem into one that we can analyze by the Fourier series method by shifting the data. The function v(x, t) = (B/k) sin (2πωt)(x − 1) satisfies the given boundary conditions; if we define w = u − v, then w satisfies the IBVP ∂2w ∂2w 4π 2 ω 2 B sin (2πωt)(x − 1), 0 < x < 1, t > 0, − c2 2 = 2 ∂t ∂x k w(x, 0) = 0, 0 < x < 1, ∂w 2πωB (x, 0) = − (x − 1), 0 < x < 1, ∂t k ∂w k (0, t) = 0, t > 0, ∂x w(1, t) = 0, t > 0. We solve this IBVP by the usual method. We write u(x, t) =

∞ X

an (t) sin (nπx).

n=1

The Fourier coefficients an are found by solving the IVP d2 an dan + c2 n2 π 2 an = cn (t), an (0) = 0, (0) = bn , dt2 dt where Z

1

cn (t) = 2 0

4π 2 ω 2 B 8πω 2 B sin (2πωt)(x − 1) sin (nπx) dx = − sin (2πωt) k kn

and bn = −

4πωB k

Z

1

(x − 1) sin (nπx) dx = 0

4ωB . kn

It follows (see Section 4.2.3 from the text) that Z t Z t bn 1 4ωB 1 an (t) = sin (cnπt) + sin (cnπ(t − s))cn (s) ds = sin (cnπt) + sin (cnπ(t − s))cn (s) ds. cnπ cnπ 0 kcn2 π cnπ 0 If ω 6= cn/2, then an (t) =

8ω 2 B(cn sin (2πωt) − 2ω sin (cnπt)) 4ωB sin (cnπt) − , 2 kcn π ckn2 π(c2 n2 − 4ω 2 )

60

CHAPTER 7. WAVES while if ω = cn/2, an (t) =

4ωB 4ω 2 B sin (cnπt) − 2 kcn π ckn2 π



sin (cnπt) − t cos (cnπt) . cnπ

(a) From the above analysis, we see that the smallest resonant frequency is ωr = c/2, the fundamental frequency of the wave equation modeling the problem. (b) If ω = ωr , then a1 (t) is given by the second formula above, while an (t), n > 1, is given by the first. We compute several snapshots of the displacement, using 20 in the Fourier series, on each of the time intervals [0, 0.005], [0.005, 0.01], [0.01, 0.015], and [0.015, 0.02]. These are shown in Figure 7.10. The graphs show that the amplitude grows as t increase, and also that the shape of the displacement is dominated by the first fundamental mode. 0.5

0.5

0

0

−0.5 0

0.5

1

−0.5 0

0.5

0.5

0

0

−0.5 0

0.5

1

−0.5 0

0.5

1

0.5

1

Figure 7.10: Resonance in Exercise 7.3.3(a). (c) With ω = ωr /4, we compute the displacements at the same times as in the previous part of the exercise. The results are shown in Figure 7.11. Here the absence of resonance is clearly seen, in that the amplitudes remain bounded as t increases. −3

5

−3

x 10

5

0

0

−5 0

0.5

1

−5 0

−3

5

x 10

x 10

5

0

−5 0

0.5

1

0.5

1

−3

x 10

0

0.5

1

−5 0

Figure 7.11: Absence of resonance in Exercise 7.3.3(b).

7.5. FINITE DIFFERENCE METHODS FOR THE WAVE EQUATION

61

5. We solve the IBVP by the usual method, taking u(x, t) =

∞ X

an (t) sin (nπx),

n=1

where the Fourier coefficient an solves the IVP dan d2 an + c2 n2 π 2 an = cn (t), an (0) = 0, (0) = 0. dt2 dt Here cn is the Fourier coefficient of the right-hand-side function f (x, t): Z 1 Z 13/5 cn (t) = 2 f (x, t) sin (nπx) dx = 2 sin (2πωt) sin (nπx) dx 0

3/5

2 sin (2πωt) = nπ



cos

3nπ 5



− cos

13nπ 20

.

It follows that 1 cnπ

t

Z

sin (cnπ(t − s))cn (s) ds − cos 13nπ 2 cos 3nπ (cn sin (2πωt) − 2ω sin (cnπt)) 5 20 = . cn2 π 2 π(c2 n2 − 4ω 2 )

an (t) =

0

Five snapshots of the solution, at t = 0.002, 0.004, 0.006, 0.008, 0.01, are shown in Figure 7.12. −6

2

x 10

1.5 1 0.5 0 −0.5 −1 −1.5 −2 0

0.2

0.4

0.6

0.8

1

Figure 7.12: Snapshots of the solution of Exercise 7.3.5.

7.5

Finite difference methods for the wave equation

1. Let φ : R → R be at least twice continuously differentiable. Then, by Taylor’s theorem, φ(x + ∆x) = φ(x) +

dφ 1 d2 φ (x)∆x + (c)∆x2 , dx 2 dx2

where c is some number between x and x + ∆x. Solving for dφ/dx yields φ(x + ∆x) − φ(x) dφ 1 d2 φ (x) = − (c)∆x dx ∆x 2 dx2 or

Since

φ(x + ∆x) − φ(x) dφ 1 d2 φ (x) − = (c)∆x. dx ∆x 2 dx2

d2 φ d2 φ (c) → (x) 2 dx dx2 as ∆x → 0, it follows that the error in the forward difference approximation to the first derivative is approximately proportional to ∆x.

62

CHAPTER 7. WAVES 3. Suppose φ : R → R is at least four times continuously differentiable. By Taylor’s theorem, dφ (x)∆x + dx dφ (x)∆x + φ(x − ∆x) = φ(x) − dx

φ(x + ∆x) = φ(x) +

1 d2 φ (x)∆x2 + 2 dx2 1 d2 φ (x)∆x2 − 2 dx2

1 d3 φ (x)∆x3 + 6 dx3 1 d3 φ (x)∆x3 + 6 dx3

1 d4 φ (c1 )∆x4 , 24 dx4 1 d4 φ (c2 )∆x4 , 24 dx4

where c1 is between x and x + ∆x and c2 is between x and x − ∆x. Adding these two equations yields d2 φ 1 d4 φ d4 φ 2 φ(x + ∆x) + φ(x − ∆x) = 2φ(x) + 2 (x)∆x + (c1 ) + 4 (c2 ) ∆x4 . dx 24 dx4 dx Solving this equation for the second derivative yields φ(x + ∆x) − 2φ(x) + φ(x − ∆x) d2 φ 1 (x) = − dx2 ∆x2 24



d4 φ d4 φ (c ) + (c ) ∆x2 . 1 2 dx4 dx4

Since

d4 φ d4 φ d4 φ (c1 ) + 4 (c2 ) → 2 4 (x) 4 dx dx dx as ∆x → 0, it follows that the error in the central difference approximation to the second derivative is approximately proportional to ∆x2 . 5. With a Dirichlet condition at the left endpoint and a Neumann condition at the right, we must compute (j)

ui

. = u(xi , tj ), i = 1, 2, . . . , k, j = 1, 2, . . . , n

(using the same notation as in the text). We apply the usual 2-2 scheme at each (xi , tj ), i = 1, 2, . . . , k − 1, (j+1) (j+1) , i = 1, 2, . . . , k − 1. To compute uk , resulting in equation (7.41) (page 301 in the text) for computing ui we approximate the Neumann condition ∂u (xk , tj ) = 0 ∂x by (j) (j) uk+1 − uk−1 (j) (j) = 0 ⇒ uk+1 = uk−1 . ∆x (j+1)

(1)

We then obtain equation (7.48) for uk . Finally, we use (7.44) to compute ui , i = 1, 2, . . . , k − 1, and (7.41) (1) with i = k and ψ(xi+1 ) replaced by ψ(xi−1 ) to compute uk . To test the scheme, we try in on the given problem, using a time interval of [0, 1] with initial values of k = 10 and n = 500 (this gives c∆t/Dt = 0.9). Doubling k and n at each iteration, we obtain the following maximum errors: ∆x 0.1 0.05 0.025 0.0125

∆t 0.002 0.001 0.0005 0.00025

These errors display the expected second-order behavior.

maximum error 1.9 · 10−2 4.8802 · 10−3 1.2250 · 10−3 3.0792 · 10−4

Chapter 8

First-order PDEs and the Method of Characteristics 8.1

The simplest PDE and the method of characteristics

1. Consider the PDE ∂u = 0. ∂y (a) Intuitively, since the solution is constant in y, the characteristics must be the vertical lines x = const. We can this formally as follows. Suppose we impose the initial condition u(f (s), g(s)) = h(s). Then the characteristics are determined by the equations ∂x = 0, x(s, 0) = f (s), ∂t ∂y = 1, y(s, 0) = g(s) ∂t which yield x(s, t) = f (s), y(s, t) = g(s) + t. For each fixed s, these equations parametrize the vertical line x = f (s). (b) The IVP ∂u = 0, u(0, y) = φ(y) ∂y is not well-posed. The general solution of the PDE is u(x, y) = ψ(x). The initial condition u(0, y) = φ(y) then yields ψ(0) = φ(y), which is only possible if φ is a constant: φ(y) = c. Moreover, in this case, there are infinitely many solutions, one for each ψ satisfying ψ(0) = c. Therefore the given IVP either has no solution (if φ is not a constant function) or infinitely many (if φ is a constant function). 3. Consider the PDE ∂u ∂u +3 = 0, 0 < x < 1, t > 0. ∂x ∂t Given the initial condition u(f (s), g(s)) = h(s), the characteristics are determined by 4

∂x = 4, x(0, τ ) = f (s), ∂τ ∂t = 3, t(0, τ ) = g(s), ∂τ which yield x(s, τ ) = 4τ + f (s), t(s, τ ) = 3τ + g(s).

63

64

CHAPTER 8. FIRST-ORDER PDES AND THE METHOD OF CHARACTERISTICS (a) Now consider the initial/boundary conditions u(x, 0) = u0 (x), 0 < x < 1, u(0, t) = φ(t), t > 0. Characteristics ing through (s, 0), 0 < s < 1, have the form x = 4τ + s, t = 3τ ⇔ τ =

4 t , s = x − t. 3 3

This characteristic can also be written as t = (3/4)(x − s). On such characteristics, u(x, t) = v(s, τ ) = u0 (s) = u0

4 x− t . 3

Characteristics through (0, s), s > 0, have the form x = 4τ, t = 3τ + s ⇔ τ =

x 3 , s = t − x. 4 4

This characteristic can also be written as t = s + (3/4)x. On such characteristics, 3 u(x, t) = v(s, τ ) = φ(s) = φ t − x . 4 The characteristic dividing the domain into two is t = (3/4)x. Thus u(x, t) =

u0 x − 43 t , φ t − 34 x ,

t ≤ 34 x, t > 34 x.

(b) Now consider the initial/boundary conditions u(x, 0) = u0 (x), 0 < x < 1, u(1, t) = ψ(t), t > 0. This set of conditions does not define a well-posed problem, since every characteristic through a point (x0 , 0), 0 < x0 < 1, also es through a point (1, t0 ), 0 < t0 < 3/4. These means that the solution is over-determined on part of the domain, and u0 , ψ cannot be specified independently of each other. 5. Consider the IBVP ∂u ∂u +c = 0, 0 < x < 1, t > 0, ∂t ∂x u(x, 0) = u0 (x), 0 < x < 1, u(0, t) = φ(t), t > 0. The solution is u(x, t) =

u0 (x − ct), φ(t − c−1 x),

t < c−1 x, t > c−1 x.

We assume that u0 and φ are both continuously differentiable. (a) We wish to determine the conditions on u0 , φ that guarantee that u is continuous. Consider any (x0 , t0 ) with t0 = c−1 x0 . We have u0 (x − ct) → u0 (x0 − ct0 ) = u0 (0) as (x, t) → (x0 , t0 ), φ(t − c−1 x) → φ(t0 − c−1 x0 ) = φ(0) as (x, t) → (x0 , t0 ), which shows that u is continuous if and only if u0 (0) = φ(0).

8.2. FIRST-ORDER QUASILINEAR PDES

65

(b) Next, we have ∂u (x, t) = ∂x



du0 (x − ct), dx −1 dφ −c dt (t − c−1 x),

t < c−1 x, t > c−1 x,

and du0 du0 du0 (x − ct) → (x0 − ct0 ) = (0) as (x, t) → (x0 , t0 ), dx dx dx dφ dφ dφ −c−1 (t − c−1 x) → −c−1 (t0 − c−1 x0 ) = −c−1 (0) as (x, t) → (x0 , t0 ), dt dt dt Therefore, ∂u/∂x is continuous if and only if du0 dφ (0) = −c−1 (0). dx dt Similar reasoning shows that this condition also ensures that ∂u/∂t is continuous. 7. Consider the IVP ∂u ∂u + = x + y, (x, y) ∈ R2 , ∂y ∂x u(x, −x) = x2 , x ∈ R. We solve the problem in characteristics variables as follows: ∂x = 1, x(s, 0) = s, ∂t ∂y = 1, y(s, 0) = −s, ∂t ∂v = x + y, v(s, 0) = s2 . ∂t The first two equations yield x = s + t, y = −s + t ⇔ s = and therefore

x−y x+y , t= , 2 2

∂v = 2t, v(s, 0) = s2 ⇒ v(s, t) = s2 + t2 . ∂t

Thus, u(x, y) =

8.2

x − y 2 2

+

x + y 2 2

=

1 2 (x + y 2 ). 2

First-order quasilinear PDEs

1. Consider the following IVP: x

∂u ∂u +y − (x + y)u = 0, ∂x ∂y u(1, y) = φ(y).

(a) The characteristics are defined by the equations ∂x = x, x(s, 0) = 1, ∂t ∂y = y, y(s, 0) = s, ∂t which yield

y . x Notice that x must be positive. Also, since s represents y0 on the initial curve x = 1, we see that the characteristics can be written as y = y0 x. Thus the characteristic through (1, y0 ) is the line with slope y0 ing through the origin. This also shows that there must be a discontinuity when x decreases to zero. x = et , y = set ⇔ t = ln x, s =

66

CHAPTER 8. FIRST-ORDER PDES AND THE METHOD OF CHARACTERISTICS (b) In the characteristic variables, the solution v(s, t) is defined by ∂v = (x + y)v = (1 + s)et v, v(s, 0) = φ(s). ∂t The ODE is separable and therefore the IVP is easily solved: v(s, t) = φ(s)e(1+s)(e Changing variables, we obtain

t

−1)

.

y

ex+y−1−y/x . x (c) The formula for u shows that it is defined for all x > 0 (or for all x < 0, but the initial curve lies in the right half-plane). This is consistent with the analysis of the characteristics. u(x, y) = φ

3. Consider the IVP u

∂u ∂u +u = 0, ∂x ∂y u(x, 0) = φ(x).

The PDE is quasi-linear and can be solved by solving the characteristic equations: ∂x = v, x(s, 0) = s, ∂t ∂y = v, y(s, 0) = 0, ∂t ∂v = 0, v(s, 0) = φ(s). ∂t The equations for v yield v(s, t) = φ(s); substituting this into the first two equations yields x(s, t) = s + φ(s)t, y = φ(s)t. From this we obtain s = x − y, t = y/φ(x − y), and the solution is u(x, y) = v(s, t) = φ(s) = φ(x − y), that is, u(x, y) = φ(x − y). If we test the initial curve using the Jacobian test, we obtain ∂x (s, 0) ∂x (s, 0) 1 ∂s ∂t ∂y = (s, 0) ∂y (s, 0) 0 ∂s ∂t

φ(s) = φ(s), φ(s)

we see that the initial curve is noncharacteristic if and only if φ(s) 6= 0. 5. Consider the IVP ∂u ∂u +u = 0, x > 0, ∂x ∂y u(x, 0) = x, x > 0. (a) The IVP can be solved by solving the characteristic equations: ∂x = 1, x(s, 0) = s, ∂t ∂y = v, y(s, 0) = 0, ∂t ∂v = 0, v(s, 0) = s. ∂t It is straightforward to obtain the solution: x(s, t) = s + t, y(s, t) = st, v(s, t) = s. Solving x = s + t, y = st for s, t (and bearing in mind that s = x when t = 0), we obtain p p x + x2 − 4y x − x2 − 4y , t= . s= 2 2

8.3. BURGERS’S EQUATION

67

Since v(s, t) = s, it follows that u(x, y) =

x+

p

x2 − 4y . 2

This solution is defined for x > 0, y ≤ x2 /4. (b) The characteristic indexed by s es through the point (s, 0), so we can regards s as x0 . Notice that x = s + t, y = st implies that t = x − s and hence y = s(x − s). This shows that the characteristics are the lines y = x0 (x − x0 ). (c) Since the slopes of the characteristics increase as x0 increases, it follows that, for x1 > x0 , the characteristics through (x0 , 0) and (x1 , 0) will intersect at some point (x, y) with x > x1 . A simple calculation shows that y = x0 (x − x0 ) and y = x1 (x − x1 ) intersect at the point (x0 + x1 , x0 x1 ). Taking the limit as x1 decreases toward x0 reveals the first point where the characteristic through (x0 , 0) intersects another characteristic: (2x0 , x20 ). Notice that this is a point on the curve y = x2 /4, which is consistent with the results obtained above.

8.3

Burgers’s equation

1. Consider the inviscid Burgers’s equation ∂u ∂u +u =0 ∂t ∂x with initial condition u(x, 0) = u0 (x) =

1 . 1 + x2

We have seen that the characteristics are the lines x = x0 + u0 (x0 )t ⇔ t = u0 (x0 )−1 (x − x0 ). (a) We wish to show that the solution of the IVP is well-defined for at each point (x, t), x < 0, t > 0. To show this, it suffices to prove that there is a unique characteristic ing through any such point. Since the slope (in the x, t plane) of every characteristic is positive, the characteristic through (x0 , 0) es through (x, t) only if x0 < x < 0 and t = (1 + x20 )(x − x0 ) ⇔ x30 − xx20 + x0 + (t − x) = 0. If we define ψ(x0 ) = x30 − xx20 + x0 + (t − x), then it is easy to that ψ(0) > 0, ψ(x0 ) < 0 for x0 sufficiently large and negative, and ψ 0 (x0 ) > 0 for all x0 < 0. It follows that ψ(x0 ) = 0 has a unique solution x0 < 0, and hence there is a unique characteristic ing through (x, t). (b) We now wish to describe the set of all (x, t), t > 0, for which the solution is uniquely defined. Referring to the analysis in the text, culminating in Figure 8.9, we see that u(x, t) is uniquely defined for all (x, t) outside the caustic, which is defined by the parametric equations x=

(1 + x20 )2 1 + 3x20 , t= , x0 > 0. 2x0 2x0

The corner of the caustic can be found by setting the derivatives of x and t (with respect to x0 ) to zero and solving, which yields √ 1 8 x0 = √ , x = 3, t = √ . 3 3 3 Moreover, we can solve 1 + 3x20 x= 2x0 for x0 in of x to get √ x ± x2 − 3 x0 = . 3 Substituting into the equation for t yields the two branches of the caustic: √ !2 √ !2 3 1+ t1 =

x−

x2 −3 3

√ 2(x − x2 − 3)

2

3 1+ , t2 =

x+

2(x +

√

x2 −3 3

x2

− 3)

2

, x≥

√

3

68

CHAPTER 8. FIRST-ORDER PDES AND THE METHOD OF CHARACTERISTICS (where t2 > t1 for all x). Comparing to Figure 8.9, we see that u(x, t) is uniquely defined if t > 0 and √ √ x < 3 or x ≥ 3 and (t < t1 or t > t2 . 3. Consider the IVP ∂u ∂u + (1 − u) = 0, −∞ < x < 0, t > 0, ∂t ∂x 1 u(x, 0) = , −∞ < x < 0. x (a) The characteristic equations are ∂x = 1 − v, x(s, 0) = s, ∂τ ∂t = 1, t(s, 0) = 0, ∂τ ∂v 1 = 0, v(s, 0) = . ∂τ s Solving yields x(s, t) = s +

s−1 1 τ, t(s, t) = τ, v(s, t) = . s s

Setting s = x0 , we obtain

x0 x0 − 1 t ⇔ t= (x − x0 ). x0 x0 − 1 These are the equations of the characteristics. x = x0 +

(b) We can solve x=s+

s−1 t s

for s as follows: x=s+

s−1 t ⇒ s2 + (t − x)s − t = 0 s p x − t ± (t − x)2 + 4t ⇒s= 2 p x − t − (t − x)2 + 4t ⇒s= , 2

where the last step uses the fact that s must be negative. We then obtain u(x, t) = v(s, τ ) = 1/s, that is, u(x, t) =

2 p . x − t − (t − x)2 + 4t

(c) Notice that the slope of the characteristic decreases as x0 increases; therefore, the characteristics through (x0 , 0) and (x1 , 0) will intersect at some point (x, t) with t < 0. Solving the equations t=

x0 x1 (x − x0 ), t = (x − x1 ) x0 − 1 x1 − 1

for (x, t) yields (x, t) = (x0 − x0 x1 + x1 , −x0 x1 ). Notice that t < 0, as expected. 5. Consider the IVP ∂u ∂u + a(u) = 0, −∞ < x < ∞, t > 0, ∂t ∂x u(x, 0) = φ(x), x > 0. The characteristics are defined by x(s, τ ) = s + a(φ(s))τ, t(s, τ ) = τ, v(s, τ ) = φ(s).

8.3. BURGERS’S EQUATION

69

We can write the characteristic ing through (x0 , 0) as x = x0 + a(φ(x0 ))t ⇔ t =

x − x0 . a(φ(x0 ))

Assuming that a(φ(x)) is an decreasing function of x, we know that characteristic ing through (x0 , 0) will intersect with the characteristic ing through (x1 , 0) if x1 > x0 . We can compute the point of intersection as follows: x = x0 + a(φ(x0 ))t, x = x1 + a(φ(x1 ))t ⇒ x0 + a(φ(x0 ))t = x1 + a(φ(x1 ))t x0 − x1 1 ⇒t= = − a(φ(x ))−a(φ(x )) . 1 0 a(φ(x1 )) − a(φ(x0 )) x1 −x0

If we take the limit as x1 decreases to x0 , we find the smallest value of t for the characteristic through (x0 , 0) intersects another characteristic: t∗ = −

1

d (a(φ(x))) x=x dx 0

1 (φ(x )) dφ (x0 ) 0 du dx

= − da

On the other hand, suppose we look for the smallest value of τ for which the change of variables from (s, τ ) to (x, t) is no longer well-defined. This occurs when the following Jacobian determinant is first zero: ∂x da dφ da (s, τ ) ∂x (s, τ ) 1 + du (φ(s)) dφ (s)τ a(φ(s)) ∂τ dx ∂s = ∂t ∂t = 1 + du (φ(s)) dx (s)τ. (s, τ ) (s, τ ) 0 1 ∂s

∂τ

Setting this expression equal to zero, solving for τ , and replacing s by x0 and τ by t yields the same result as before.

Chapter 9

Green’s Functions 9.1

Green’s functions for BVPs in ODEs: Special cases

1. u(x) = xe−x . 3.

(a) This BVP models a bar whose top end (originally at x = 0) is free and whose bottom end is fixed at x = `. If we apply a unit force to the cross-section at x = s, then the part of the bar originally between x = s and x = ` will compress, and the part of the bar originally above x = s will just move rigidly with u(x) = u(s) for 0 ≤ x < s. The compression of the bottom part of the bar will satisfy Hooke’s law: k

u(x) `−x = 1 ⇒ u(x) = `−x k

(the uncompressed length of the part of the bar between x = s and x = ` is ` − x). Therefore we obtain ( `−s , 0 < x < s, k u(x) = `−x , s < x < `. k The Green’s function is

(

`−s , k `−x , k

G(x; s) =

0 < x < s, s < x < `.

(b) The Green’s function is Z

`

G(x; s) = max{x,s}

dz . k(z)

(c) If k is constant, then ` − max{x, s} G(x; s) = = k

(

`−s , k `−x , k

0 < x < s, s < x < `.

(d) u(x) =

ln (1 + x) ln 2 1 x2 x − − + − . 2 4 4 2 2

5. By direct integration, we obtain the following solution to (9.11): Z x ds u(x) = p . 0 k(s) Since x = min{x, `} for x ∈ [0, `], we can write min{x,s}

Z u(x) = p 0

ds , k(s)

that is, u(x) = g(x; `)p.

71

72

CHAPTER 9. GREEN’S FUNCTIONS 2 2 7. Let the operators K : Cm [0, `] → C[0, `] and M : C[0, `] → Cm [0, `] be defined by

d dx

Ku = −

du k(x) dx

and `

Z (M f )(x) =

G(x; y)f (y) dy, 0

where min{x,y}

Z G(x; y) = 0

1 dz. k(z)

2 We wish to prove directly that (M K)u = u for all u ∈ Cm [0, `]. Given the definition of G, we will write x

Z (M f )(x) =

Z

`

G(x; y)f (y) dy + 0

0

1 f (y) dz dy + k(z)

Z

x

1 f (y) dy dz + k(z)

Z

0 x

Z

= 0

`

y

= Z

G(x; y)f (y) dy x

xZ

Z

z

x

y

1 f (y) dz dy k(z)

`

1 f (y) dy dz, k(z)

0 x

0

Z Z

x

where the last step follows from changing the order of integration. Therefore, Z xZ ` du du 1 d 1 d k(y) (y) dy dz − k(y) (y) dy dz dx dx 0 0 z k(z) dx x k(z) dx Z x Z x 1 du du 1 du du =− k(x) (x) − k(z) (z) dz − k(`) (`) − k(x) (x) dz dx dx dx dx 0 k(z) 0 k(z) Z x du 1 −k(z) (z) dz =− dx 0 k(z) Z x du = (z) dz = u(x) 0 dx Z

x

Z

x

(M Ku)(x) = −

(notice how we used both boundary conditions in the course of this calculation). This is the desired result.

9.2

Green’s functions for BVPs in ODEs: the symmetric case

1. The Green’s function is ( G(x; y) =

sin (θy)((tan θ) cos (θx)−sin (θx)) , θ tan θ sin (θx)((tan θ) cos (θy)−sin (θy)) , θ tan θ

0 ≤ y ≤ x, x ≤ y ≤ 1.

3. The Green’s function is G(x; y) =

−kv1 (y)v2 (x), −kv1 (x)v2 (y),

0 ≤ y ≤ x, x ≤ y ≤ 1,

where v1 (x) = eθx − e−θx , u2 (x) = eθx − eθ(2−x) , and k=

1 . 2(e2θ − 1)θ

5. Let G = G(x, y) be any continuous function of two variables, a ≤ x ≤ b, a ≤ y ≤ b, and define M : C[a, b] → C[a, b] by Z b (M f )(x) = G(x, y)f (y) dy. a

9.2. GREEN’S FUNCTIONS FOR BVPS IN ODES: THE SYMMETRIC CASE

73

Suppose G is symmetric in the sense that G(y, x) = G(x, y) for all x, y ∈ [a, b]. We then have Z b (M f, g)L2 = (M f )(x)g(x) dx a

Z

b

Z

b

Z

b

= a

Z

G(x, y)f (y) dy g(x) dx

a b

=

G(x, y)f (y)g(x) dy dx a

Z

a b

b

Z

=

G(x, y)f (y)g(x) dx dy (changing the order of integration) a

Z

a b

=

Z

b

G(x, y)g(x) dx dy

b

G(y, x)g(x) dx dy (since G is symmetric)

f (y) a

a

Z

b

=

Z f (y)

a

Z =

a b

f (y)(M g)(y) dy = (f, M g)L2 . a

Thus M is a symmetric operator with respect to the L2 (a, b) inner product. 7. Define H : R → R by H(x) = We wish to prove that Z

0, 1.

∞

x < 0, 0<x

Z

∞

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

H(x) −∞

dφ (x) dx dx

holds for all φ ∈ D(R). That is, we wish to prove that if φ ∈ D(R), then Z ∞ dφ H(x) (x) dx = −φ(0). dx −∞ This is a direct calculation. Choose b > 0 sufficiently large that φ(x) = 0 for all x ≥ b. Then Z ∞ Z ∞ Z b dφ dφ dφ H(x) (x) dx = (x) dx = (x) dx = φ(b) − φ(0) = −φ(0), dx dx −∞ 0 0 dx as desired. 9. Let G be the Green’s function for the BVP du d − P (x) + R(x)u = f (x), a < x < b, dx dx du α1 u(a) + α2 (a) = 0, dx du β1 u(b) + β2 (b) = 0. dx Then

( G(x; y) =

(y)v2 (x) − vP1(x)W , a ≤ y ≤ x, (x) (x)v2 (y) − vP1(x)W , x ≤ y ≤ b, (x)

where v1 and v2 are certain solutions of the homogeneous ODE d du − P (x) + R(x)u = 0, a < x < b. dx dx We wish to show that, for fixed y ∈ (a, b), u(x) = G(x; y) is the solution of the above BVP with the right-hand-side function f (x) replaced by δ(x − y). To do this, we must interpret the BVP in its weak sense, which is b Z b Z b du du dv −P (x) (x)v(x) + P (x) (x) (x) + R(x)u(x)v(x) dx = δ(x − y)v(x) dx = v(y). dx dx dx a a a

74

CHAPTER 9. GREEN’S FUNCTIONS Here u must satisfy the given boundary conditions and the above variational equation must hold for all test functions v satisfying the same boundary conditions. Now, the Green’s function is continuous but its derivative du/dx has a jump discontinuity at x = y. We will write v 0 (y)v2 (y) du du (y−) = lim (x) = − 1 , dx P (y)W (y) x→y − dx v1 (y)v20 (y) du du (y+) = lim (x) = − . dx P (y)W (y) x→y + dx We now substitute u into the left-hand side of the variational equation to obtain b Z b du dv du (x)v(x) + P (x) (x) (x) + R(x)u(x)v(x) dx dx dx dx a a Z y du du dv du = P (a) (a)v(a) − P (b) (b)v(b) + P (x) (x) (x) + R(x)u(x)v(x) dx dx dx dx dx a Z b dv du + P (x) (x) (x) + R(x)u(x)v(x) dx. dx dx y −P (x)

The integrands of both integrals on smooth on the intervals of integration, and hence integration by parts is valid: b Z b du du dv (x)v(x) + P (x) (x) (x) + R(x)u(x)v(x) dx dx dx dx a a du du du du du du = P (a) (a)v(a) − P (b) (b)v(b) + P (y) (y − )v(y) − P (a) (a)v(a) + P (b) (b)v(b) − P (y) (y + )v(y) dx dx dx dx dx dx Z b Z y du du d d P (x) (x) v(x) + R(x)u(x)v(x) dx + P (x) (x) v(x) + R(x)u(x)v(x) dx + − − dx dx dx dx y a Z b du du du d = P (y) (y − )v(y) − P (y) (y + )v(y) + P (x) (x) + R(x)u(x) v(x) dx − dx dx dx dx a du − du + = P (y) (y )v(y) − P (y) (y )v(y), dx dx −P (x)

where the last step follows because, on both (a, y) and (y, b), u(x) = G(x; y) is a smooth solution of the homogeneous ODE d du − P (x) + R(x)u = 0. dx dx We now have b Z b du du dv −P (x) (x)v(x) + P (x) (x) (x) + R(x)u(x)v(x) dx dx dx dx a a du − du + = P (y)v(y) (y ) − (y ) dx dx ! 2 1 (y) − dv (y)v2 (y) v1 (y) dv dx dx = P (y)v(y) = v(y) P (y)W (y) 2 since W (y) = v1 (y) dv (y) − dx

dv1 (y)v2 (y). dx

This completes the proof.

11. We wish to find a formula for the solution of d du − P (x) + R(x)u = 0, a < x < b, dx dx du α1 u(a) + α2 (a) = va , dx du β1 u(b) + β2 (b) = vb , dx

9.3. GREEN’S FUNCTIONS FOR BVPS IN ODES: THE GENERAL CASE

75

in of the Green’s function G for the BVP d du − P (x) + R(x)u = f (x), a < x < b, dx dx du α1 u(a) + α2 (a) = 0, dx du β1 u(b) + β2 (b) = 0. dx We proceed as in Section 9.2.2; in particular, we apply dv du dv du (u, Lv) − (Lu, v) = P (a) u(a) (a) − (a)v(a) − P (b) u(b) (b) − (b)v(b) dx dx dx dx as in the text (where L is defined by (9.16)), with v(x) = G(x; y) (for a fixed y) and u the solution of the above BVP with inhomogeneous boundary conditions. We then have (Lv)(x) = δ(x − y), and hence (u, Lv) = u(y). Also, Lu = 0 since u satisfies the homogeneous ODE, and therefore (Lu, v) = 0. We thus obtain dv du dv du u(y) = P (a) u(a) (a) − (a)v(a) − P (b) u(b) (b) − (b)v(b) . dx dx dx dx We now apply the boundary condition satisfied by u and v: du (a) = va ⇒ dx du β1 u(b) + β2 (b) = vb ⇒ dx dv α1 v(a) + α2 (a) = 0 ⇒ dx dv β1 v(b) + β2 (b) = 0 ⇒ dx

α1 u(a) + α2

du 1 α1 (a) = va − u(a), dx α2 α2 du 1 β1 (b) = vb − u(b), dx β2 β2 dv α1 (a) = − v(a), dx α2 dv β1 (b) = − v(b). dx β2

We thus obtain du dv du dv (a)v(a) − P (b) u(b) (b) − (b)v(b) u(y) = P (a) u(a) (a) − dx dx dx dx α1 1 α1 β1 1 β1 = P (a) − u(a)v(a) − va v(a) + u(a)v(a) − P (b) − u(b)v(b) − vb v(b) + u(b)v(b) α2 α2 α2 β2 β2 β2

= β2−1 P (b)v(b)vb − α2−1 P (a)v(a)va = β2−1 P (b)G(b; y)vb − α2−1 P (a)G(a; y)va = β2−1 P (b)G(y; b)vb − α2−1 P (a)G(y; a)va . In the last step we used the symmetry of G. This is the desired formula: u(y) = β2−1 P (b)G(y; b)vb − α2−1 P (a)G(y; a)va .

9.3

Green’s functions for BVPs in ODEs: the general case

1. The Green’s function is

( G(x; y) =

(y)v2 (x) , − vP1(x)W (x) v1 (x)v2 (y) − P (x)W (x) ,

0 ≤ y ≤ x, x ≤ y ≤ 1,

where v1 (x) = ex − e2x , v2 (x) = ex − e2x−1 , and P (x)W (x) is the constant 1 − e−1 . The solution to the BVP is 1

Z u(x) =

G(x; y)f (y)w(y) dy, 0

where w(x) = e−3x .

76

CHAPTER 9. GREEN’S FUNCTIONS 3. The Green’s function is

( G(x; y) =

(y)v2 (x) , − vP1(x)W (x) v1 (x)v2 (y) − P (x)W (x) ,

0 ≤ y ≤ x, x ≤ y ≤ 1,

where v1 (x) = ex sin (2x), v2 (x) = ex (sin(2x) − (tan (2)) cos (2x)), and P (x)W (x) is the constant 2 tan (2). The solution to the BVP is Z 1

G(x; y)f (y)w(y) dy,

u(x) = 0

where w(x) = e−2x . 5. The Green’s function is

G(x; y) =

ey 6 + 6x + 3x2 + x3 , x 2 3 e 6 + 6y + 3y + y ,

1 ≤ y ≤ x, 2 x ≤ y ≤ 1.

The solution to the BVP when f (x) = x is u(x) =

23 23 + x + 6x2 + 2x3 − 7ex−1 . 2 2

Rb 7. Let G be the Green’s function derived in this section, and let M be defined by (M f )(x) = a G(x; y)f (y) dy. (Here 2 M maps complex C[a, b] into complex C [a, b].) We wish to prove that (M f, g)w = (f, M g)w for all f, g ∈ C[a, b]. In the following calculation, we use the facts that G is symmetric (G(x; y) = G(y; x) for all x, y ∈ [a, b]) and G is real-valued, so that G(x; y) = G(x; y). We have Z b Z b Z b (M f, g)w = (M f )(x)g(x)w(x) dx = G(x; y)f (y)w(y) dy g(x)w(x) dx a

a

b

Z

a b

Z

=

G(x; y)f (y)w(y)g(x)w(x) dy dx a

a bZ

Z

b

=

G(x; y)f (y)w(y)g(x)w(x) dx dy a

a b

Z =

b

G(x; y)g(x)w(x) dx w(y) dy

b

G(x; y)g(x)w(x) dx w(y) dy

b

G(y; x)g(x)w(x) dx w(y) dy

Z f (y)

a

a b

Z =

Z f (y)

a

a b

Z =

Z f (y)

a

a b

Z

f (y)(M g)(y)w(y) dy = (f, M g)w .

= a

(Notice how the key step was the interchanging of the order of integration.) This proves the desired result.

9.4

Introduction to Green’s functions for IVPs

1. We wish to show that the solution of d2 u du + θ2 u = f (t), u(0) = u0 , (0) = v0 dt2 dt is u(t) =

∂G (t; 0)u0 + G(t; 0)v0 + ∂t

where G(t; s) =

∞

Z

sin (θ(t−s)) , θ

0,

G(t; s)f (s) ds, 0

t > s, t < s.

We can write the proposed solution explicitly as u(t) = cos (θt)u0 +

sin (θt) v0 + θ

Z 0

t

sin (θ(t − s)) f (s) ds. θ

9.4. INTRODUCTION TO GREEN’S FUNCTIONS FOR IVPS

77

Then Z t sin (θ(t − t)) du (t) = −θ sin (θt)u0 + cos (θt)v0 + f (t) + cos (θ(t − s))f (s) ds dt θ 0 Z t cos (θ(t − s))f (s) ds = −θ sin (θt)u0 + cos (θt)v0 + 0

Z t d2 u 2 (t) = −θ cos (θt)u − θ sin (θt)v + cos (θ(t − t))f (t) − θ sin (θ(t − s))f (s) ds 0 0 dt2 0 Z t θ sin (θ(t − s))f (s) ds. = −θ2 cos (θt)u0 − θ sin (θt)v0 + f (t) − 0

From these formulas, we can easily that u satisfies the ODE and the initial conditions: Z 0 sin (0) sin (θ(t − s)) u(0) = cos (0)u0 + v0 + f (s) ds = u0 , θ θ 0 Z 0 du (0) = −θ sin (0)u0 + cos (0)v0 + cos (θ(t − s))f (s) ds = v0 , dt 0 and Z t d2 u 2 2 (t) + θ u(t) = −θ cos (θt)u − θ sin (θt)v + f (t) − θ sin (θ(t − s))f (s) ds+ 0 0 dt2 0 Z t θ2 cos (θt)u0 + θ sin (θt)v0 + θ sin (θ(t − s))f (s) ds 0

= f (t). 3. The Green’s function is G(t; s) = 5. The Green’s function is

G(t; s) =

sin (2(t−s)) , 2

0,

t > s, t < s.

e−(t−s) − e−2(t−s) , 0,

7. The Green’s function is

G(t; s) =

e0.02(t−s) , 0,

t > s, t < s.

t > s, t < s,

and the solution to the IVP is Z ∞ P (t) = G(t; 0)P0 + G(t; s)f (s) ds = 55.5e0.02t + 450 + 10t − 450e0.02t = 450 + 10t − 394.5e0.02t . 0

9. Let the solution of du d2 u + a1 (t) + a0 (t)u = 0, dt2 dt u(s) = 0, du (s) = f (s) dt be u(t) = S(t; s)f (s), and define the Green’s function G(t; s) by S(t; s), t > s, G(t; s) = 0, t < s, We wish to show that G satisfies ∂2G ∂G (t; s) + a1 (t) (t; s) + a0 (t)G(t; s) = δ(t − s), ∂t2 ∂t G(t0 ; s) = 0, ∂G (t0 ; s) = 0, ∂t

78

CHAPTER 9. GREEN’S FUNCTIONS provided t0 < s. We notice that G(t, s) = H(t − s)S(t, s), where H is the Heaviside function. We will show that ∂S ∂2G ∂G ∂2S (t, s) = H(t − s) (t, s) and (t, s) = δ(t − s) + H(t − s) 2 (t, s) 2 ∂t ∂t ∂t ∂t (in the weak sense). Let φ be any smooth function having compact in [t0 , ∞), that is, any function for which φ(t) = 0 for all t 6∈ (a, b), where t0 < a < b < ∞. We first show that Z ∞ Z ∞ ∂S dφ H(t − s) G(t; s) (t) dt = (t, s)φ(t) dt, − dt ∂t t0 t0 which implies that ∂G ∂S (t, s) = H(t − s) (t, s) ∂t ∂t in the weak sense. We have Z ∞ Z ∞ dφ dφ − G(t; s) (t) dt = − H(t − s)S(t, s) (t) dt dt dt t0 t0 Z b dφ =− S(t, s) (t) dt dt s Z b ∂S = − S(t, s)φ(t)|bs + (t, s)φ(t) dt (by integration by parts) s ∂t Z b ∂S (t, s)φ(t) dt (since S(s, s) = 0 and φ(b) = 0) = ∂t Zs ∞ ∂S = H(t − s) (t, s)φ(t) dt. ∂t t0 This proves the desired formula for ∂G/∂t. Next, we will show that Z ∞ Z ∞ ∂G dφ ∂2S − (t; s) (t) dt = φ(s) + H(t − s) 2 (t, s)φ(t) dt, ∂t dt ∂t t0 t0 which implies that ∂2G ∂2S (t, s) = δ(t − s) + H(t − s) 2 (t, s) 2 ∂t ∂t in the weak sense. We have Z ∞ Z ∞ ∂G dφ ∂S dφ − (t; s) (t) dt = − H(t − s) (t, s) (t) dt ∂t dt ∂t dt t0 t0 Z b ∂S dφ =− (t, s) (t) dt ∂t dt s b Z b 2 ∂S ∂ S =− (t, s)φ(t) + (t, s)φ(t) dt (by integration by parts) ∂t ∂t2 s s Z b 2 ∂ S = φ(s) + (t, s)φ(t) dt (since ∂S (s, s) = 1 and φ(b) = 0) ∂t 2 s ∂t Z ∞ ∂2S = φ(s) + H(t − s) 2 (t, s)φ(t) dt. ∂t t0 This proves the desired formula for ∂G2 /∂t2 . Since S(t, s), as a function of t, satisfies the homogeneous version of the ODE, we see that ∂2G ∂G (t; s) + a1 (t) (t; s) + a0 (t)G(t; s) = δ(t − s), ∂t2 ∂t as desired. Also, ∂G ∂S (t0 , s) = H(t0 − s) (t0 , s) = 0 ∂t ∂t since H(t0 − s) = 0 because t0 − s < 0 by assumption. This completes the proof. G(t0 , s) = H(t0 − s)S(t0 , s) = 0,

9.5. GREEN’S FUNCTIONS FOR THE HEAT EQUATION

9.5

79

Green’s functions for the heat equation

1. Given that

2 1 e−x /(4kt) , S(x, t) = √ 2 πkt

we have

∞

Z ∞ 2 1 S(x, t) dx = √ e−x /(4kt) dx. 2 πkt −∞ −∞ √ Making the change of variables u = x/(2 kt), we obtain Z ∞ Z ∞ 2 1 1 √ S(x, t) dx = √ e−u du = √ π = 1. π −∞ π −∞ Z

Also, for any x, y ∈ R,

∞

Z

Z

∞

S(x − y, t) dx = −∞

S(v, t) dv = 1. −∞

3. For t = 660, six of the Fourier series are sufficient to give an accurate graph. For t = 61, even 100 are not sufficient to give a correct graph. 5. Following the example in Section 9.5.2, we obtain the following Green’s function: ( P∞ −k(2n−1)2 π2 (t−s)/(4`2 ) 2 sin (2n−1)πx , e sin (2n−1)πy n=1 `ρc 2` 2` G(x, t; y, s) = 0,

0 ≤ s ≤ t, s > t.

Here we used the correct eigenfunctions for the given mixed boundary conditions and also adjust for the the fact that the constants appear differently in the PDE in this exercise than in the example in Section 9.5.2 (notice the extra factor of 1/(ρc) in the formula for this Green’s function). As before, k = κ/(ρc). Snapshots of G(x, t; 75, 60) are shown in Figure 9.1. Graph of G(x,t;75,60) for various values of t

temperature

0.02 t=120 t=240 t=360 t=480 t=600

0.015

0.01

0.005

0 0

20

40

60

80

100

x

Figure 9.1: Snapshots of the Green’s function in Exercise 9.5.5 at times t = 120, 240, 360, 480, 600 seconds. Twenty of the Fourier series were used to create these graphs. 7. We have `

Z u(x, t) =

G(x, t; y, 0)φ(y) dy 0

∞ nπy nπx 2 X −kn2 π2 t/`2 e sin sin φ(y) dy ` ` 0 ` n=1 ∞ Z ` nπy nπx X 2 2 2 2 = sin φ(y) dy e−kn π t/` sin ` 0 ` ` n=1

Z

`

=

=

∞ X n=1

2

bn e−kn

π 2 t/`2

sin

nπx `

,

80

CHAPTER 9. GREEN’S FUNCTIONS where 2 `

bn =

Z

`

sin 0

nπy φ(y) dy, n = 1, 2, 3, . . . `

are the Fourier sine coefficients of the function φ. We recognize u as the solution of the IBVP ∂u ∂2u − k 2 = 0, 0 < x < `, t > 0, ∂t ∂x u(x, 0) = φ(x), 0 < x < `, u(0, t) = 0, t > 0, u(`, t) = 0, t > 0 (see the derivation at the beginning of Section 6.1 in the text).

9.6

Green’s functions for the wave equation

1. Let 1 2c

u(x, t) =

t

Z

Z

x+c(t−s)

f (y, s) dy ds. 0

x−c(t−s)

We wish to show that u satisfies the inhomogeneous wave equation (with right-hand side f (x, t)) on the real line, subject to zero boundary conditions. We have 1 ∂u (x, t) = ∂t 2c

Z

1 2c

Z

= =

1 2c

x+c(t−t)

f (y, t) dy + x−c(t−t) x

f (y, t) dy + x Z t

1 2c

1 2c

t

Z

{cf (x + c(t − s), s) + cf (x − c(t − s), s)} ds 0

t

Z

{cf (x + c(t − s), s) + cf (x − c(t − s), s)} ds 0

{cf (x + c(t − s), s) + cf (x − c(t − s), s)} ds, 0

Z t ∂2u 1 1 2 df 2 df (cf (x + c(t − t), t) + cf (x − c(t − t), t)) + c (x + c(t − s), s) − c (x − c(t − s), s) ds (x, t) = ∂t2 2c 2c 0 dx dx Z t 1 df 1 df = (cf (x, t) + cf (x, t)) + c2 (x + c(t − s), s) − c2 (x − c(t − s), s) ds 2c 2c 0 dx dx Z c t df df = f (x, t) + (x + c(t − s), s) − (x − c(t − s), s) ds, 2 0 dx dx Z t ∂u 1 (x, t) = {f (x + c(t − s), s) − f (x − c(t − s), s)} ds, ∂x 2c 0 Z t 1 ∂2u df df (x, t) = (x + c(t − s), s) − (x − c(t − s), s) ds. ∂x2 2c 0 dx dx Thus ∂2u ∂2u c (x, t) − c2 2 (x, t) = f (x, t) + 2 ∂t ∂x 2

t

df df (x + c(t − s), s) − (x − c(t − s), s) ds dx dx 0 Z t c df df − (x + c(t − s), s) − (x − c(t − s), s) ds 2 0 dx dx Z



= f (x, t), as desired. It is obvious from the formulas for u and ∂u/∂t that both are identically zero when t = 0. 3. We wish to show that H(c(t − s) − |x − y|) = (H((x − y) + c(t − s)) − H((x − y) − c(t − s))) H(t − s)

9.6. GREEN’S FUNCTIONS FOR THE WAVE EQUATION

81

for all x, y, t, s ∈ R such that c(t − s) 6= |x − y|. (The two expressions differ when c(t − s) = |x − y| > 0, at least when H is defined as in the text: H(t) = 1 if t > 0 and H(t) = 0 if t ≤ 0.) We have H(c(t − s) − |x − y|) = 1 ⇔ c(t − s) − |x − y| > 0 ⇔ c(t − s) > |x − y| ⇔ −c(t − s) < |x − y| < c(t − s) and t − s > 0 ⇔ x − y + c(t − s) > 0 and x − y − c(t − s) < 0 and t − s > 0 ⇒ H(x − y + c(t − s)) = 1 and H(x − y − c(t − s)) = 0 and H(t − s) = 1 ⇔ (H((x − y) + c(t − s)) − H((x − y) − c(t − s))) H(t − s) = 1. The reasoning can be reversed if we add the condition c(t − s) 6= |x − y|. This completes the proof. 5. Let u be the solution to ∂2u ∂2u − c2 2 = f˜(x, t), −∞ < x < ∞, t > 0, 2 ∂t ∂x u(x, 0) = 0, −∞ < x < ∞, ∂u (x, 0) = 0, −∞ < x < ∞, ∂t where f˜ is odd (f˜(−x, t) = −f˜(x, t)). We wish to show that u is also odd. We will do this by defining v(x, t) = −u(−x, t) and proving that v satisfies the same IVP; then, by uniqueness, it will follow that v = u, that is, u(x, t) = −u(−x, t), as desired. We have ∂v ∂u (x, t) = − (−x, t), ∂t ∂t ∂2v ∂2u (x, t) = − 2 (−x, t), ∂t2 ∂t ∂v ∂u (x, t) = (−x, t), ∂x ∂x 2 ∂ v ∂2u (x, t) = − 2 (−x, t). 2 ∂x ∂x Therefore, ∂2v ∂2v ∂2u ∂2u (x, t) − c2 2 (x, t) = − 2 (−x, t) + c2 2 (−x, t) 2 ∂t ∂x ∂t ∂x 2 2 ∂ u 2∂ u =− (−x, t) − c (−x, t) ∂t2 ∂x2 = −f˜(−x, t) = f˜(x, t) (since f˜ is odd). Thus v satisfies the same wave equation as does u. Moreover, it is obvious that v also satisfies the same initial conditions; hence v = u and the proof is complete. 7. We outline the derivation of the Green’s function for ∂2u ∂2u − c2 2 = f (x, t), 0 < x < `, t > 0, 2 ∂t ∂x u(x, 0) = 0, 0 < x < `, ∂u (x, 0) = 0, 0 < x < `, ∂t ∂u (0, t) = 0, t > 0, ∂x ∂u (`, t) = 0, t > 0. ∂x

82

CHAPTER 9. GREEN’S FUNCTIONS We define f˜ to be the even, periodic extension of f , and define u ˜ to be the solution of ∂2u ˜ ∂2u ˜ − c2 2 = f˜(x, t), −∞ < x < ∞, t > 0, 2 ∂t ∂x u ˜(x, 0) = 0, −∞ < x < ∞, ∂u ˜ (x, 0) = 0, −∞ < x < ∞. ∂t From earlier results, we know that Z t Z x+c(t−s) Z tZ ∞ 1 1 H(c(t − s) − |x − y|)f˜(y, s) dy ds = f˜(y, s) dy ds. u ˜(x, t) = 2c 0 x−c(t−s) 0 −∞ 2c Then

t

∂u ˜ (x, t) = ∂x

Z

∂u ˜ (x, 0) = ∂x

Z

∂u ˜ (x, `) = ∂x

Z

n o f˜(x + c(t − s), s) − f˜(x − c(t − s), s) ds.

0

It follows that t

n o f˜(c(t − s), s) − f˜(−c(t − s), s) ds,

0 t

n o f˜(` + c(t − s), s) − f˜(` − c(t − s), s) ds.

0

Since f˜ is even, we have f˜(c(t − s), s) = f˜(−c(t − s), s), and hence ∂u ˜ (x, 0) = 0. ∂x Also, f˜(` − c(t − s), s) = f˜(−` + c(t − s), s) (since f˜ is even) = f˜(` + c(t − s), s) (since f˜ is 2`-periodic). Therefore, ∂u ˜ (x, `) = 0. ∂x It follows that u ˜, restricted to the interval 0 < x < `, solves the original IBVP. The derivation of the Green’s function for the IBVP is almost exactly the same as that given in the text (for Dirichlet conditions). The result is G(x, t; y, s) =

∞ X n=−∞

1 {H(c(t − s) − |x − y − 2n`|) + H(c(t − s) − |x + y − 2n`|)} . 2c

Chapter 10

Sturm-Liouville Eigenvalue Problems 10.2

Properties of the Sturm-Liouville operator

1. The verification that L is symmetric is exactly the same as that given (for the general case) on pages 389–390, except for the treatment of the boundary term. We have b du du du −P (x) (x)v(x) = −P (b) (b)v(b) + P (a) (a)v(a) dx du dx a du = −P (b) (b)v(b) (since v(a) = 0) du β1 β1 du = (b) = − u(b) . P (b)u(b)v(b) since β2 dx β2 Therefore, (Lu, v)w =

β1 P (b)u(b)v(b) + β2

b

Z

P (x) a

du dv (x) (x) dx + dx dx

b

Z

R(x)u(x)v(x) dx. a

The expression on the right is symmetric in u and v, which shows that (u, L(v))w = (L(v), u)w = (L(u), v)w , ing that L is symmetric with respect to (·, ·)w . 3. Consider a Robin condition of the form

du (a) = 0 dx (relative to the interval [a, b]). We wish to determine the physically meaningful sign for α1 in the case that this boundary condition models heat flow through the left end of a bar. Recall that the heat flux at x = a (a vector quantity) is given by du −κ (a) dx (the direction of the vector is indicated by the sign). Since the normal direction to the left end of the bar is −1, the rate at which heat flows out of the bar at x = a is α1 u(a) + κ

−κ

du du (a)(−1) = κ (a). dx dx

Thus consider the equation du (a) = γ(u(a) − T ), dx where T is the temperature of the surroundings. This equation states that the rate at which heat flows out of the bar is proportional to the difference between the temperature at the left end of the bar and the temperature of the bar. The above boundary condition is equivalent to κ

−γu(a) + κ

83

du (a) = −γT, dx

84

CHAPTER 10. STURM-LIOUVILLE EIGENVALUE PROBLEMS which shows that we should take α1 = −γ in the Robin boundary condition. Thus α1 < 0 is the physically meaningful case for a Robin boundary condition at the left endpoint. At the right endpoint x = b, the normal direction is 1, and therefore the rate at which heat flows out of the bar is κ

du (b). dx

The same reasoning as above then shows that the boundary condition is γu(b) + κ

du (b) = γT. dx

Therefore, we should take β1 = γ, which shows that β1 > 0 is the physically meaningful case. 5. The eigenpairs are λn = 3 + 7.

nπ 2 ln 2

, un (x) = sin

nπ ln x , n = 1, 2, 3, . . . . ln 2

(a) The eigenpairs are λn = 4 + β +

nπ 2 ln 2

, un (x) = x2 sin



nπ ln x , n = 1, 2, 3, . . . . ln 2

(b) 0 is an eigenvalue if and only if β = −4 −

nπ 2 ln 2

for some positive integer n. (c) There is exactly one negative eigenvalue if and only if −4 −

4π 2 π2 ≤ β < −4 − , 2 (ln 2) (ln 2)2

while there are exactly k negative eigenvalues if and only if −4 −

(k + 1)2 π 2 k2 π2 ≤ β < −4 − . (ln 2)2 (ln 2)2

9. The eigenpairs are λn = 2 +

10.3

nπ 2 ln 2

, un (x) = x−1



nπ cos ln 2



nπ ln x ln 2



+ sin

nπ ln x ln 2

, n = 1, 2, 3, . . . .

Numerical Methods for Sturm-Liouville problems

1. We wish to apply the finite element method to du d P (x) + R(x)u = λw(x)u, a < x < b, − dx dx u(a) = 0, du (b) = 0. dx Using the usual uniform mesh on the interval [a, b], the finite element space of consists of all piecewise linear functions v (relative to the given mesh) satisfying v(a) = 0 (note that the Neumann condition is a natural boundary condition). The standard basis for this space is {φ1 , . . . , φn } (where each φj is the usual “hat” function). Applying the Galerkin method to the weak form (Equation (10.12) in the text), we obtain Au = λWu, where A and W are the n × n matrices defined by Z b dφj dφi Aij = P (x) (x) (x) + R(x)φj (x)φi (x) dx, dx dx a Z b Wij = w(x)φj (x)φi (x) dx. a

(These are the same formulas derived in the text, except now i and j vary from 1 to n.)

10.4. EXAMPLES OF STURM-LIOUVILLE PROBLEMS

85

3. The result is exactly the same as in Exercise 1, but now the basis for the finite element space is {φ0 , φ1 , . . . , φn }, the matrices A and W are (n + 1) × (n + 1), and i, j vary from 0 to n in the formulas for Aij , Wij . 5. Let λ ∈ C and x ∈ Cn satisfy Ax = λWx, and let (x, y)W = (Wx) · y. We can assume x satisfies (x, x)W = 1. We then obtain λ = λ(x, x)W = (λWx) · x = (Ax) · x = x · (Ax) = x · (λWx) = λ x · (Wx) = λ ((Wx) · x) = λ. (Notice how the symmetry of both A and W is used to move these matrices from one side of the inner products to the other.) Since λ = λ, it follows that λ ∈ R.

10.4

Examples of Sturm-Liouville problems

1. Let u(x) =

n X

αi φi (x)

i=0

be a continuous piecewise linear function defined on [0, `], relative to the usual uniform mesh. Here {φ0 , φ1 , . . . , φn } is the standard basis. The average value of u on [0, `] is then Z Z Z ` n n 1 ` 1X 1 `X αi φi (x) dx = u(x) dx = αi φi (x) dx. ` 0 ` 0 i=0 ` i=0 0 Since the graph of each φi , i = 1, . . . , n − 1, forms a triangle of height 1 and base 2h (h = `/n), it follows that Z ` φi (x) dx = h, i = 1, 2, . . . , n − 1; 0

similarly, Z

`

`

Z φ0 (x) dx =

φn (x) dx =

0

0

h , i = 1, 2, . . . , n − 1. 2

Thus the average value of u on [0, `] is h `



1 1 α0 + α1 + · · · + αn−1 + αn . 2 2

3. The fundamental frequency is about 3254 Hz.

10.5

Robin boundary conditions

1. The solution to the IBVP is u(x, t) =

∞ X

an (t) cos

n=1

(2n − 1)πx , 2`

where 2

an (t) = cn e−κ(2n−1)

π 2 t/(4`2 ρc)

, cn =

40(−1)n+1 , n = 1, 2, . . . . (2n − 1)π

Some snapshots of the solution are shown in Figure 10.1. 3. Consider the following Sturm-Liouville problem: d2 u = λu, 0 < x < `, dx2 u(0) = 0,

−

κ

du (`) + αu(`) = 0 dx

where κ > 0, α < 0. We can write the second boundary condition as du (`) − αu(`) = 0, dx where α = −α/κ > 0. The analysis of the eigenvalues turns out to depend on the value of α`. We do the analysis by considering three cases:

86

CHAPTER 10. STURM-LIOUVILLE EIGENVALUE PROBLEMS 11 10 9 8 7 6 5 4 3 2 1 0 0

20

40

60

80

100

Figure 10.1: The temperature in the bar of Exercise 10.5.1 after 30, 60, 90, 120, 150, 180 minutes. Case 1 (λ > 0) The positive eigenvalues are the solutions of the equation √ √ λ tan ( λ`) = α or s , s > 0, tan (s) = α` √ . where s = λ`. There are solutions sk = (2k + 1)π/2, k = 1, 2, 3, . . .. The corresponding eigenpairs are s x s2 k , k = 1, 2, 3, . . . . λk = 2k , ψk (x) = sin ` ` These eigenpairs exist regardless of the value of α`. In addition, if α` < 1, then there is an additional solution of tan (s) = s/(α`), namely, s0 ∈ (0, π/2), with corresponding eigenvalue λ0 = s20 /`2 and eigenfunction ψ0 (x) = sin (s0 x/`). Case 2 (λ = 0) If α` = 1, then λ0 = 0 is an eigenvalue with eigenfunction ψ0 (x) = x. If α` 6= 1, then 0 is not an eigenvalue. Case 3 (λ < 0) Negative eigenvalues must satisfy the equation √ √ −λ tanh ( −λ`) = α or s , s > 0, tanh (s) = α` √ where s = −λ`. If α` ≤ 1, then this equation has no solution and hence there are no negative eigenvalues. If α` > 1, then there is a single solution s0 , leading to a single negative eigenvalue λ0 = −s20 /`2 and eigenfunction ψ0 (x) = sinh (s0 x/`). Therefore, when α < 0, there is a sequence of eigenpairs s x s2 k λk = 2k , ψk (x) = sin , k = 1, 2, 3, . . . , ` ` where 2 2 . (2k + 1) π λk = . 2 4` √ There is also an eigenvalue λ0 , which is positive if α` < 1 (with eigenfunction ψ0 (x) = sin √ ( λ0 x)), zero if α` = 1 (with eigenfunction ψ0 (x) = x), and negative if α` > 1 (with eigenfunction ψ0 (x) = sinh ( −λ0 x)).

10.6

Finite element methods for Robin boundary conditions

1. The finite element formulation of −

d dx

du P (x) + R(x)u = F (x), a < x < b, dx u(a) = 0, du β (b) + αu(b) = 0, dx

10.6. FINITE ELEMENT METHODS FOR ROBIN BOUNDARY CONDITIONS

87

where P (x) > 0 for all x ∈ [a, b], α, β > 0, is nearly the same as for the BVP (10.26) given in the text. The Dirichlet condition at x = a is an essential boundary condition (whereas the Neumann condition in (10.26) is a natural boundary condition). This implies that we only consider finite element functions satisfying the Dirichlet condition, so we use the basis {φ1 , φ2 , . . . , φn } instead of {φ0 , φ1 , . . . , φn }. The weak form of the BVP is unchanged from (10.27), although now du P (a) (a)v(a) dx vanishes because v(a) = 0 rather than because du/dx(a) = 0. The result is the linear system Au = f , where A and f are given by the same formulas as in the text (page 420), except that the entries corresponding to φ0 are omitted. Thus A is now n × n and f is now an n-vector. 3. The weak form of −

d dx

P (x)

du dx

+ R(x)u = F (x), a < x < b,

du (a) = 0, dx du β1 u(b) + β2 (b) = 0 dx

−α1 u(a) + α2

is P (a)

du du (a)v(a) − P (b) (b)v(b) + dx dx

b

Z

P (x) a

du dv (x) (x) dx + dx dx

b

Z

b

Z R(x)u(x)v(x) dx =

a

F (x)v(x) dx a

for all test functions v. The boundary conditions imply du α1 du β1 (a) = u(a), (b) = − u(b), dx α2 dx β2 and therefore the weak form can be written as α1 β1 P (a)u(a)v(a) + P (b)u(b)v(b) + α2 β2

Z

b

P (x) a

dv du (x) (x) dx + dx dx

b

Z

b

Z R(x)u(x)v(x) dx =

a

F (x)v(x) dx. a

Notice that the Robin conditions at the endpoints are natural boundary conditions, and hence no boundary conditions are imposed on the test functions. We therefore take {φ0 , φ1 , . . . , φn } as the basis for the finite element space (using the usual notation), and obtain the linear system Au = f , where A = K + M + G. The matrices K and M are the same stiffness and mass matrices defined in the text (page 420), and G ∈ R(n+1)×(n+1) is defined by  α1  α2 P (a), i = j = 0, β1 Gij = P (b)β, i = j = n,  β2 0, otherwise. 5. Suppose u satisfies −κ

du = F (x), 0 < x < `, dx

du (0) = 0, dx κ

du (`) + αu(`) = 0. dx

Then Z

`

`

Z F (x) dx = −κ

0

0

d2 u (x) dx = −κ dx2



du du (`) − (0) dx dx du du = −κ (`) since (0) = 0 dx dx = αu(`) (applying the Robin condition at x = `).

88

CHAPTER 10. STURM-LIOUVILLE EIGENVALUE PROBLEMS

10.7

The theory of Sturm-Liouville problems: an outline

1. Let L : U → V be a compact linear operator, where U , V are normed linear spaces. We wish to show that L is bounded. We will argue by contradiction and assume L is unbounded. Then there exists a sequence {un } ⊂ U with kun kU = 1 for all n and kLun kV → ∞ as n → ∞. (Here we are using the alternate definition of the norm of L: kLk = sup{kLukV : u ∈ U, kukU = 1}.) Since L is compact and {un } is bounded in U , there must exist a subsequence {unk } of {un } such that {Lunk } is convergent in V . But kLun kV → ∞ implies that kLunk kV → ∞, and hence that {Lunk } is not convergent. This is a contradiction, which shows that the assumption that L is unbounded is not possible. Therefore, L must be bounded, which is what we wanted to prove. 3. Suppose that, for a given λ ∈ R, u = φ, u = ψ both satisfy du d − P (x) + R(x)u = λw(x)u, a < x < b, dx dx du α1 u(a) + α2 (a) = 0, dx du β1 u(b) + β2 (b) = 0, dx where P and w are assumed to be positive on the interval [a, b]. Let φ(x) dφ dψ ψ(x) W (x) = dφ dψ = φ(x) dx (x) − dx (x)ψ(x) (x) (x) dx dx be the Wronskian of φ, ψ. Since φ, ψ satisfy the linear homogeneous ODE d du − P (x) + (R(x) − λw(x)u = 0 dx dx on [a, b], we can show that {φ, ψ} is linearly independent by proving that W (x) = 0 for any given x ∈ [a, b]. We have dφ dφ α1 α1 φ(a) + α2 (a) = 0 ⇒ (a) = − φ(a), dx dx α2 and similarly dψ α1 (a) = − ψ(a). dx α2 Therefore, dφ α1 dψ α1 (a) − (a)ψ(a) = − φ(a)ψ(a) + W (a) = φ(a) φ(a)ψ(a) = 0. dx dx α2 α2 This shows that {φ, ψ} is linearly dependent, as desired.

Chapter 11

Problems in Multiple Spatial Dimensions 11.1

Physical models in two or three spatial dimensions

1. Let Ω ⊂ R2 be the rectangular domain Ω = x ∈ R2 : a < x1 < b, c < x2 < d , and let F : R2 → R2 be a smooth vector field. Then Z Z bZ d ∂F1 ∂F2 ∇·F= (x1 , x2 ) + (x1 , x2 ) dx2 dx1 ∂x1 ∂x2 Ω a c Z bZ d Z bZ d ∂F1 ∂F2 = (x1 , x2 ) dx2 dx1 + (x1 , x2 ) dx2 dx1 ∂x 1 a a c c ∂x2 Z dZ b Z bZ d ∂F1 ∂F2 = (x1 , x2 ) dx1 dx2 + (x1 , x2 ) dx2 dx1 c a ∂x1 a c ∂x2 Z d Z b = {F1 (b, x2 ) − F1 (a, x2 )} dx2 + {F2 (x1 , d) − F2 (x2 , c)} dx1 c

a d

Z = Zc =

Z b Z b F1 (b, x2 ) dx2 − F1 (a, x2 ) dx2 + F2 (x1 , d) dx1 − F2 (x2 , c) dx1 c a a Z Z Z Z F·n+ F·n+ F·n+ F·n= F · n, Z

Γ2

d

Γ4

Γ3

Γ1

∂Ω

as desired. Here Γ1 , Γ2 , Γ3 , Γ4 denote the four sides of ∂Ω, labeled as in Figure 11.3 on page 452 of the text. 3. We have ∇ · F(x) = 1, and therefore Z ∇ · F = area(Ω) = π. Ω

On the other hand, on ∂Ω, n = x and therefore F · n = 2x1 x2 + x22 . In polar coordinates, F · n = sin (2θ) + sin2 (θ), and hence Z Z 2π sin (2θ) + sin2 (θ) dθ = π. F·n= ∂Ω

0

This verifies the divergence theorem for this domain and vector field. 5. Let L : C 2 (Ω) → C(Ω) be defined by Lu = −∇ · (k(x)∇u). By the product rule, ∇ · (k(∇u)v) = k∇u · ∇v + ∇ · (k∇u)v; therefore, Z

Z

−

Z

∇ · (k∇u)v = Ω

k∇u · ∇v − Ω

89

∇ · (k(∇u)v). Ω

90

CHAPTER 11. PROBLEMS IN MULTIPLE SPATIAL DIMENSIONS By the divergence theorem, Z

Z

Z

∇ · (k(∇u)v) =

k(∇u)v) · n =

Ω

∂Ω

k ∂Ω

∂u v. ∂n

The boundary term disappears if u and v both satisfy Dirichlet condition (since then v = 0 on ∂Ω) and also if u and v both satisfy Neumann conditions (since then ∂u/∂n = 0 on ∂Ω). Thus, in either case, Z Z (Lu, v) = − ∇ · (k∇u)v = k∇u · ∇v. Ω

Ω

Since the right-hand side is symmetric in u and v, so is the left-hand side; that is, (Lu, v) = (Lv, u) = (u, Lv), and we see that L is a symmetric operator. 7. We have ∂ ∂ ∂ (F1 (u(x))) + (F2 (u(x))) + (F3 (u(x))) ∂x1 ∂x2 ∂x3 ∂u ∂u ∂u = F10 (u(x)) (x) + F20 (u(x)) (x) + F30 (u(x)) (x) ∂x1 ∂x2 ∂x3 = F0 (u(x)) · ∇u(x).

∇ · F(u(x)) =

9. Using the product rule for scalar functions, we obtain ∇ · (v∇u) =

X 3 3 X ∂ ∂v ∂u ∂u ∂2u v = +v 2 ∂xi ∂xi ∂xi ∂xi ∂xi i=1 i=1

=

3 3 X X ∂2u ∂v ∂u +v ∂xi ∂xi ∂x2i i=1 i=1

=

∇v · ∇u + v∆u.

2 11. Suppose u, v ∈ Cm (Ω). Then

Z

Z

(Lm u, v) = −

v∆u

Z ∇u · ∇v −

=

Ω

Ω

v ∂Ω

∂u (Green’s first identity) ∂n

Z ∇u · ∇v

= ZΩ =

u ∂Ω

∂v − ∂n

Z u∆v (Green’s first identity) Ω

Z =

−

u∆v = (u, Lm v). Ω

The boundary vanish because the product that forms the integrand is zero over the entire boundary. For example, Z ∂u v =0 ∂Ω ∂n since v = 0 on Γ1 and ∂u/∂n = 0 on Γ2 .

11.2 1.

Fourier series on a rectangular domain (a) cmn = −

4 (5 + 7(−1)m ) (−1 + (−1)n ) m3 n3 π 6

(b) The graphs of the error are given in Figure 11.1. 3. The solution is given by u(x) =

∞ X ∞ X cmn sin (mπx1 ) sin (nπx2 ), λ mn m=1 n=1

where the cmn are the coefficients from Exercise 1 and λmn = (m2 + n2 )π 2 .

11.2. FOURIER SERIES ON A RECTANGULAR DOMAIN

91

−3

x 10 2 0 −2 1

0.5 0

0

0.2

x2

0.4

0.8

1

0.6

0.8

1

0.6

x1

−3

x 10 2 0 −2 1

0.5 0

0

0.2

x2

0.4 x1

Figure 11.1: The error in approximating f (x, y) by the first 4 (top) and the first 25 (bottom) of the double Fourier sine series. (See Exercise 11.2.1.)

5.

(a) The IBVP is ρc

∂u − κ∆u ∂t u(x, 0)

The domain Ω is the rectangle x ∈ R2

0.02, x ∈ Ω, t > 0,

=

5, x ∈ Ω,

0, x ∈ ∂Ω, t > 0. : 0 < x1 < 50, 0 < x2 < 50 .

u(x, t)

=

=

(b) The solution is ∞ X ∞ X

u(x, t) =

amn (t) sin

m=1 n=1

mπx 1

50

sin

nπx 2

50

,

where amn (t)

=

bmn

=

cmn

=

λmn

=

cmn cmn e−κλmn t/(ρc) + , κλmn κλmn 20 (−1 + (−1)m ) (−1 + (−1)n ) , mnπ 2 2 (−1 + (−1)m ) (−1 + (−1)n ) , 25mnπ 2 (m2 + n2 )π 2 . 2500 bmn −

(c) The steady-state temperature us (x) satisfies the BVP

The solution is us (x) =

−κ∆u

=

0.02, x ∈ Ω,

u(x)

=

0, x ∈ ∂Ω.

∞ X ∞ mπx nπx X cmn 1 2 sin sin . κλ 50 50 mn m=1 n=1

(d) The maximum difference between the temperature after 10 minutes and the steady-state temperature is about 1 degree. The difference is graphed in Figure 11.2. 7. The minimum temperature in the plate reaches 4 degrees Celsius after 825 seconds.

√ 9. The leading edge of the wave is initially 2/5 units from the boundary, wave travels at a speed of 261 2 √ and the . units per second. Therefore, the wave reaches the boundary after 2/1305 = 0.00108 seconds. Figure 11.6 from the text shows the wave about to reach the boundary after 10−3 seconds.

92

CHAPTER 11. PROBLEMS IN MULTIPLE SPATIAL DIMENSIONS

1 0.8 0.6 0.4 0.2 0 60 60

40 40 20

20 0

x2

0

x1

Figure 11.2: The difference between the temperature in the plate after 10 minutes and the steady-state temperature. (See Exercise 11.2.5.) 11. The difficult task is to compute the Fourier coefficients of the initial displacement ψ. If we write ψ(x) =

∞ X ∞ X

bmn sin (mπx1 ) sin (nπx2 ),

m=1 n=1

then we find that

( bmn =

We then have that u(x, t) =

0, sin2 (mπ/2) , 1250m2 π 2

∞ X ∞ X

m 6= n, m = n.

amn (t) sin (mπx1 ) sin (nπx2 ),

m=1 n=1

where amn (t) satisfies d2 amn + λmn c2 amn dt2 amn (0) damn (0) dt The result is

amn (t) =

Thus the solution is u(x, t) =

∞ X

=

0,

=

bmn ,

=

0.

0, √ bmm cos (c λmm t),

m 6= n, m = n.

√ bmm cos (c λmm t) sin (mπx1 ) sin (mπx2 ).

m=1

Four snapshots of u are shown in Figure 11.3. 13. The solution is u(x) =

∞ X ∞ X 4 sin (mπ/2) sin (nπ/2) sin (mπx1 ) sin (nπx2 ). (m2 + n2 ) π 2 m=1 n=1

The graph of u is shown in Figure 11.4. 15.

(a) The IBVP is ρc

∂u − κ∆u ∂t u(x, 0) u(x, t) ∂u (x, t) ∂n

=

0.02, x ∈ Ω, t > 0,

=

5, x ∈ Ω,

=

0, x ∈ Γ1 ∪ Γ4 , t > 0,

=

0, x ∈ Γ2 ∪ Γ3 , t > 0.

The domain Ω is the rectangle

x ∈ R2 : 0 < x1 < 50, 0 < x2 < 50 .

11.2. FOURIER SERIES ON A RECTANGULAR DOMAIN −4

−4

x 10

x 10

1

1

0

0

−1 1 0.5 x

2

−1 1

1 0 0

0.5 x

0.5 x1

2

−4

0.5 x1

0 0

0.5 x1

x 10 1

0

0

−1 1

2

1 0 0

−4

x 10 1

0.5 x

93

−1 1

1

0.5 x

0.5 0 0

x1

2

1

Figure 11.3: Four snapshots of the solution to Exercise 11.2.11: t = 0 (top left), t = T /8 (top right), t = T /2 (bottom left), t = 5T /8 (bottom right). The solution is periodic with period T .

1

0.5

0

−0.5 1 1 0.5 0.5 0

x2

0

x1

Figure 11.4: The solution to Exercise 11.2.13. (b) The solution is u(x, t) =

∞ X ∞ X

amn (t) sin

m=1 n=1

(2m − 1)πx1 100



sin

(2n − 1)πx2 , 100

where amn (t)

=

bmn

=

cmn

=

λmn

=

cmn cmn , bmn − e−κλmn t/(ρc) + κλmn κλmn 80 (2m − 1)(2n − 1)π 2 8 25(2m − 1)(2n − 1)π 2

((2m − 1)2 + (2n − 1)2 )π 2 . 10000

(c) The steady-state temperature us (x) satisfies the BVP

The solution is us (x) =

−κ∆u

=

0.02, x ∈ Ω,

u(x) ∂u (x) ∂n

=

0, x ∈ Γ1 ∪ Γ4 ,

=

0, x ∈ Γ2 ∪ Γ3 .

∞ X ∞ X (2m − 1)πx1 (2n − 1)πx2 cmn sin sin . κλmn 100 100 m=1 n=1

94

CHAPTER 11. PROBLEMS IN MULTIPLE SPATIAL DIMENSIONS (d) After 10 minutes, the temperature u is not very close to the steady-state temperature; the difference u(x, y, 600) − us (x, y) is graphed in Figure 11.5. (Note: The temperature variation in this problem may be outside the range in which the linear model is valid, so these results may be regarded with some skepticism.)

0 −2 −4 −6 −8 −10 60 60

40 40 20 x2

20 0 0

x1

Figure 11.5: The difference between the temperature in the plate after 10 minutes and the steady-state temperature. (See Exercise 11.2.15.)

11.3

Fourier series on a disk

1. The next three coefficients in the series for g are . . . c40 = −0.0209908, c50 = 0.0116362, c60 = −0.00722147. 3. The solution is u(r, θ) =

∞ X

am0 J0 (s0m r),

m=1

where am0 =

cm0 , s20m

the cm0 are the coefficients of f (r, θ) = 1 − r, and the s0m are the positive roots of J0 . A direct calculation shows that a10 a20 a30 a40 a50 a60

. = . = . = . = . = . =

0.226324, −0.0157747, 0.00386078, −0.00148156, 0.000716858, −0.000399554.

The solution (approximated by six of the series) is graphed in Figure 11.6. 5. The solution is u(r, θ) =

∞ X

am0 J0 (s0m r),

m=1

where am0 =

cm0 , s20m

11.3. FOURIER SERIES ON A DISK

95

0.25 0.2 0.15 0.1 0.05 0 1 1

0.5 0.5

0 0

−0.5

−0.5 −1 −1

x2

x1

Figure 11.6: The approximation to solution u, computed with six of the (generalized) Fourier series (see Exercise 11.3.3).

the cm0 are the coefficients of f (r, θ) = r, and the s0m are the positive roots of J0 . A direct calculation shows that . = . = . = . = . = . =

a10 a20 a30 a40 a50 a60

0.141350, −0.0371986, 0.0106599, −0.00537260, 0.00283289, −0.00182923.

The solution (approximated by six of the series) is graphed in Figure 11.7.

0.15

0.1

0.05

0 1 1

0.5 0.5

0 0

−0.5 x2

−0.5 −1 −1

x1

Figure 11.7: The approximation to solution u, computed with six of the (generalized) Fourier series (see Exercise 11.3.5). 7. If we write φ(r, θ) =

∞ X

cm0 J0 (αm0 r) +

m=1

∞ X ∞ X

(cmn cos (nθ) + dmn sin (nθ)) Jn (αmn r),

m=1 n=1

where the cmn and dmn are known, and u(r, θ, t) =

∞ X

am0 (t)J0 (αm0 r) +

m=1

∞ X ∞ X

(amn (t) cos (nθ) + bmn (t) sin (nθ)) Jn (αmn r),

m=1 n=1

then amn (t)

=

cmn e−κλmn t/(ρc) , m = 1, 2, 3, . . . , n = 0, 1, 2, . . . ,

bmn (t)

=

dmn e−κλmn t/(ρc) , m, n = 1, 2, 3, . . . .

96

CHAPTER 11. PROBLEMS IN MULTIPLE SPATIAL DIMENSIONS However, since φ(r, θ) = r(1 − r) cos (θ)/5 (a function of r times cos (θ)), all of the coefficients of φ are zero except for cm1 , m = 1, 2, 3, . . .. Therefore, u(r, θ, t) =

∞ X

am1 (t) cos (θ)J1 (αm1 r),

m=1

with am1 (t) given above. After 30 seconds, the temperature distribution can be approximated accurately using a single eigenfunction (corresponding to the largest eigenvalue, λ11 ), so we need only . c11 = 9.04433. The temperature distribution after 30 seconds is graphed in Figure 11.8.

0.04 0.02 0 −0.02 −0.04 10 10

5 5

0 0

−5 x2

−5 −10 −10

x1

Figure 11.8: The approximation to solution u at t = 30, computed with one term of the (generalized) Fourier series (see Exercise 11.3.7). 9. A circular drum of radius A has area πA2 , the same area as a square drum of side length frequency of such a square drum is q 2 2 c (√ππA)2 + (√ππA)2 c c 1 . = 0.398942 . = √ 2π A 2π A

√ πA. The fundamental

Comparing to Example 11.9, we see that the circular drum sounds a lower frequency than a square drum of equal area. . . . . 11. s41 = 7.588342434503804, s42 = 11.06470948850118, s43 = 14.37253667161759, s44 = 17.61596604980483

11.4

Finite elements in two dimensions

1. The mesh for this problem is shown in Figure 11.9, in which the free nodes are labeled. The stiffness matrix is   4.4815 −1.1759 −1.1759 0.0000 4.9259 0.0000 −1.3426  .  −1.1759 , K=  −1.1759 0.00000 4.9259 −1.3426  0.0000 −1.3426 −1.3426 5.8148 while the load vector is



 0.11111 .  0.11111   F=  0.11111  . 0.11111

The resulting weights for the finite element approximation are given by   0.048398 .  0.044979   u = K−1 F =   0.044979  . 0.039879

11.4. FINITE ELEMENTS IN TWO DIMENSIONS

97

1

0.8 3

4

1

2

x2

0.6

0.4

0.2

0 0

0.2

0.4

0.6

0.8

1

x1

Figure 11.9: The mesh for Exercise 11.4.1.

3. The mesh for this problem is shown in Figure 11.10, in which the free nodes are labeled. The stiffness matrix is 16 × 16 and neither it nor the load vector will be reproduced here. The matrix is singular, as is expected for a Neumann problem. The unique solution to KU = F with its last component equal to zero is              .  u=             

1

−0.071138 −0.047315 −0.012656 0.001458 −0.067934 −0.047349 −0.014145 0.001566 −0.065432 −0.046198 −0.015160 −0.000329 −0.062914 −0.046391 −0.016065 0.000000

              .             

13

14

15

16

9

10

11

12

5

6

7

8

0.8

x2

0.6

0.4

0.2

0 0

1

2

0.2

3

0.4

0.6

4

0.8

1

x1

Figure 11.10: The mesh for Exercise 11.4.3.

5. Let Vn be the space of continuous piecewise linear functions relative to a given triangulation of the domain Ω, and 2 let {φ1 , . . . , φn } be the standard basis for Vn . Suppose that u : Ω → PR belongs to L (Ω). We wish to find v ∈ Vn that minimizes kw − ukL2 (Ω) over all w ∈ Vn . We can write v = n α φ and use the projection theorem: v i=1 i i must satisfy (u − v, w)L2 (Ω) = 0 for all w ∈ Vn , which is equivalent to (u − v, φi )L2 (Ω) = 0 for all i = 1, 2, . . . , n.

98

CHAPTER 11. PROBLEMS IN MULTIPLE SPATIAL DIMENSIONS We obtain the following equations: (v, φi ) = (u, φi ), i = 1, 2, . . . , n ! n X ⇒ αj φj , φi = (u, φi ), i = 1, 2, . . . , n j=1

⇒

n X

αj (φj , φi ) = (u, φi ), i = 1, 2, . . . , n.

j=1

This last system of equations is equivalent to Mα = f , where M is the usual mass matrix, Z Mij = φj φi , i, j = 1, 2, . . . , n, Ω

and f is the vector defined by Z fi =

vφi , i = 1, 2, . . . , n. Ω

7.

(a) A direct calculation shows that the inhomogeneous Dirichlet problem −∇ · (k(x)∇u)

=

f (x), x ∈ Ω,

u(x)

=

g(x), x ∈ ∂Ω,

has weak form 2 u ∈ G + V, a(u, v) = (f, v) for all v ∈ V = CD (Ω),

where G is any function satisfying the condition G(x) = g(x) for x ∈ ∂Ω, and G + V = {G + v : v ∈ V } . Substituting u = G + w, where w ∈ V , into the weak form yields w ∈ V, a(w, v) = (f, v) − a(G, v) for all v ∈ V. In the Galerkin method, we replace V by a finite-dimensional subspace Vn and solve w ∈ Vn , a(w, v) = (f, v) − a(G, v) for all v ∈ Vn . When using piecewise linear finite elements, we can satisfy the boundary conditions approximately by taking G to be a continuous piecewise linear function whose values at the boundary nodes agree with the given boundary function g. For simplicity, we take G to be zero at the interior nodes. The resulting load vector is then given by fi = (f, φi ) − a(G, φi ), i = 1, 2, 3, . . . , n. Since G is zero on interior nodes, the quantity a(G, φi ) is nonzero only if (free) node i belongs to a triangle adjacent to the boundary. (b) The regular triangulation of the unit square having 18 triangles has only four interior (free) nodes, and each one belongs to triangles adjacent to the boundary. This means that every entry in the load vector is modified (which is not the typical case). Since f = 0, the load vector f is defined by fi = −a(G, φi ), i = 1, 2, 3, 4. The load vector is

while the solution to Ku = f is



 0.22222222222222 .  1.55555555555556   f =  0.22222222222222  , 1.55555555555556 

 0.27777777777778 .  0.61111111111111   u=  0.27777777777778  . 0.61111111111111

11.4. FINITE ELEMENTS IN TWO DIMENSIONS 9.

99

(a) Consider the BVP −∇ · (k(x)∇u) = f (x), in Ω, ∂u (x) = h(x), x ∈ ∂Ω. ∂n

(11.1)

Multiplying the PDE by a test function v ∈ V˜ = C 2 (Ω) and applying Green’s first identity yields Z Z Z k∇u · ∇v = fv + kvh for all v ∈ V˜ . Ω

Ω

∂Ω

We will apply Galerkin’s method with V˜m equal to the space of all continuous piecewise linear function relative to a given triangulation T . We label the standard basis of V˜m as φ1 , φ2 , . . . , φn , and the nodes as z1 , z2 , . . . , zn . (Every node in the mesh is now free.) Thus V˜n = span{φ1 , φ2 , . . . , φn }. We then define the Galerkin approximation wn by Z wn ∈ V˜n , a(wn , φi ) = (f, φi ) + khφi , i = 1, 2, . . . , m. ∂Ω

Writing wn =

Pm

j=1 βj φj and



 β1  β2    w =  . ,  ..  βn ˜ ˜ ˜ we obtain the system Kw = F. The stiffness matrix K ∈ Rm×m is given by ˜ ij = a(φj , φi ), i, j = 1, 2, . . . , n, K ˜ ∈ Rm by and the load vector F F˜i = (f, φi ) +

Z khφi , i = 1, 2, . . . , n. ∂Ω

The boundary integral in the expression for F˜i is zero unless zi is a boundary node. (b) The compatibility condition for (11.1) is determined by the following calculation: Z Z Z Z ∂u f =− ∇ · (k∇u) = − k =− kh. ∂Ω Ω ∂Ω ∂n Ω (c) Now consider (11.1) with 3 3x2 , h(x) = − 1 , 4 5 where Ω is the unit square. We will produce the finite element solution using the regular grid with 18 ˜ ∈ R16×16 is singular, as should be expected, and there are triangles (16 nodes). The stiffness matrix K infinitely many solutions. The unique solution u with last component equal to zero is   0.3622  0.2984     0.1344     −0.0843     0.3728     0.3302     0.2081     0.0255  .  u=   0.3803   0.3447     0.2334     0.0591     0.3665     0.3257     0.1878  k(x) = 1, f (x) = x1 x2 +

0.0000

100

CHAPTER 11. PROBLEMS IN MULTIPLE SPATIAL DIMENSIONS

11.5

The free-space Green’s function for the Laplacian

1. Following the hint and using the trigonometric identities cos (a − b) = cos (a) cos (b) + sin (a) sin (b), sin (a − b) = sin (a) cos (b) − cos (a) sin (b), we obtain y1 = r cos (θ − α) = r cos (θ) cos (α) + r sin (θ) sin (α) = x1 cos (α) + x2 sin (α), y2 = r sin (θ − α) = r sin (θ) cos (α) − r cos (θ) sin (α) = −x1 sin (α) + x2 cos (α),

as desired. 3. We have −

1 dψ d2 ψ dψ d2 ψ − = φ(r) ⇒ −r − = rφ(r) dr2 r dr dr2 dr d dψ ⇒− r (r) = rφ(r). dr dr

This last ODE can be solved by integrating twice to yield the general solution Z r r ψ(r) = c1 + c2 ln (r) − ln φ(s)s ds. s 0 5. Let V be the space of all smooth functions defined on R2 and having compact . We wish to show that the variational problem Z ∇u(x) · ∇v(x) dx = v(y) for all v ∈ V R2

is equivalent to Z 0

2π

Z 0

∞



1 ∂u ∂v ∂u ∂v + 2 ∂r ∂r r ∂θ ∂θ

v rdrdθ = v(y) for all v ∈ V

in polar coordinates. (a) We first convert Z ∇u(x) · ∇v(x) dx R2

to polar coordinates, where y ∈ R2 is fixed and represents the origin for the polar coordinates. By the chain rule, we have sin (θ) ∂u ∂u cos (θ) ∂u ∂u ∂u ∂u = cos (θ) − , = sin (θ) + ∂x1 ∂r r ∂θ ∂x2 ∂r r ∂θ (cf. the derivation in Section 11.3.1) and similarly for v. Therefore, sin (θ) ∂u sin (θ) ∂v ∂u ∂v ∇u · ∇v = cos (θ) − cos (θ) − + ∂r r ∂θ ∂r r ∂θ cos (θ) ∂u cos (θ) ∂v ∂u ∂v sin (θ) + sin (θ) + ∂r r ∂θ ∂r r ∂θ sin (θ) cos (θ) ∂u ∂v sin (θ) cos (θ) ∂u ∂v sin2 (θ) ∂u ∂v ∂u ∂v − − + + ∂r ∂r r ∂θ ∂r r ∂r ∂θ r2 ∂θ ∂θ 2 sin (θ) cos (θ) ∂u ∂v sin (θ) cos (θ) ∂u ∂v cos (θ) ∂u ∂v ∂u ∂v sin2 (θ) + + + ∂r ∂r r ∂θ ∂r r ∂r ∂θ r2 ∂θ ∂θ 2 2 sin (θ) cos (θ) ∂u ∂v ∂u ∂v 2 2 = (sin (θ) + cos (θ)) + + ∂r ∂r r2 r2 ∂θ ∂θ ∂u ∂v 1 ∂u ∂v = + 2 . ∂r ∂r r ∂θ ∂θ = cos2 (θ)

11.5. THE FREE-SPACE GREEN’S FUNCTION FOR THE LAPLACIAN

101

It follows that Z

2π

Z

∞

Z



∇u(x) · ∇v(x) dx = R2

0

0

∂u ∂v 1 ∂u ∂v + 2 ∂r ∂r r ∂θ ∂θ

v rdrdθ,

as desired. (b) We can derive the same result by starting with the PDE −∆u = δ(x − y) in polar coordinates, multiplying by a test function, and integrating. The right-hand side becomes Z δ(x − y)v(x) dx = v(y). R2

On the left side, we have 2π

∞

∂2u 1 ∂2u 1 ∂u − vr dr dθ − ∂r2 r ∂r r2 ∂θ2 0 0 Z 2π Z ∞ 2 Z 2π Z ∞ Z ∞ Z 2π 2 ∂ u ∂u 1 ∂ u − vr dr dθ − v dr dθ − v dθ r dr ∂r2 ∂r r2 0 ∂θ2 0 0 0 0 0 ∞ Z ∞ Z 2π Z ∞ Z 2π Z ∞ ∂u ∂u ∂v ∂u ∂u −r v + r dr + v dr dθ − v dr dθ− ∂r ∂r ∂r ∂r ∂r 0 0 0 0 0 0 ( ) 2π Z 2π Z ∞ ∂u 1 ∂u ∂v v − dθ r dr r2 ∂θ 0 ∂θ ∂θ 0 0 Z 2π Z 2π Z ∞ ∂u ∂v 1 ∂u ∂v r dr dθ + r dr dθ ∂r ∂r r2 ∂θ ∂θ 0 0 0 Z 2π ∂u ∂v 1 ∂u ∂v + 2 r dr dθ. ∂r ∂r r ∂θ ∂θ 0 Z

Z



−

= =

= =

Once again, we have derived the desired result. 7. The derivation is long and tedious. As of the time of this writing, it can be found at planetmath.org/encyclopedia/DerivationOfTheLaplacianFromRectangularToSphericalCoordinates.html. 9. Let g(ρ) = c1 /ρ. Then c1 ∂g (ρ) = − 2 ∂ρ ρ and therefore Z 0

2π

π

Z 0

∞

Z 0

∂g ∂v 2 ρ sin (φ) dρ dφ dθ = −c1 ∂ρ ∂ρ

Z

2π

Z

2π

Z

0

Z

π

Z

π

Z

0

∞

0

1 ∂v 2 ρ sin (φ) dρ dφ dθ ρ2 ∂ρ

∞

∂v sin (φ) dρ dφ dθ ∂ρ Z ∞ Z 2π Z π ∂v dρ dφ dθ = −c1 sin (φ) ∂ρ 0 0 0 Z 2π Z π = −c1 sin (φ) (−v(y)) dφ dθ = −c1

0

0

0

Z

0

0 2π Z π

= c1 v(y)

sin (φ) dφ dθ 0

0

Z = 2c1 v(y)

2π

dθ 0

= 4πc1 v(y). From this we see that g is a solution of (11.83) if and only if c1 = 1/(4π).

102

CHAPTER 11. PROBLEMS IN MULTIPLE SPATIAL DIMENSIONS

11.6

The Green’s function for the Laplacian on a bounded domain

1. Suppose u satisfies −∆u = f (x) in Ω, α1 u + α2

∂u = 0 on ∂Ω, ∂n

and let v be any test function in C 2 (Ω). Then Z Z Z fv ⇒ − − (∆u)v = Ω

∂Ω

Ω

∂u v+ ∂n

Z

Z ∇u · ∇v =

f v. Ω

Ω

The boundary condition yields ∂u α1 = u on ∂Ω, ∂n α2 Z Z uv + ∇u · ∇v = f v,

− and hence we obtain

α1 α2

Z ∂Ω

Ω

Ω

as desired. 3. Let G be the free-space Green’s function for the Laplacian and define w(x; y) = α1 G(x; y) + α

∂G (x; y) ∂nx

for all x ∈ ∂Ω and y ∈ Ω. Assume that u(x) = vΩ (x; y) satisfies −∆u = 0 in Ω α1 u + α2

∂u = w(x; y) on ∂Ω, ∂n

where w(x; y) = α1 G(x; y) + α2 Define

∂G (x; y). ∂nx

Z u2 (x) =

vΩ (x; y)f (y) dy. Ω

Then, for all v ∈ C 2 (Ω), Z − (∆x vΩ v = 0 ZΩ Z ∂vΩ ⇒ − (x; y)v(x) dσx + ∇x vΩ (x; y) · ∇v(x) dx = 0 ∂nx Ω Z∂Ω Z Z α1 1 ⇒ vΩ (x; y)v(x) dσx − w(x; y)v(x) dσx + ∇x vΩ (x; y) · ∇v(x) dx = 0 α2 ∂Ω α2 ∂Ω Ω Z Z Z α1 1 ⇒ vΩ (x; y)v(x) dσx + ∇x vΩ (x; y) · ∇v(x) dx = w(x; y)v(x) dσx . α2 ∂Ω α2 ∂Ω Ω Notice how we used the boundary condition satisfied by vΩ to substitute for ∂vΩ /∂nx . We now multiply both sides of this last equation by f (y) and integrate over Ω to obtain Z Z Z Z α1 vΩ (x; y)v(x) dσx f (y) dy+ ∇x vΩ (x; y) · ∇v(x) dx f (y) dy α2 Ω ∂Ω Ω Ω Z Z 1 = w(x; y)v(x) dσx f (y) dy. α2 Ω ∂Ω Changing the order of integration in the two integrals on the left yields Z Z Z Z α1 vΩ (x; y)f (y) dy v(x) dσx + ∇x vΩ (x; y)f (y) dy · ∇v(x) dx α2 ∂Ω Ω Ω Ω Z Z 1 = w(x; y)v(x) dσx f (y) dy. α2 Ω ∂Ω

11.6. THE GREEN’S FUNCTION FOR THE LAPLACIAN ON A BOUNDED DOMAIN Since



Z

Z

103

∇x vΩ (x; y)f (y) dy = ∇

vΩ (x; y)f (y) dy

= ∇u2 (x),

Ω

Ω

we obtain α1 α2

Z

Z u2 (x)v(x) dσx + ∂Ω

∇u2 (x) · ∇v(x) dx Z Z 1 w(x; y)v(x) dσx f (y) dy, = α2 Ω ∂Ω Ω

as desired. 5. Let G be the free-space Green’s function for the Laplacian, and suppose vΩ (x; y) satisfies the BVP −∆x vΩ = 0 in Ω, vΩ (x; y) = G(x; y) on ∂Ω for all y ∈ Ω. (Notice that G(x; y) is a continuous function of x ∈ ∂Ω for all y ∈ Ω, and therefore vΩ is welldefined.) We define GΩ (x; y) = G(x; y) − vΩ (x; y), and we will show that GΩ is the Green’s function for the BVP −∆u = f (x) in Ω, u = 0 on ∂Ω. The weak form of the last BVP is 2 u ∈ CD (Ω),

Z

Z

2 f (x)v(x) dx for all v ∈ CD (Ω).

∇u(x) · ∇v(x) dx = Ω

Ω

We prove that GΩ is the desired Green’s function by proving that Z GΩ (x; y)f (y) dy u(x) = Ω

satisfies the above variational problem. We have Z Z G(x; y)f (y) dy − vΩ (x; y)f (y) dy = u1 (x) − u2 (x), u(x) = Ω

Ω

where

Z

Z

u1 (x) =

G(x; y)f (y) dy, u2 (x) = Ω

vΩ (x; y)f (y) dy. Ω

Letting V be the space of all smooth functions with compact , we have Z Z ∇u1 (x) · ∇v(x) dx = f (x)v(x) dx for all v ∈ V, R2

R2

where f has been extended to be zero outside of Ω. This last variational equation holds because G is the free-space 2 Green’s function. For any v ∈ CD (Ω), we can extend v to be identically zero outside of Ω and thereby obtain an element v of V . We then have Z Z 2 ∇u1 (x) · ∇v(x) dx = f (x)v(x) dx for all v ∈ CD (Ω). Ω

Ω

On the other hand, we have −∆u2 (x) = 0 for all x ∈ Ω; multiplying by v ∈

2 CD (Ω)

and integrating yields Z 2 − ∆u2 (x)v(x) dx = 0 for all v ∈ CD (Ω). Ω

Integrating by parts on the left, we obtain Z 2 ∇u2 (x) · ∇v(x) dx = 0 for all v ∈ CD (Ω). Ω

104

CHAPTER 11. PROBLEMS IN MULTIPLE SPATIAL DIMENSIONS (Notice that the boundary term vanishes because v = 0 on ∂Ω.) Putting together the variational equations satisfies by u1 and u2 , we see that u = u1 − u2 satisfies Z Z 2 ∇u(x) · ∇v(x) dx = f (x)v(x) dx for all v ∈ CD (Ω). Ω

Ω

Also, for x ∈ ∂Ω, we have Z Z Z Z u(x) = G(x; y)f (y) dy − vΩ (x; y)f (y) dy = G(x; y)f (y) dy − G(x; y)f (y) dy = 0. Ω

Ω

Ω

Ω

2 This shows that u ∈ CD (Ω); hence u satisfies the given BVP, which shows that GΩ is the desired Green’s function.

7. Green’s second identity is Z



Z (v∆u − u∆v) = Ω

v ∂Ω

∂v ∂u −u ∂n ∂n

.

Let GΩ be the Green’s function for the BVP −∆u = f (x) in Ω α1 u + α2

∂u = 0 on ∂Ω, ∂n

and, for y, z ∈ Ω, define u(x) = GΩ (x; y), v(x) = GΩ (x; z). Since GΩ (x; y) satisfies −∆x GΩ = δ(x − y) in Ω α1 GΩ + α2

∂GΩ = 0, on ∂Ω, ∂nx

we see that Z v∆u = −v(y) = −G(y; z), ZΩ u∆v = −u(z) = −G(z; y). Ω

Therefore, Z (v∆u − u∆v) = G(z; y) − G(y; z). Ω

We also have

∂u α1 = − u on ∂Ω, ∂n α2

or

α1 ∂GΩ (x; y) = − GΩ (x; y). ∂nx α2

Similarly, ∂GΩ α1 (x; z) = − GΩ (x; z). ∂nx α2 Therefore,

Z

v ∂Ω

∂u ∂v −u ∂n ∂n





Z

−

= ∂Ω

α1 α1 GΩ (x; z)GΩ (x; y) + GΩ (x; y)GΩ (x; z) dσx = 0. α2 α2

Thus, by Green’s second identity, G(z; y) − G(y; z) = 0, that is, G(z; y) = G(y; z), as desired. 9. Let Ω ⊂ R2 be given, and let u be the solution of −∆u = 0 in Ω, u = φ(x) on ∂Ω.

11.6. THE GREEN’S FUNCTION FOR THE LAPLACIAN ON A BOUNDED DOMAIN

105

Let v(x) = GΩ (x; y), where GΩ is the Green’s function for the BVP −∆u = f (x) in Ω, u = 0 on ∂Ω. Finally, define Ω = {x ∈ Ω : kx − yk > }, where y ∈ Ω is given and > 0 is small enough that B (y) ⊂ Ω. Notice that ∆u = 0 in Ω, ∆v = 0 in Ω , and v(x) = 0 for all x ∈ ∂Ω. Applying Green’s second identity to u and v yields Z Z ∂u ∂v (v∆u − u∆v) = v −u dσx ∂n ∂n ∂Ω Ω Z Z ∂v ∂u ∂v ∂u (x) − u(x) (x) dσx + v(x) (x) − u(x) (x) dσx ⇒0= v(x) ∂n ∂n ∂n ∂n S (y) ∂Ω Z Z ∂u ∂u ∂GΩ ∂GΩ ⇒0= 0 GΩ (x; y) (x) − φ(x) (x; y) dσx + (x) − u(x) (x; y) dσx ∂n ∂nx ∂n ∂nx ∂Ω S (y) Z Z Z ∂GΩ ∂u ∂GΩ ⇒ φ(x) (x; y) dσx = GΩ (x; y) (x) dσx − u(x) (x; y) dσx . ∂nx ∂n ∂nx ∂Ω S (y) S (y) Now, recall that GΩ = G − vΩ , where G is the free-space Green’s function for the Laplacian and vΩ (x; y) satisfies −∆x vΩ = 0 in Ω, vΩ (x; y) = G(x; y) on ∂Ω. We therefore have Z φ(x) ∂Ω

Z

∂GΩ (x; y) dσx = ∂nx

G(x; y) S (y)

∂u (x) dσx − ∂n

Z vΩ (x; y) S (y)

Z

∂u (x) dσx + ∂n

u(x) S (y)

∂G (x; y) dσx − ∂nx

Z u(x) S (y)

∂vΩ (x; y) dσx . ∂nx

The integrands for the last two integrals are continuous, with no singularities as → 0+ . Therefore, these two integrals tend to zero as → 0+ . We have G(x; y) = −

1 ln (kx − yk) , 2π

and therefore G(x; y) is constant for x ∈ S (y). It follows that Z ∂u G(x; y) (x) dσx ∂n S (y) is just a multiple of Z S (y)

∂u (x) dσx = ∂n

Z ∆u(x) dx = 0 B (y)

(applying the divergence theorem and using the fact that ∆u = 0 in Ω). It remains only to calculate Z ∂G u(x) (x; y) dσx . ∂n x S (y) Now, ∇x G(x; y) = − and, on S (y), nx = −

x−y 2πkx − yk2

x−y kx − yk

(notice that this is the outward-pointing normal to Ω ). Therefore, ∂GΩ 1 (x; y) = ∇x G(x; y) · nx = . ∂nx 2π

106

CHAPTER 11. PROBLEMS IN MULTIPLE SPATIAL DIMENSIONS Thus

Z u(x) S (y)

∂G 1 (x; y) dσx = ∂nx 2π

Z u(x) dσx , S (y)

which is the average value of u on S (y), and therefore Z ∂G u(x) (x; y) dσx → u(y) ∂nx S (y) as → 0+ . Therefore, it finally follows, by taking the limit as → 0+ in Green’s second identity, that Z ∂GΩ φ(x) (x; y) dσx = −u(y), ∂nx ∂Ω which, given the symmetry of GΩ , yields the desired result. 11. In the text it is shown that, if u solves −∆u = 0 in BR (0), u = φ(x) on SR (0), then Z u(x) = SR (0)

R2 − kxk2 φ(y) dσy . 2πRkx − yk2

Now suppose that u is to represented as a function of polar coordinates (r, θ). Let the circle SR (0) be represented by y = (R cos (ψ), R sin (ψ)), 0 ≤ ψ ≤ 2π. Then s 2 2 q ∂y1 ∂y2 dσy = + dψ = R2 cos2 (ψ) + R2 sin2 (ψ) dψ = Rdψ. ∂ψ ∂ψ With x represented by the polar coordinates (r, θ), we have kxk2 = r2 and kx − yk2 = (r cos (θ) − R cos (ψ))2 + (r sin (θ) − R sin (ψ))2 = r2 + R2 − 2rR cos (ψ − θ) (where we used the trignometric identity cos (θ) cos (ψ) + sin (θ) sin (ψ) = cos (ψ − θ)). Changing variables from y to (r, θ) thus yields Z Z 2π R2 − kxk2 R2 − r 2 1 φ(ψ) dψ, φ(y) dσy = 2 2 2 2π 0 R + r − 2Rr cos (ψ − θ) SR (0) 2πRkx − yk as desired. 13. In Chapter 9, it was shown that the solution to du d − P (x) + R(x)u = 0, a < x < b, dx dx u(a) = ua , u(b) = ub is

∂G ∂G (y; a) − P (b)ub (y; b), ∂y ∂y where G is the Green’s function for the inhomogeneous differential equation with homogeneous boundary conditions. If we take P (x) = 1, we obtain ∂G ∂G u(y) = − − (y; a)ua + (y; b)ub . ∂y ∂y u(y) = P (a)ua

Thus u is obtained by “integrating” (that is, summing) the normal derivative of the Green’s function, times the boundary “function” over the boundary of Ω = (a, b) (that is, the set {a, b}). Notice that the unit normal to Ω at x = b is n = 1, while at x = a it is n = −1, so that −

∂G (y; a) ∂y

11.7. GREEN’S FUNCTION FOR THE WAVE EQUATION

107

is the normal derivative of G at x = a, while ∂G (y; b) ∂y is the normal derivative of G at x = b. Thus formula (9.24) is exactly analogous to formula (11.103). 15. If Ω = BR (z) for some z ∈ R2 , then formula (11.107) can be modified to show that the solution of −∆u = f (x), x ∈ BR (z), u(x) = 0, x ∈ SR (z) is Z u(x) = SR (z)

R2 − kx − zk2 φ(y) dσy . 2πRkx − yk2

Substituting x = z, we see that kz − yk2 = R2 (since y ∈ SR (z), and of course kz − zk2 = 0; therefore, Z Z R2 1 u(z) = φ(y) dσy = φ(y) dσy . 3 2πR SR (z) SR (z) 2πR Thus u(z) is the average value of u over the circle SR (z) (recall that u(y) = φ(y) on SR (z)).

11.7

Green’s function for the wave equation

1. Let ψ : R3 → R, c ∈ R, and y ∈ R3 be given. Then, by the chain rule, 3

X ∂ψ ∂ ∂ (ψ(x + cty)) = (x + cty) (xi + ctyi ) ∂t ∂x ∂t i i=1 =

3 X ∂ψ (x + cty)(cyi ) ∂x i i=1

=c

3 X ∂ψ (x + cty)yi = ∇ψ(x + cty) · y. ∂x i i=1

Similarly, ! 3 X ∂ψ c (x + cty)yi ∂xi i=1 3 X ∂ ∂ψ =c (x + cty)yi ∂t ∂xi i=1

∂2 ∂ (ψ(x + cty)) = ∂t2 ∂t

=c

3 X 3 X ∂2ψ ∂ (x + cty)yi (xj + ctyj ) ∂x ∂t j ∂xi i=1 j=1

=c

3 X 3 X ∂2ψ (x + cty)yi cyj ∂xj ∂xi i=1 j=1

2

=c

3 X i=1

= c2

3 X

yi

3 X ∂2ψ (x + cty)yj ∂xj ∂xi j=1

yi ∇2 ψ(x + cty)y

i

i=1

= c2 y · ∇2 ψ(x + cty)y. 3. We wish to derive the solution of ∂2u − c2 ∆u = 0, x ∈ R2 , t > 0, ∂t2 u(x, 0) = 0, x ∈ R2 , ∂u (x, 0) = ψ(x), x ∈ R2 . ∂t

!

108

CHAPTER 11. PROBLEMS IN MULTIPLE SPATIAL DIMENSIONS We can follow the same reasoning as in 11.7.2. The solution of the same IVP in three dimensions is Z t u(x, t) = ψ(x + cty) dσy . 4π ∂B1 (0) If, in this formula, ψ does not depend on x3 , then neither does u, and it is easy to see that u, regarded as a function on R2 , is the solution of the IVP given above. Therefore, it remains only to write u as an integral over a two-dimensional region, namely, B1 (0). If S represents the upper hemisphere of ∂B1 (0), then Z Z ψ(x + cty) dσy ψ(x + cty) dσy = 2 S

∂B1 (0)

1

Z

Z √1−x2

=2

ψ(x1 + cty1 , x2 + cty2 ) p dy2 dy1 1 − y12 − y22 ψ(x + cty) p dy. 1 − kyk2 √

−

−1

Z =2 B1 (0)

It follows that

1−x2

Z

ψ(x + cty) p dy. 1 − kyk2 B1 (0) Performing the change of variables z = x + cty yields the alternate formula Z ψ(z) 1 p dy. u(x, t) = 2πc Bct (x) c2 t2 − kz − xk2 u(x, t) =

t 2π

5. The solution of ∂2u − c2 ∆u = 0, x ∈ R3 , t > 0, ∂t2 u(x, 0) = 0, x ∈ R3 , ∂u (x, 0) = ψ(x), x ∈ R3 . ∂t is u(x, t) =

t 4π(ct)2

Let

ψ(x) =

1, 0,

Z ψ(z) dσz . ∂Bct (x)

kxk < 1, otherwise,

and let u be the corresponding solution of the above IVP. Suppose first that x ∈ R3 , kxk > R, and t < (kvk−R)/c. This last condition implies that ct < kxk − R, and therefore, if z ∈ Sct (x), then kz − xk = ct < kxk − R. It is easy to see that this condition implies that kzk > R (since z is closer to x than x is to SR (0)). Therefore, z ∈ Sct (x) ⇒ ψ(z) = 0, and it follows from the formula for u(x, t) that u(x, t) = 0. Now suppose t > (kxk + R)/c, that is, ct > kxk + R. Then z ∈ Sct (x) ⇒ kz − xk = ct > kxk + R, which implies that kzk > R (since kz − xk ≤ kzk + kxk). Therefore, in this case also, z ∈ Sct (x) ⇒ ψ(z) = 0, and we see that u(x, t) = 0. 7. The causal Green’s function is G(x, t; y, s) =

( P ∞

n=1

sin (c

√ λn (t−s)) √ ψn (x)ψn (y), c λn

0,

t > s, t < s.

11.8. GREEN’S FUNCTIONS FOR THE HEAT EQUATION

11.8

109

Green’s functions for the heat equation

1. We wish to show that the Gaussian kernel, 2

S(x, t) = (4kπt)−n/2 e−kxk

/(4kt)

satisfies the homogeneous heat equation in Rn . We have 2 kxk2 ∂S n (x, t) = (4kπt)−n/2 e−kxk /(4kt) − , ∂t 4kt2 2t 2 xi ∂S (x, t) = (4kπt)−n/2 e−kxk /(4kt) − , ∂xi 2kt ∂2S x2i 1 −n/2 −kxk2 /(4kt) (x, t) = (4kπt) e − , ∂x2i 4k2 t2 2kt 2 kxk2 n ∆S(x, t) = (4kπt)−n/2 e−kxk /(4kt) − . 4k2 t2 2kt From these formulas, it follows immediately that ∂S (x, t) − k∆S(x, t) = 0. ∂t 3. We have seen that the IBVP ∂v − k∆v = f (x, T − t), x ∈ Ω, 0 < t < T, ∂t v(x, 0) = φ(x), x ∈ Ω, v(x, t) = 0, x ∈ ∂Ω, 0 < t < T has a unique solution v. Let us define u(x, t) = v(x, T − t). Then ∂u ∂v (x, t) = − (x, T − t), ∂t ∂t ∆u(x, t) = ∆v(x, T − t), and hence −

∂u ∂v (x, t) − k∆u(x, t) = (x, T − t) − k∆v(x, T − t) = f (x, T − (T − t)) = f (x, t). ∂t ∂t

Also, u(x, T ) = v(x, 0) = φ(x) for all x ∈ Ω, and u(x, t) = v(x, T − t) = 0 for all x ∈ ∂Ω and 0 < t < T . Thus u satisfies −

∂u − k∆u = f (x, t), x ∈ Ω, 0 < t < T, ∂t u(x, T ) = φ(x), x ∈ Ω, u(x, t) = 0, x ∈ ∂Ω, 0 < t < T,

which shows that this IBVP has a solution. 5. Let u, v be smooth functions defined for x ∈ Ω and 0 ≤ t ≤ T . We have Z TZ ∂v ∂u u − k∆v − v − k∆u dx dt ∂t ∂t 0 Ω Z Z T Z Z T Z TZ Z TZ ∂v ∂u = u dt dx − k u∆v dx dt + v dt dx + k v∆u dx dt ∂t ∂t Ω 0 0 Ω Ω 0 0 Ω Z Z T Z T Z Z Z Z T ∂u ∂v ∂u dt dx − k u dσx − ∇u · ∇v dx dt + v dx dt+ = uv|T0 − v ∂t ∂n ∂t Ω 0 0 ∂Ω Ω Ω 0 Z T Z Z ∂u k v dσx − ∇u · ∇v dx dt 0 ∂Ω ∂n Ω Z Z TZ ∂u ∂v = (u(x, T )v(x, T ) − u(x, 0)v(x, 0)) dx + k v −u dσx . ∂n ∂n Ω 0 ∂Ω

Chapter 12

More about Fourier Series 12.1

The complex Fourier series

1. The complex Fourier series of f is ∞ X

cn eiπnx ,

n=−∞

where c0 = 2/3 and cn = −

2(−1)n , n = ±1, ±2, . . . . n2 π 2

The errors in approximating f by a partial Fourier series are shown in Figure 12.1. 0.04 0.02 0 −0.02 −0.04 −1 0.04

−0.5

0 x

0.5

1

−0.5

0 x

0.5

1

−0.5

0 x

0.5

1

0.02 0 −0.02 −0.04 −1 0.04 0.02 0 −0.02 −0.04 −1

Figure 12.1: The error in approximating f (x) = 1 − x2 by its partial complex Fourier series with N = 10 (top), N = 20 (middle), and N = 40 (bottom). (See Exercise 12.1.1.)

3. The complex Fourier coefficients of f are cn =

(−1)n+1 sin (1) . nπ − 1

The magnitudes of the error in approximating f by a partial Fourier series are shown in Figure 12.2. 5. Let g(x) = a0 +

∞ nπx nπx X an cos + bn sin , ` ` n=1

111

112

CHAPTER 12. MORE ABOUT FOURIER SERIES 1

0.5

0 −1 1

−0.5

0 x

0.5

1

−0.5

0 x

0.5

1

−0.5

0 x

0.5

1

0.5

0 −1 1

0.5

0 −1

Figure 12.2: The magnitude of the error in approximating f (x) = eix by its partial complex Fourier series with N = 10 (top), N = 20 (middle), and N = 40 (bottom). (See Exercise 12.1.3.) where a0

=

an

=

bn

=

Z ` 1 g(x) dx, 2` −` Z nπx 1 ` g(x) cos dx, n = 1, 2, 3, . . . , ` −` ` Z ` nπx 1 g(x) sin dx, n = 1, 2, 3, . . . , ` −` `

(see (6.25) in the text), and similarly let h(x) = p0 +

∞ nπx nπx X pn cos + qn sin . ` ` n=1

The complex Fourier coefficient of f is given by cn =

1 2`

`

Z

f (x)e−iπnx/` dx.

−`

For n > 0, we have cn

= =

=

Z ` nπx nπx 1 (g(x) + ih(x)) cos − i sin dx 2` −` ` ` Z Z nπx nπx 1 1 ` 1 ` g(x) cos dx − i g(x) sin dx 2 ` −` ` ` −` ` Z ` Z nπx nπx 1 1 ` +i h(x) cos dx + h(x) sin dx ` −` ` ` −` ` 1 (an + qn + i(pn − bn )) . 2

Similarly, if n > 0, then c−n =

1 (an − qn + i(pn + bn )) . 2

Finally, c0 = a0 + ip0 . 7. We have N −1 X

eπijn/N

2

=

n=0

N −1 X

e2πij/N

n=0 2πij

=

e −1 , e2πij/N − 1

n

12.2. FOURIER SERIES AND THE FFT

113

using the formula for the sum of a finite geometric series: m X

rn =

n=0

rm+1 − 1 . r−1

Moreover, since eiθ is a 2π-periodic function of θ, we see that e2πij = 1 and hence that N −1 X

eπijn/N

2

= 0.

n=0

On the other hand, by Euler’s formula, N −1 X

eiπjn/N

2

N −1 X

2 jnπ jnπ + i sin N N n=0 N −1 X jnπ jnπ jnπ jnπ − sin2 + 2i cos sin cos2 N N N N n=0 N −1 N −1 X X jnπ jnπ − cos2 sin2 N N n=0 n=0 N −1 X jnπ jnπ +2i sin . cos N N n=0

=

n=0

=

=



cos

Since this expression equals zero, we must have N −1 X

cos2



n=0

jnπ N

N −1 X

−

N −1 X

sin2



n=0

cos

n=0

jnπ N



sin

jnπ N



jnπ N



=

0,

=

0.

The second equation is one of the results we set out to prove. The first equation can be written, using the trigonometric identity sin2 (θ) = 1 − cos2 (θ), as N −1 X

cos2

n=0



jnπ N

−

N −1 X

1 − cos2



n=0

jnπ N

which simplifies to 2

N −1 X

cos2



n=0

jnπ N

= N,

or N −1 X

cos2



n=0

jnπ N



jnπ N



=

N , 2

=

N 2

which is the desired result. The result N −1 X

sin2

n=0



then follows.

12.2

Fourier series and the FFT

1. The graph of {|Fn |} is shown in Figure 12.3. 3. As shown in the text,

. cn = Fn , n = −16, −15, . . . , 15,

= 0,

114

CHAPTER 12. MORE ABOUT FOURIER SERIES 0.5

0.4

|Fn|

0.3

0.2

0.1

0 0

5

10 n

15

20

Figure 12.3: The magnitude of the sequence {Fn } from Exercise 12.2.1. 0.5 0.4 0.3 0.2 0.1 0 −0.1 −0.2 −0.3 −0.4 −0.5 −1

−0.5

0 x

0.5

1

Figure 12.4: The function f (x) = x(1 − x2 ) (the curve) and the estimates of f (xj ) (the circles) computed from the complex Fourier coefficients of f (see Exercise 12.2.3.

15 where {Fn }15 n=−16 is the DFT of {fj }j=−16 ,

fj = f (xj ), xj = −1 +

j , j = −16, −15, . . . , 15 16

(note that f (−1) = f (1)). Therefore, we can estimate {f (xj )} by taking the inverse DFT (using the inverse FFT) of {cn }. The results are shown in Figure 12.4. 5. We have N −1 X

An e2πijn/N

=

n=0

=

=

N −1 N −1 1 X X am e−2πimn/N e2πijn/N N n=0 m=0 N −1 N −1 1 X X am e2πi(j−m)n/N N n=0 m=0 (N −1 ) N −1 X 2πi(j−m)n/N 1 X am e . N m=0 n=0

If j = m, then N −1 X n=0

e2πi(j−m)n/N =

N −1 X n=0

1 = N,

12.3. RELATIONSHIP OF SINE AND COSINE SERIES TO THE FULL FOURIER SERIES

115

while if j 6= m, then this sum is a finite geometric series: N −1 X

e2πi(j−m)n/N

N −1 X

=

n=0

e2πi(j−m)/N

n

n=0

=

e2πi(j−m)/N

N

−1

−1 e2πi(j−m) − 1 e2πi(j−m)/N − 1 0. e2πi(j−m)/N

= = Therefore, am

N −1 X

e2πi(j−m)n/N =



N aj , 0,

n=0

m = j, m 6= j,

and so we obtain N −1 X

2πijn/N

An e

n=0

N −1 1 X am = N m=0

(N −1 X

) 2πi(j−m)n/N

e

=

n=0

1 N aj = aj , N

as desired. 7. Below we show the exact Fourier sine coefficients of f , the coefficients estimated by the DST, and the relative error. Since both an and the estimated an are zero for n even, we show only an for n odd. n

an −1

estimated an

relative error

−1

1

2.5801 · 10

2.5801 · 10

6.2121 · 10−6

3

9.5560 · 10−3

9.5511 · 10−3

5.1571 · 10−4

5

−3

2.0641 · 10

−3

2.0555 · 10

4.1824 · 10−3

7

7.5222 · 10−4

7.3918 · 10−4

1.7340 · 10−2

9

−4

3.5393 · 10

−4

3.3531 · 10

5.2608 · 10−2

11

1.9385 · 10−4

1.6778 · 10−4

1.3448 · 10−1

13

1.1744 · 10−4

8.0874 · 10−5

3.1135 · 10−1

15

−5

−5

6.8241 · 10−1

7.6448 · 10

2.4279 · 10

9. We have 2

N −1 X

Fn sin

n=1



jnπ N

=

4

N −1 N −1 X X

fm sin

mnπ N

n=1 m=1

=

4

N −1 X

fm

m=1

N −1 X

sin

mnπ N

n=1



jnπ N



jnπ . N

sin

sin



By the hint, N −1 X

sin

mnπ

n=1

N

sin

jnπ N



is either N/2 or 0, depending on whether m = j or not. It then follows immediately that

2

N −1 X n=1

which is what we wanted to prove.

Fn sin

jnπ N

=4

N fj = 2N fj , 2

116

CHAPTER 12. MORE ABOUT FOURIER SERIES 1.5

1

0.5

0

−0.5

−1

−1.5 −3

−2

−1

0 x

1

2

3

Figure 12.5: The partial Fourier sine series (50 ) of f (x) = x.

12.3

Relationship of sine and cosine series to the full Fourier series

1. The partial Fourier sine series, with 50 , is shown in Figure 12.5. It appears that the series converges to the function F satisfying F (x) = x − 2k, −2k − 1 < x < 2k + 1, k = 0, ±1, ±2, . . . . The period of this function is 2. 3.

(a) The odd extension fodd is continuous on [−`, `] only if f (0) = 0. Otherwise, fodd has a jump discontinuity at x = 0.

(b) The even extension feven is continuous on [−`, `] for every continuous f : [0, `] → R. 5. The quarter-wave sine series of f : [0, `] → R is the full Fourier series of the function f˜ : [−2`, 2`] → R obtained by first reflecting the graph of f across the line x = ` (to obtain a function defined on [0, 2`]) and then taking the odd extension of the result. That is, f˜ is defined by  f (x), 0 ≤ x ≤ `,    f (2` − x), ` < x ≤ 2`, f˜(x) = −f (−x), −` ≤ x < 0,    −f (x + 2`), −2` ≤ x < −`. Since this function is odd, its full Fourier series has only sine (all of the cosine coefficients are zero), and because of the other symmetry, the even sine also drop out.

12.4

Pointwise convergence of Fourier series

1. The Fourier series of f converges pointwise to the function F defined by 0, x = ±(2k − 1)π, k = 1, 2, 3, . . . , F (x) = (x − 2kπ)3 , (2k − 1)π < x < (2k + 1)π, k = 0, ±1, ±2, . . . . 3. There are many such functions f , but they all have one or more discontinuities. An example is x, 0 ≤ x ≤ 12 , f (x) = x + 1, 21 < x ≤ 1. 5. If {fN } converges uniformly to f on [a, b], then it obviously converges pointwise (after all, the maximum difference, over x ∈ [a, b], between fN (x) and f (x) converges to zero, so each individual difference must converge to zero). To prove that uniform convergence implies mean-square convergence, define MN = max{|f (x) − fN (x)| : a ≤ x ≤ b}. Then s

Z

s Z

b

|f (x) − fN (x)|2 dx

b 2 MN dx

≤

a

a

=

√ MN b − a.

By hypothesis, MN → 0 as N → ∞, which gives the desired result.

12.5. UNIFORM CONVERGENCE OF FOURIER SERIES

117

7. Suppose f : [0, `] → R is continuous. Then feven , the even, periodic extension of f to R, is defined by f (x − 2k`), 2k` ≤ x < (2k + 1)`, k = 0, ±1, ±2, . . . , feven (x) = f (−x + 2k`), (2k − 1)` ≤ x < 2k`, k = 0, ±1, ±2, . . . . Obviously, then, feven is continuous on every interval (2k`, (2k + 1)`) and ((2k − 1)`, 2k`), that is, except possibly at the points 2k` and (2k − 1)`, k = 0, ±1, ±2, . . .. We have lim

feven (x) =

lim

feven (x) =

x→2k`−

lim

f (−x + 2k`) = lim f (x) = f (0)

lim

f (x − 2k`) = lim f (x) = f (0).

x→2k`−

x→0+

and x→2k`+

x→2k`+

x→0+

Since feven (2k`) = f (0) by definition, this shows that feven is continuous at x = 2κ`. A similar calculation shows that lim

x→(2k−1)`−

feven (x) =

lim

x→(2k−1)`+

feven (x) = feven ((2k − 1)`) = f (`),

and hence that feven is continuous at x = (2k − 1)`. 9. Given f : [0, `] → R, define f˜ : [−2`, 2`] → R by  f (x),    f (2` − x), ˜ f (x) = −f (−x),    −f (x + 2`),

0 ≤ x ≤ `, ` < x ≤ 2`, −` ≤ x < 0, −2` ≤ x < −`.

Then define F : R → R to be the periodic extension of f˜ to R. Assuming f is piecewise smooth, the quarter-wave sine series of f converges to F (x) if F is continuous at x, and to 1 [F (x−) + F (x+)] 2 if F has a jump discontinuity at x. If f is continuous, then its quarter-wave sine series converges to f except possibly at x = 0. If f is continuous and f (0) = 0, then the quarter-wave sine series of f converges to f at every x ∈ [0, `].

12.5

Uniform convergence of Fourier series

1. The function h(x) fails to satisfy h(−1) = h(1), so its Fourier coefficients decay like 1/n and its Fourier series is the slowest to converge. The function f satisfies f (−1) = f (1), but df /dx has a discontinuity at x = 0 (and the derivative of fper also has discontinuities at x = ±1). Therefore, the Fourier coefficients of f decay like 1/n2 . Finally, gper and its first derivative are continuous, but its second derivative has a jump discontinuity at x = ±1, so its Fourier coefficients decay like 1/n3 . The Fourier series of g is the fastest to converge. 3. 5. Figure 12.6 shows the f , its partial Fourier series with 21, 41, 81, and 161 , and the line y = 1.09. Zooming in near x = 0 shows that the overshoot is indeed about 9%; see Figure 12.7.

12.6 1.

Mean-square convergence of Fourier series (a) The infinite series ∞ X 1 2 n n=1

converges to a finite value.1 Therefore, {1/n} ∈ `2 and so, by Theorem 12.36, the sine series converges to a function in L2 (0, 1). 1A

standard result in the theory of infinite series is that ∞ X 1 k n n=1

118

CHAPTER 12. MORE ABOUT FOURIER SERIES 1.5

y

1

0.5

0 21 41 81 161 −0.5 −1

−0.5

0 x

0.5

1

Figure 12.6: The function f from Exercise 12.5.5, its partial Fourier series with 21, 41, 81, and 161 , and the line y = 1.09. 1.105 21 41 81 161

1.1

y

1.095

1.09

1.085

1.08

1.075 −0.02

0

0.02

0.04 x

0.06

0.08

0.1

Figure 12.7: Zooming in on the overshoot in Figure 12.6. (b) Figure 12.8 shows the sum of the first 100 of the sine series. The graph suggests that the limit f is of the form f (x) = m(1 − x). The Fourier sine series of such an f are Z 1 2m cn = 2 m(1 − x) sin (nπx) = , n = 1, 2, 3, . . . . nπ 0 Therefore, in order that cn = 1/n, we must have m = π/2. Therefore, f (x) =

π (1 − x). 2

2

1.5

1

0.5

0

−0.5 0

0.2

0.4

0.6

0.8

1

x

Figure 12.8: The partial sine series, with 100 , from Exercise 12.6.1. converges if k is greater than 1. This can be proved, for example, by comparison with the improper integral Z ∞ dx . xk 1

12.7. A NOTE ABOUT GENERAL EIGENVALUE PROBLEMS 3. The infinite series

119

2 X ∞ ∞ X 1 1 √ = n n n=1 n=1

√ is the harmonic series, a standard example of a divergent series. Therefore, {1/ n} 6∈ `2 , and so the series does 2 not converge to an L (0, 1) function. R1 5. (a) The function f (x) = xk belongs to L2 (0, 1) if and only if the integral 0 x2k dx is finite. Provided k 6= −1/2, we have 1 Z 1 Z 1 , k > − 21 , 1 − r2k+1 1+2k x2k dx = lim x2k dx = lim = ∞, k < − 12 . 1 + 2k r→0+ r→0+ r 0 For k = −1/2, Z 1 Z 1 Z 1 dx dx x2k dx = = lim = lim − ln (r) = ∞. + x x r→0 r→0+ 0 0 r 2 Thus we see that f ∈ L (0, 1) if and only if k > −1/2. (b) The first few Fourier sine coefficients of f (x) = x−1/4 , as computed by the DST, are approximately 1.5816, 0.25353, 0.63033, 0.17859, 0.41060, 0.14157, . . . . (c) The graphs of f , the partial Fourier sine series, with 63 , and the difference between the two, are shown in Figure 12.9. 6

y=x−1/4 sine series (63 )

y

4

2

0 0

0.2

0.4

0.6

0.8

1

0.6

0.8

1

x 0.1

0

−0.1 0

0.2

0.4 x

Figure 12.9: The function f (x) = x−1/4 , together with its partial Fourier sine series (first 63 ) (top), and the difference between the two (bottom). See Exercise 12.6.5. 7. Since vn → v, there exists a positive integer N such that n ≥ N ⇒ kv − vn k <

. 2

Then, if n, m ≥ N , we have kvn − vm k ≤ kvn − vk + kv − vm k <

+ = . 2 2

Therefore, {vn } is Cauchy.

12.7

A note about general eigenvalue problems

2 1. The proof is a direct calculation, using integration by parts (Green’s identity): Suppose u, v ∈ CD (Ω). Then Z Z Z ∂u (KD u, v) = − ∇ · (k(x)∇u) v = k(x)∇u · ∇v − k(x)v ∂n Ω ∂Ω ZΩ = k(x)∇u · ∇v (since v = 0 on ∂Ω) Ω Z Z ∂v = − ∇ · (k(x)∇v) u + k(x)u ∂n ∂Ω ZΩ = − ∇ · (k(x)∇v) u (since u = 0 on ∂Ω) Ω

=

(u, KD v).

120

CHAPTER 12. MORE ABOUT FOURIER SERIES

3. If t is very large, then we can approximate u(x, t) by the first term in its generalized Fourier series. Using the same notation as before for the eigenvalues and eigenfunctions of the negative Laplacian on Ω, we have . u(x, t) = c1 e−κλ1 t/(ρc) ψ1 (x), t large. The constant c1 is the first generalized Fourier coefficient of the initial temperature 5: R 5ψ1 c1 = RΩ 2 . ψ Ω 1 The approximation is valid provided λ1 is a simple eigenvalue; that is, there is only one linearly independent eigenvector corresponding to λ1 . Then all of the other in the generalized Fourier series decay to zero much more rapidly than does the first term.

Chapter 13

More about Finite Element Methods 13.1

Implementation of finite element methods

1. We must check that the one-point rule gives the exact values for the integral f (x) = x2 . With f (x) = 1, we have Z 1 f = (area of T R ) · 1 = 2 TR

R

TR

f when f (x) = 1, f (x) = x1 , or

and 1 f 2



1 1 , 3 3



1 . 2

=

With f (x) = x1 , we have Z

1

Z f

1−x1

Z

=

x1 dx2 dx1

TR

0

0 1

Z = 0

1 1 1 x1 − x21 dx1 = − = . 2 3 6

On the other hand,

1 f 2

1 1 , 3 3

=

1 1 1 · = , 2 3 6

so the rule is exact in this case as well. By symmetry, the rule must hold for f (x) = x2 , which shows that the rule has degree of precision at least 1. To show that the degree of precision is not greater than 1, we note that, for f (x) = x21 , Z Z 1 Z 1−x1 Z 1 1 1 1 f= x21 dx2 dx1 = x21 − x31 dx1 = − = , 3 4 12 R 0 0 0 T while 1 f 2



1 1 , 3 3

=

1 1 1 · = . 2 9 18

3. We have Z

3/2

Z

2x1 −2

Z

f= T

1

x1 x22

2

Z

4−2x1

Z

dx2 dx1 +

0

3/2

x1 x22 dx2 dx1 =

0

On the other hand, the mapping from T R to T is x1

=

1 + z1 +

x2

=

z2 .

1 z2 , 2

The Jacobian is J=

121

1 0

1 2

1

,

7 1 1 + = . 120 15 8

122

CHAPTER 13. MORE ABOUT FINITE ELEMENT METHODS and its determinant is 1. Therefore, with g(z1 , z2 ) = f

1 1 1 + z1 + z2 , z2 = 1 + z1 + z2 z22 , 2 2

we have

1−z1

Z



g=

f= T

1

Z

Z

Z

TR

0

0

1 1 + z1 + z2 z22 dz2 dz1 . 2

This last integral also equals 1/8. 5. Recall that the nodes in the mesh are enumerated from 1 to M . Let the constrained boundary nodes be those with numbers c1 , c2 , . . . , cK . We define a new K × 1 array CNodePtrs whose ith entry is ci . (Thus CNodePtrs contains pointers into the NodeList array, allowing one to find the coordinates of any given constrained node.) We need the inverse of the mapping i 7→ ci ; let us define j 7→ dj by i if j = ci , dj = 0 if j 6= ci for all i = 1, 2, . . . , K. Now, under the original description of NodePtrs, the ith entry is zero if the ith node is constrained. We now redefine the ith entry to be −di if the ith node, while the ith entry is still gi if the ith node is free. (The use of the negative sign is just a trick to allow us to distinguish a free node from a constrained node; it avoids the need to store a separate flag.) Suppose we wish to solve a BVP with inhomogenous Dirichlet conditions (on all or part of the boundary). Let G be the piecewise linear function whose nodal value is zero everywhere except at the constrained boundary nodes. At the nonfree boundary nodes, the value of G is the assigned Dirichlet data. To solve the inhomogeneous Dirichlet problem, we merely need to modify the load vector by replacing Z Fi = f φi Ω

by Z f φi − a(G, φi ).

Fi = Ω

The quantity a(G, φi ) is nonzero only if the free node nfi belongs to a triangle that also contains a constrained boundary node. The point here is simply that this can easily be determined from the information in the data structure. As we loop over the triangles, we can determine, for each, whether it contains both a free node and a constrained node. If it does, we modify the load vector accordingly.

13.2

Solving sparse linear systems

1. Computing A = LU requires n − 1 steps, and step i uses (n − i)(1 + 2(n − i)) = 2(n − i)2 + (n − i) arithmetic operations (this is determined by simply counting the operations in the pseudo-code on page 604). The total number of operations required is n−1 X i=1



n−1 X 2 n−1 X n2 n 2n3 2(n − i)2 + (n − i) = 2 j + − − . j= 3 2 6 j=1 j=1

This agrees with the O(2n3 /3) operation count given in the text. The computation of L−1 b requires n − 1 steps, with 2(n − i) arithmetic operations per step. The total is 2

n−1 X i=1

(n − 1) = 2

n−1 X

j = n2 − n.

j=1

The final step of back substitution requires n steps, with 2(n − i) + 1 operations per step, for a total of n X

{2(n − i) + 1} = 2

i=1

n−1 X

j + n = n2 .

j=1 2

The total for the two triangular solves is 2n − n, which also agrees with the count given in the text.

13.2. SOLVING SPARSE LINEAR SYSTEMS

123

3. With φ(x) =

1 x · Ax − b · x, 2

we have φ(x + y) =

1 (x + y) · A(x + y) − b · (x + y) 2

1 1 x · Ax + y · Ax + y · Ay − b · x − b · y 2 2 1 φ(x) + (Ax − b) · y + y · Ay. 2

= =

The fact that A is symmetric allows the simplification 1 1 x · Ay + y · Ax = y · Ax, 2 2 which was used above. We thus see that φ(x + y) = φ(x) + ∇φ(x) · y + O kyk2 , with ∇φ(x) = Ax − b. 5. An inner product has three properties: (a) (x, x) ≥ 0 for all x, and (x, x) = 0 if and only if x = 0. (b) (x, y) = (y, x) for all x, y. (c) (αx + βy, z) = α(x, z) + β(y, z) for all x, y, z and all α, β. These properties hold for (x, y)A = x · Ay. The third property would be true for any matrix A, the second property obviously requires that A be symmetric, and the first property requires that A be positive definite. 7. (a) With x(0) = (4, 2), the negative gradient is −∇φ x(0) = (−2, −1), and so we must minimize φ x(0) − α∇φ x(0) = 7α2 − 5α − 18. The minimizer is α = 5/14, and so x(1) =



23 23 , 7 14

.

(b) At x(1) = (23/7, 23/14), the negative gradient is (−3/14, 3/7), so we must minimize 45 529 27 2 φ x(1) − α∇φ x(1) = α − α− . 196 196 28 The minimizer is α = 5/6, and so x(2) =



87 ,2 . 28

The desired solution is (3, 2). (c) With x(0) = 0, the CG algorithm produces . x(1) = (2.67456, 2.34024), with x(2) = (3, 2), the exact solution. 9. If y = x − x(0) , then x = y + x(0) and so Ax = b ⇒ Ay + Ax(0) = b ⇒ Ay = b − Ax(0) . We can therefore replace b by b − Ax(0) , apply the CG algorithm to estimate y, and add x(0) to get the estimate of x.

124

CHAPTER 13. MORE ABOUT FINITE ELEMENT METHODS

13.3

An outline of the convergence theory for finite element methods

1. Let the vector-valued function F be defined by F(x) = (f (x), 0). Then ∇ · (φF) = φ∇ · F + ∇φ · F = φ

∂f ∂φ + f, ∂x1 ∂x1

so, by the divergence theorem, Z Z Z Z Z ∂f ∂φ φ + f= ∇ · (φF) = φF · n = φf n1 . Ω ∂x1 Ω ∂x1 Ω ∂Ω ∂Ω If φ ∈ C0∞ (Ω), then the boundary integral vanishes, and we obtain Z Z ∂φ ∂f φ =− f. Ω ∂x1 Ω ∂x1 The derivation for ∂/∂x2 is exactly analogous. 3. A direct computation shows that kf − gkL2 =

1 for all integers m, n, 2

while

1p 1 + m2 π 2 + n2 π 2 for all m, n. 2 Thus the L2 error is independent of m, n, while the H 1 error becomes arbitrarily large as m, n → ∞. 5. The following table shows the L2 and energy norm errors. The initial mesh has eight triangles, and the side lengths of the triangles are divided by two each time the mesh is refined. kf − gkH 1 =

√h √2/2 √2/4 √ 2/8 √2/16 √2/32 2/64

L2 error 1.7803 · 10−2 5.6705 · 10−3 1.5112 · 10−3 3.8399 · 10−4 9.6392 · 10−5 2.4123 · 10−5

energy error 1.0190 · 10−1 5.3955 · 10−2 2.7318 · 10−2 1.3700 · 10−2 6.8554 · 10−3 3.4283 · 10−3

We see that the L2 error goes down by a factor of approximately four each time the mesh is refined, which is consistent with O(h2 ) convergence. Similarly, the energy norm of the error goes down by a factor of approximately two each time the mesh is refined, which is consistent with O(h) convergence.

13.4 1.

Finite element methods for eigenvalue problems (a) Let u, v ∈ Rn contain the nodal values of piecewise linear functions uh , vh , respectively, so that uh =

n X

ui φi , vh =

i=1

n X

v i φi .

i=1

It follows that Z (uh , vh ) =

uh vh = Ω

Z X n n X (ui φi ) (vi φi ) Ω i=1

=

j=1

n X n Z X Ω

i=1 j=1

=

n X n X

φi φi ui vi

Mij ui vi

i=1 j=1

=

u · Mv.

(b) As shown in the text, if u and v are generalized eigenvectors for the problem Ku = λMu, then LT u and LT v are orthogonal in the Euclidean norm. But then 0 = LT u · LT v = u · LLT v = u · Mv = (uh , vh ). Thus uh and vh are orthogonal in the L2 norm. 3. The smallest eigenvalue is approximately 9.80 (correct to three digits).