Contents 1 Classification of Differential Equations
1
2 Mo dels in one dimension
3
2.1 2.2 2.3 2.4
Heat Flow in a bar; Fourier’s law . . . . The hanging bar . . . . . . . . . . . . . The wave equation for a vibrating string Advection; kinematic waves . . . . . . .
. . . .
. . . .
. . . .
. . . .
. . . .
. . . .
. . . .
. . . .
. . . .
. . . .
. . . .
. . . .
. . . .
. . . .
. . . .
. . . .
. . . .
. . . .
. . . .
. . . .
. . . .
. . . .
. . . .
. . . .
. . . .
. . . .
. . . .
. . . .
. . . .
. . . .
. . . .
Linear systems as line inear oper perator equations . . . . . Existence Existence and and uniquene uniqueness ss of solutions solutions to to Ax = b . . Basis and dimension . . . . . . . . . . . . . . . . . . Orthogonal bases and proje ojections . . . . . . . . . . . Eige Eigen nvalue aluess and and eigen igenv vect ectors ors of a sym symmetri etricc ma matr trix ix
. . . . .
. . . . .
. . . . .
. . . . .
. . . . .
. . . . .
. . . . .
. . . . .
. . . . .
. . . . .
. . . . .
. . . . .
. . . . .
. . . . .
. . . . .
. . . . .
. . . . .
. . . . .
. . . . .
. . . . .
. . . . .
. . . . .
. . . . .
. . . . .
3 Essential Linear Algebra
3.1 3.2 3.3 3.4 3.5 3.5
7
4 Essential Ordinary Differential Equations
4.1 4.2 4.3 4.4 4.5
Background . . . . . . . . . . . . . . . . . . . Solutions to some simple ODEs . . . . . . . . Linear systems with constant coeffi oefficients . . . Numerica ical methods ods for init initia iall value problems Stiff systems of ODEs . . . . . . . . . . . . .
. . . . .
. . . . .
. . . . .
. . . . .
. . . . .
. . . . .
. . . . .
. . . . .
. . . . .
. . . . .
. . . . .
. . . . .
. . . . .
. . . . .
. . . . .
. . . . .
. . . . .
. . . . .
. . . . .
. . . . .
. . . . .
. . . . .
. . . . .
. . . . .
. . . . .
. . . . .
. . . . .
. . . . .
The The analo nalogy gy bet between BV BVPs Ps and and lin linea earr alg algeebra braic syst system emss Introd troduc ucti tion on to the spec spectr traal met meth hod; od; eig eigeenfun nfunct ctio ions ns . . . Solving the BVP using Fourier series . . . . . . . . . . . Finite element methods ods for BVPs . . . . . . . . . . . . . The Galerkin method . . . . . . . . . . . . . . . . . . . Piec Piecew ewis isee poly polyno nom mials ials and the finit finitee elem lement ent metho ethod d . .
. . . . . .
. . . . . .
. . . . . .
. . . . . .
. . . . . .
. . . . . .
. . . . . .
. . . . . .
. . . . . .
. . . . . .
. . . . . .
. . . . . .
. . . . . .
. . . . . .
. . . . . .
. . . . . .
. . . . . .
. . . . . .
. . . . . .
. . . . . .
. . . . . .
. . . . . .
25 27 29 32 33 35 39
Fourier series ies methods ods for the heat equation ion . . . . . . . Pure Pure Neum Neuman ann n condi onditi tion onss and and the the Fouri ourieer cosi cosine ne ser series ies . Period riodic ic boun bound dary ary cond condit itio ions ns and and the the full full Fourie urierr seri series es Finite element methods ods for the heat equation . . . . . . Finite elements and Neumann conditio ition ns . . . . . . . . .
. . . . .
. . . . .
. . . . .
. . . . .
. . . . .
. . . . .
. . . . .
. . . . .
. . . . .
. . . . .
. . . . .
. . . . .
. . . . .
. . . . .
. . . . .
. . . . .
. . . . .
. . . . .
. . . . .
. . . . .
. . . . .
. . . . .
. . . .
. . . .
. . . .
. . . .
. . . .
. . . .
. . . .
. . . .
. . . .
. . . .
. . . .
. . . .
. . . .
. . . .
. . . .
. . . .
. . . .
. . . .
. . . .
. . . .
. . . .
. . . .
7 Waves
7.1 7.1 7.2 7.3 7.4
15 17 19 21 23 25
6 Heat Flow and Diffusion
6.1 6.2 6.2 6.3 6.3 6.4 6.5
7 9 10 11 12 15
5 Boundary Value Problems in Statics
5.1 5.1 5.2 5.2 5.3 5.4 5.5 5.6 5.6
3 4 5 5
39 41 44 47 49 53
The The hom omog ogeeneou neouss wave equ equatio ation n with withoout boun bounda dari ries es Fourier series ies methods ods for the wave equation . . . . Finite element methods ods for the wave equation . . . . Resonance . . . . . . . . . . . . . . . . . . . . . . . . i
. . . .
. . . .
53 54 57 58
ii
7.5
CONTENTS
Finit initee diffe ifference methods ods for the wave equation . . . . . . . . . . . . . . . . . . . . . . . . . . .
8 First-order PDEs and the Method of Characteristics
8.1 8.1 8.2 8.3
63
The The sim simples plestt PDE PDE and and the meth ethod of chara aracte cteris ristics tics . . . . . . . . . . . . . . . . . . . . . . . . First-order quasilinear PDEs . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . Burgers’s equation . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . .
9 Green’s Functions
9.1 9.2 9.2 9.3 9.3 9.4 9.5 9.6
Green’s functions for BVPs in ODEs: Spe Special cases . . . . Gre Green’s en’s func unction tionss for BV BVPs Ps in ODEs ODEs:: the the sym symmetri etricc case ase Gre Green’s en’s func unction tionss for BV BVPs Ps in ODEs ODEs:: the the gene genera rall case case . . Introdu oduction ion to Green’s functions for IVPs . . . . . . . . . Green’s functions for the heat equation ion . . . . . . . . . . . Green’s functions for the wave equation . . . . . . . . . .
. . . . . .
. . . . . .
. . . . . .
. . . . . .
. . . . . .
. . . . . .
. . . . . .
. . . . . .
. . . . . .
. . . . . .
. . . . . .
. . . . . .
. . . . . .
. . . . . .
. . . . . .
. . . . . .
. . . . . .
. . . . . .
. . . . . .
. . . . . .
. . . . . .
Proper perties ies of the Sturm-Liou iouville lle oper perator . . . . . . . . Num Numeric ericaal Me Meth thod odss for Stur Sturm m-Lio Liouvil uville le pro proble blems . . . . Examples of Sturm-Liou iouville lle problems . . . . . . . . . . R ob obin boundary conditions . . . . . . . . . . . . . . . . Finit Finitee ele elem ment ent me method thodss for for Robi Robin n bou boun ndary dary cond condit itio ions ns . The The the theory ory of Stur Sturm m-Lio -Liouv uvil ille le prob proble lem ms: an an outl outlin inee . . .
. . . . . .
. . . . . .
. . . . . .
. . . . . .
. . . . . .
. . . . . .
. . . . . .
. . . . . .
. . . . . .
. . . . . .
. . . . . .
. . . . . .
. . . . . .
. . . . . .
. . . . . .
. . . . . .
. . . . . .
. . . . . .
. . . . . .
. . . . . .
. . . . . .
. . . . . .
Phys Physic icaal mod model elss in two or thr threee sp spatia atiall dim dimen ensi sioons . . . . . . Fourier ier series on a rectangular domain . . . . . . . . . . . . . F Foourier series on a disk . . . . . . . . . . . . . . . . . . . . . . Finite elements in two dimensions . . . . . . . . . . . . . . . . The The fr freeee-spac spacee G Gre reen en’s ’s fun functio ction n fo for the the La Lapl plaacian cian . . . . . . . The The Gre Green en’s ’s func functi tion on for for the the La Lapl plac acia ian n on on a boun bounde ded d dom domai ain n Green’s fu function for the wave equation . . . . . . . . . . . . . Green’s fu functions for the heat equation ion . . . . . . . . . . . . .
. . . . . . . .
. . . . . . . .
. . . . . . . .
. . . . . . . .
. . . . . . . .
. . . . . . . .
. . . . . . . .
. . . . . . . .
. . . . . . . .
. . . . . . . .
. . . . . . . .
. . . . . . . .
. . . . . . . .
. . . . . . . .
. . . . . . . .
. . . . . . . .
. . . . . . . .
. . . . . . . .
. . . . . . . .
The The complex Fourier series . . . . . . . . . . . F Foourier series and the FFT . . . . . . . . . . . Rela Relati tion onsh ship ip of sine sine and and cos cosin inee ser serie iess to to the the full full Poin Pointtwise wise con convergen rgence ce of Fouri ourier er ser series ies . . . . . Uniform convergence of Fourier series ies . . . . . Mean Mean--squa square re con converge ergenc ncee of Fourie urierr seri series es . . . A no note abou bout general ei eigenvalue problems . . .
89 90 94 96 100 100 10 1022 107 109 111
. . . . . . . . . . . . . . . . Fou Fouri rier er seri series es . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . .
. . . . . . .
. . . . . . .
. . . . . . .
. . . . . . .
. . . . . . .
. . . . . . .
. . . . . . .
. . . . . . .
. . . . . . .
. . . . . . .
. . . . . . .
. . . . . . .
. . . . . . .
. . . . . . .
. . . . . . .
. . . . . . .
. . . . . . .
. . . . . . .
. . . . . . .
13 More ab out Finite Element Methods
13.1 13.2 13.3 13.3 13.4 13.4
83 84 85 85 86 88 89
12 More ab out Fourier Series
12.1 12.2 12.3 12.3 12.4 12.4 12.5 12.6 12.6 12.7
71 72 75 76 79 80 83
11 Problems in Multiple Spatial Dimensions
11.1 11.1 11.2 11.3 11.4 11.5 11.5 11.6 11.6 11.7 11.8
63 65 67 71
10 Sturm-Liouville Eigenvalue Problems
10.2 10.3 10.3 10.4 10.5 10.6 10.6 10.7 10.7
61
Implementation ion of of finite elem lement methods ods . . . . . . . . . . . . Solving sparse linear systems . . . . . . . . . . . . . . . . . . . An outl outlin inee of of the the con converg vergen ence ce theo theory ry for for fini finite te elem elemen entt met method hodss Finit Finitee ele elem ment ent me method thodss for for eige eigen nvalue alue prob proble lem ms . . . . . . . . .
111 113 11 1166 116 116 117 117 117 119 121
. . . .
. . . .
. . . .
. . . .
. . . .
. . . .
. . . .
. . . .
. . . .
. . . .
. . . .
. . . .
. . . .
. . . .
. . . .
. . . .
. . . .
. . . .
121 122 124 124 124 124
ii
7.5
CONTENTS
Finit initee diffe ifference methods ods for the wave equation . . . . . . . . . . . . . . . . . . . . . . . . . . .
8 First-order PDEs and the Method of Characteristics
8.1 8.1 8.2 8.3
63
The The sim simples plestt PDE PDE and and the meth ethod of chara aracte cteris ristics tics . . . . . . . . . . . . . . . . . . . . . . . . First-order quasilinear PDEs . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . Burgers’s equation . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . .
9 Green’s Functions
9.1 9.2 9.2 9.3 9.3 9.4 9.5 9.6
Green’s functions for BVPs in ODEs: Spe Special cases . . . . Gre Green’s en’s func unction tionss for BV BVPs Ps in ODEs ODEs:: the the sym symmetri etricc case ase Gre Green’s en’s func unction tionss for BV BVPs Ps in ODEs ODEs:: the the gene genera rall case case . . Introdu oduction ion to Green’s functions for IVPs . . . . . . . . . Green’s functions for the heat equation ion . . . . . . . . . . . Green’s functions for the wave equation . . . . . . . . . .
. . . . . .
. . . . . .
. . . . . .
. . . . . .
. . . . . .
. . . . . .
. . . . . .
. . . . . .
. . . . . .
. . . . . .
. . . . . .
. . . . . .
. . . . . .
. . . . . .
. . . . . .
. . . . . .
. . . . . .
. . . . . .
. . . . . .
. . . . . .
. . . . . .
Proper perties ies of the Sturm-Liou iouville lle oper perator . . . . . . . . Num Numeric ericaal Me Meth thod odss for Stur Sturm m-Lio Liouvil uville le pro proble blems . . . . Examples of Sturm-Liou iouville lle problems . . . . . . . . . . R ob obin boundary conditions . . . . . . . . . . . . . . . . Finit Finitee ele elem ment ent me method thodss for for Robi Robin n bou boun ndary dary cond condit itio ions ns . The The the theory ory of Stur Sturm m-Lio -Liouv uvil ille le prob proble lem ms: an an outl outlin inee . . .
. . . . . .
. . . . . .
. . . . . .
. . . . . .
. . . . . .
. . . . . .
. . . . . .
. . . . . .
. . . . . .
. . . . . .
. . . . . .
. . . . . .
. . . . . .
. . . . . .
. . . . . .
. . . . . .
. . . . . .
. . . . . .
. . . . . .
. . . . . .
. . . . . .
. . . . . .
Phys Physic icaal mod model elss in two or thr threee sp spatia atiall dim dimen ensi sioons . . . . . . Fourier ier series on a rectangular domain . . . . . . . . . . . . . F Foourier series on a disk . . . . . . . . . . . . . . . . . . . . . . Finite elements in two dimensions . . . . . . . . . . . . . . . . The The fr freeee-spac spacee G Gre reen en’s ’s fun functio ction n fo for the the La Lapl plaacian cian . . . . . . . The The Gre Green en’s ’s func functi tion on for for the the La Lapl plac acia ian n on on a boun bounde ded d dom domai ain n Green’s fu function for the wave equation . . . . . . . . . . . . . Green’s fu functions for the heat equation ion . . . . . . . . . . . . .
. . . . . . . .
. . . . . . . .
. . . . . . . .
. . . . . . . .
. . . . . . . .
. . . . . . . .
. . . . . . . .
. . . . . . . .
. . . . . . . .
. . . . . . . .
. . . . . . . .
. . . . . . . .
. . . . . . . .
. . . . . . . .
. . . . . . . .
. . . . . . . .
. . . . . . . .
. . . . . . . .
. . . . . . . .
The The complex Fourier series . . . . . . . . . . . F Foourier series and the FFT . . . . . . . . . . . Rela Relati tion onsh ship ip of sine sine and and cos cosin inee ser serie iess to to the the full full Poin Pointtwise wise con convergen rgence ce of Fouri ourier er ser series ies . . . . . Uniform convergence of Fourier series ies . . . . . Mean Mean--squa square re con converge ergenc ncee of Fourie urierr seri series es . . . A no note abou bout general ei eigenvalue problems . . .
89 90 94 96 100 100 10 1022 107 109 111
. . . . . . . . . . . . . . . . Fou Fouri rier er seri series es . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . .
. . . . . . .
. . . . . . .
. . . . . . .
. . . . . . .
. . . . . . .
. . . . . . .
. . . . . . .
. . . . . . .
. . . . . . .
. . . . . . .
. . . . . . .
. . . . . . .
. . . . . . .
. . . . . . .
. . . . . . .
. . . . . . .
. . . . . . .
. . . . . . .
. . . . . . .
13 More ab out Finite Element Methods
13.1 13.2 13.3 13.3 13.4 13.4
83 84 85 85 86 88 89
12 More ab out Fourier Series
12.1 12.2 12.3 12.3 12.4 12.4 12.5 12.6 12.6 12.7
71 72 75 76 79 80 83
11 Problems in Multiple Spatial Dimensions
11.1 11.1 11.2 11.3 11.4 11.5 11.5 11.6 11.6 11.7 11.8
63 65 67 71
10 Sturm-Liouville Eigenvalue Problems
10.2 10.3 10.3 10.4 10.5 10.6 10.6 10.7 10.7
61
Implementation ion of of finite elem lement methods ods . . . . . . . . . . . . Solving sparse linear systems . . . . . . . . . . . . . . . . . . . An outl outlin inee of of the the con converg vergen ence ce theo theory ry for for fini finite te elem elemen entt met method hodss Finit Finitee ele elem ment ent me method thodss for for eige eigen nvalue alue prob proble lem ms . . . . . . . . .
111 113 11 1166 116 116 117 117 117 119 121
. . . .
. . . .
. . . .
. . . .
. . . .
. . . .
. . . .
. . . .
. . . .
. . . .
. . . .
. . . .
. . . .
. . . .
. . . .
. . . .
. . . .
. . . .
121 122 124 124 124 124
Chapter 1
Classification of Differential Equations 1.
(a) First-order First-order,, homogeneous, homogeneous, linear linear,, nonconstan nonconstant-coeffici t-coefficient, ent, scalar scalar ODE. ODE. (b) Second-order, inhomogeneous, linear, constant-coefficient, scalar PDE. (c) First-order First-order,, nonlinear, scalar scalar PDE.
3.
(a) Second-orde Second-order, r, nonlinear, nonlinear, scalar scalar ODE. ODE. (b) First-order, inhomogeneous, linear, constant-coefficient, constant-coefficient, scalar PDE. (c) Second-order, inhomogeneous, linear, nonconstant-coefficient, scalar PDE.
5.
(a) (a) No No. (b) Yes. (c) No. No.
7. There is only one one such f such f :: f (t) = t cos(t cos(t). 9. For any constant constant C = 1, the function w function w((t) = C u(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/(cmK). 3. The units of Aρcu∆x are g J cm2 3 Kcm = J. cm g K r 5. ∂u (, t) = Aκ , t > t 0 . ∂x 7. Measuring x in centimeters, the steady-state temperature distribution is u(x) = 0.1x + 20, the solution of 2
− ddxu = 0, u(0) = 20, u(100) = 30. 2
The heat flux is
−κ du (x) = −0.802 · 0.1 = −0.0802 W/cm . dx 2
. Since the area of the bar’s cross-section is π cm 2 , 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 (, t) = α (u(, t) − u ) , −κ ∂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)
(0, t) = αu . − κ ∂u ∂x 0
(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 ∂u ∂ 2 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
−
3
4
CHAPTER 2. MODELS IN ONE DIMENSION 13.
(a) The rate is 2
−πr D
u
−u
0
. (b) The rate of mass transfer varies inversely with and directly with the square of r.
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
−
() + pA = 0, −Ak() du dx or
du () = p. dx The other boundary condition, of course, is u(0) = 0. 3. (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 195 pA 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 k()
2 9d u dx2
−195 · 10
=
0, 0 < x < 1,
u(0) = 0, du 195 109 (1) = 109 . 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 ∂ 2 u 11 ∂ u 7.9 103 1.95 10 = f (x, t), 0 < x < 1. ∂t 2 ∂x 2 The units of f are N/m3 (force per unit volume). The units of the first term on the left are
·
·
−
·
kg m kgm/s2 N = = 3, m3 s2 m3 m and the units of the second term on the left are N 1 N = 3. 2 m m m 7. (a) We have ∂ 2 u (x, t) = c2 θ2 u(x, t), ∂t 2 while ∂ 2 u (x, t) = θ2 u(x, t). ∂x 2 Therefore, 2 ∂ 2 u 2 ∂ u c = c2 θ2 u + c2 θ2 u = 0. ∂t 2 ∂x 2 (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, . . . .
± ±
5
2.3. THE WAVE EQUATION FOR A VIBRATING STRING
2.3
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 ∂u T (, 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 (x ct), ∂x ∂v (x, t)= cf (x ct), ∂t
−
−
−
∂ 2 v (x, t)=f (x ct), ∂x 2 ∂ 2 v (x, t)=c2 f (x ct), 2 ∂t
−
−
Therefore,
2 ∂ 2 u 2 ∂ u (x, (x, t) = c2 f (x t) c ∂t 2 ∂x 2 as desired. Similarly, if v(x, t) = f (x + ct), then
−
∂v (x, t)=f (x + ct), ∂x ∂v (x, t)=cf (x + ct), ∂t and hence
2.4
∂ 2 v (x, t) ∂t 2
2
−c
2
− ct) − c f (x − ct) = 0,
∂ 2 v (x, t)=f (x + ct), ∂x 2 ∂ 2 v (x, t)=c2 f (x + ct), 2 ∂t
∂ 2 v (x, t) = c2 f (x + ct) 2 ∂x
2
− c f (x + ct) = 0.
Advection; kinematic waves
1. The solution u of ∂u ∂u +c = 0, 0 < x < , t > 0, ∂t ∂x u(x, 0) = u 0 (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) = u 0 (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 t 0 , where c(t t0 ) = x, that is, t 0 = t (1/c)x. Thus u(x, t) = φ(t (1/c)x) if x ct < 0. It is easy to verify by a direct calculuation that
−
−
u(x, t) = solves the given IBVP.
−
− −
− − u0 (x φ t
ct), ,
1 x c
−
x x
−
− ct > 0, − ct < 0
3. Suppose u solves ∂u ∂u +c = 0, < x < , t > 0, ∂t ∂x u(x, 0) = φ(x),
−∞
∞
−∞
where c > 0, and suppose that φ(x) = 0 for all x x ct a, that is, if t > (x a)/c.
− ≤
−
≤ a.
∞
The solution is u(x, t) = φ(x
− ct), so u(x, t) = 0 if
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 where a = f (1). 3. Let f : R 2
2
→R
be defined by f (x) =
x21 + x22 x2 x21
−
→ R is linear, then f (x) = ax,
If we take x = (1, 1), then f (x) = (2, 0), while 2x = (2, 2) and hence f (2x) = (8, 2) = 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
7. Suppose u
P , but not a subspace .)
1
∈ C [a, b] is any nonzero function. Then L(2u) = 2
n
du du + (2u)3 = 2 + 8u3 , dx dx
while 2Lu = 2
du + 2u3 . dx
Since L(2u) = 2Lu, L is not linear.
9. Define K : C 2 [a, b]
→ C [a, b] by
Ku = x 2
d2 u dx2
+ 3u. − 2x du dx
Then K is linear: d2 d (αu) 2x (αu) + 3(αu) 2 dx dx 2 d u du = αx2 2 2αx + 3αu dx dx 2 d u du = α x2 2 2x + 3u = αKu, dx dx
K (αu) = x 2
−
−
−
d2 d (u + v) 2x (u + v) + 3(u + v) dx2 dx d2 u d2 v du dv = x 2 + 2x + + 3(u + v) 2 2 dx dx dx dx
K (u + v) = x 2
=
−
d2 u x2 2 dx
−
− −
du d2 v 2x + 3u + x2 2 dx dx
7
dv 2x + 3v = K u + Kv . dx
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 )
2
→ C(R ) by
Lu = ρc
2
∂u ∂t
− κ ∂ ∂xu . 2
We prove that L is linear as follows: ∂ L(αu) = ρc (αu) ∂t ∂ L(u + v) = ρc (u + v) ∂t
−
∂ 2 ∂u κ 2 (αu) = αρc ∂x ∂t
−
∂ 2 κ 2 (u + v) = ρc ∂x
−
= 13.
∂ 2 u ∂u ακ 2 = α ρc ∂x ∂t
∂u ∂ v + ∂t ∂t
∂u ρc ∂t
−
−
∂ 2 u ∂ 2 v + ∂x 2 ∂x 2
− κ
∂ 2 u κ 2 ∂x
∂ 2 u κ 2 = αLu, ∂x
∂v + ρc ∂t
−
∂ 2 v κ 2 = Lu + Lv. ∂x
(a) The proof is a direct calculation: A(αx + β z)
= = = = = =
a11 a21
a12 a22
a11 a21
a12 a22
α
z 1 z 2
+ β
αx1 + βz 1 αx2 + βz 2
a11 (αx1 + βz 1 ) + a12 (αx2 + βz 2 ) a21 (αx1 + βz 1 ) + a22 (αx2 + βz 2 )
α(a11 x1 + a12 x2 ) + β (a11 z 1 + a12 z 2 ) α(a21 x1 + a22 x2 ) + β (a21 z 1 + a22 z 2 )
a11 x1 + a12 x2 a21 x1 + a22 x2 αAx + β Az. α
x1 x2
a11 z 1 + a12 z 2 a21 z 1 + a22 z 2
+ β
This shows that the operator defined by A is linear. (b) We have n
(A(αx + β y))i =
aij (αxj + βy j )
j=1
and (αAx + β Ay)i = α
n
aij xj + β
j=1
Now,
aij yj .
j=1
n
(A(αx + β y))i
=
aij (αxj + βy j )
j=1 n
=
(αaij xj + βa ij yj )
j=1 n
=
αaij xj +
j=1
n
=
α
j=1
=
n
βa ij yj
j=1
n
aij xj + β
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)
which is a line in the plane. See Figure 3.1.
∈ R} ,
: α
3
2
1
2
x
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 =
Ax
=
1 1
1 1
x1 + 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 = b ). 7. The null space of L is the set of all first-degree polynomials:
N (L) = {u : [a, b] → R
9.
: u(x) = mx + c for some m, c
∈ R} .
(a) The only functions in C 2 [a, b] that satisfy
2
− ddxu = 0 2
are functions of the form u(x) = C 1 x + C 2 . The Neumann boundary condition at the left endpoint implies that C 1 = 0, and the Dirichlet boundary condition at the right endpoint implies that C 2 = 0. Therefore, only the zero function is a solution to L m ˜ u = 0, and so the null space of L m ˜ is trivial. 2 (b) Suppose f C [a, b] is given and u C m ˜ [a, b] satisfies
∈
∈
2
− ddxu (x) = f (x), a < x < b. 2
Integrating once yields, by the fundamental theorem of calculus, x
x d2 u (s) ds = f (s) ds 2 a dx a x du du (x) + (a) = f (s) ds dx dx a 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: du (x) = dx b
⇒ ⇒
u(b)
− u(x) = a
⇒
f (s) ds
a
b
du (z ) dz = dx
x
x
−
x b
x
f (s) dsdz
a z
f (s) dsdz
a
z
u(x) =
z
− −
x
f (s) dsdz.
a
This shows that L m ˜ u = f has a unique solution for each f
∈ C [a, b].
11. By the fundamental theorem of calculus,
x
u(x) =
f (s) ds
a
satisfies Du = f . However, this solution is not unique; for any constant C , x
u(x) =
f (s) ds + C
a
is another solution.
3.3 1.
Basis and dimension (a) Both equal (14, 4, 5). (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: 2
c1 (1
2
2
− x + 2x ) + c (1 − 2x ) + c (1 − 3x + 7x ) = 0 ⇔(c + c + c ) + (−c − 3c )x + (2c − 2c + 7c )x = 0 1
⇔
2
3
c1 + c2 + c3 c1 3c3 2c1 2c2 + 7c3
− − −
2
1
3
3
1
2
2
3
= 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.
P
7. We know that 2 has dimension 3, and therefore it suffices to show either that the given set of three vectors spans 2 or that it is linearly independent. If p 2 , then we write
P
∈ P
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 ) = c 1 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) = c 1 L1 (x) + c2 L2 (x) + c3 L3 (x), and we have shown that every p proof.
∈ P can be written as a linear combination of L , L , L . 2
1
2
3
This completes the
11
3.4. ORTHOGONAL BASES AND PROJECTIONS
2 9. Let L : C N [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 C N [a, b] belongs to the null space of L if and only if it satisfies
→
∈
d2 u (x) = 0 for all x dx2
du (0) = () = 0. ∈ [a, b], du 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 u 1 (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
=
√ 13 √ 13 + √ 13 √ 13 + √ 13 √ 13
=
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 + 0 v2 + v3 . 3 6
·
·− √
√
·
3. We have 2
x + y
=
(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)
=
x
2
+ 2(x, y) + y 2 .
2
= x
Therefore,
x + y
if and only if ( x, y) = 0.
2
2
+ y
5. First of all, if (y, z) = 0 for all z then since w i
∈ W for each i, we have, in particular,
∈ W,
(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 n
z =
i=1
αi wi .
12
CHAPTER 3. ESSENTIAL LINEAR ALGEBRA We then have n
(y, z)
·
=
y,
αi wi
i=1
n
=
αi (y, wi )
i=1 n
=
α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
P is
or
2
√
2 2 5(π 2 q 1 (x) + 0q 2 (x) + π π3 2 10(π2 12) + π π3
−
− 12) q (x) 3
2
2
− 60 π π− 12 x + 60 π π− 12 x . 3
2
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 u1 =
− 0.8 0.6
, u2 =
0.6 0.8
,
13
3.5. EIGENVALUES AND EIGENVECTORS OF A SYMMETRIC MATRIX respectively. The solution is x =
·
b u1
λ1
u1 +
·
b u2
λ2
u2 =
3. The eigenvalues and eigenvectors of A are λ 1 = 2, λ2 = 1, λ3 =
u1 =
The solution of the Ax = b is
√ 13 √ 13 √ 13
·
b u1
x =
, u2 =
u1 +
λ1
·
−
b u2
λ2
√ 12
0
√ 12
u2 +
·
1 1
.
−1 and
, u3 =
b u3
λ3
u3 =
√ 16 √ 26 √ 16
−
.
− 3/2 1/2 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
∗
7.
(a) For k = 2, 3, . . . , n
2
− 1, we have
Ls(j)
k
1 ( sin((k 1) jπh) + 2sin (kjπh) sin((k + 1) jπh)) h2 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 2cos( jπh))sin(kjπh) h2 1 (2 2cos( jπh)) s(j) k . 2 h
=
− −
=
−
−
−
−
− −
= = Thus
2
− 1) = 4n − n = O(4n ).
1 (2 2cos( jπh)) s(j) 1. k , k = 2, 3, . . . , n h2 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 2cos( jπh) λj = . h2
(j)
Ls
k
=
−
−
−
(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 =
· · (1)
s
(1)
b s
+ s(2) b s(2) +
while the solution x of Lx = b is
x =
s(1) x
·
λ1
s(1) +
s(2) x
·
λ2
s(2) +
··· +
(n)
s
(n)
·b
· ·· + s λ · x s n
(n)
(n)
s
,
.
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
du . dt
x1 = u, x2 = Then dx1 dt dx2 dt 3. Define x1 = u, x2 =
=
x2 ,
=
− 12 x . 1
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
5. Define
− x + sin t. 1
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 ,
=
f 1 (x1 , x3 , x2 , x4 ),
=
x4 ,
=
f 2 (x1 , x3 , x2 , x4 ).
7. Let u 1 (t) = e −t , u 2 (t) = e −2t .
15
16
CHAPTER 4. ESSENTIAL ORDINARY DIFFERENTIAL EQUATIONS (a) We have du1 (t) = dt
t
−e− ,
d2 u1 (t) = e−t , dt2
and therefore d2 u1 du1 (t) + 3 (t) + 2u(t) = e−t dt2 dt
t
− 3e−
+ 2e−t = 0.
Similarly, du2 (t) = dt
2t
−2e−
,
d2 u2 (t) = 4e−2t , dt2
and thus d2 u2 du2 (t) + 3 (t) + 2u(t) = 4e−2t 2 dt dt
− 6e−
2t
+ 2e−2t = 0.
(b) The Wronskian of u 1 , u 2 is W (t) =
u1 (t) du1 (t) dt
u2 (t) du2 (t) dt
{
}
=
e−t e−t
e−2t 2e−2t
−
−
Since W (t) = 0 for all t, we see that u1 , u2 is linearly independent. 9.
=
−e−
3t
.
(a) Let u 1 , u2 be two functions defined on an interval, let t 0 be a point in that interval, and suppose
u1 (t0 ) du1 (t0 ) dt
u2 (t0 ) du2 (t0 ) dt
= 0.
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 c 1 du (t) + c2 du (t) = 0 for all t. But then, taking dt dt t = t 0 , we obtain
{
}
{
∈
}
c1 u1 (t0 ) + c2 u2 (t0 ) = 0, du1 du2 c1 (t0 ) + c2 (t0 ) = 0. 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 ) = (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 u 1 (t) = cos(t), u 2 (t) = cos(2t). Then
u1 (t) du1 (t) dt
u2 (t) du2 (t) dt
−
=
−
cos(t) sin(t)
cos (2t) 2sin(2t)
−
= cos(2t)sin(t)
{
}
− 2cos(t)sin(2t).
We see that W (π/2) = 1 = 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) = 0 for all t, unless u 1 and u 2 are both solutions to the same second-order linear differential equation.
{
}
17
4.2. SOLUTIONS TO SOME SIMPLE ODES
4.2 1.
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 w = αu + βv satisfies
∈ S and α, β ∈ R, then
a
d2 w dw +b + cw 2 dt dt
= = = =
d 2 d [αu + βv] + b [αu + βv] + c(αu + βv) 2 dt dt 2 2 d u d v du dv a α 2 + β 2 + b α + β + c(αu + βv) dt dt dt dt a
d2 u du d2 v dv + cu + β a 2 + b + cv α a 2 +b dt dt dt dt α 0 + β 0 = 0.
·
·
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 = 0), the solution space is spanned by two functions. Thus S is finite-dimensional and it has dimension at most 2. Moreover, it is easy to see that each of the sets er1 t , er2 t (r1 = 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 c 1 , c2 . We wish to show that, given k 1 , k2 R, there is a unique choice of c1 , c2 such that u(0) = k1 , du/dt(0) = k 2 . We have
{
}
−
−
∈
du (t) = rc 1 ert + c2 ert + rc2 tert . dt Therefore, the equations u(0) = k 1 , du/dt(0) = k 2 simplify to
1 r
0 1
c = k.
Since the coefficient matrix is obviously nonsingular (regardless of the value of r), there is a unique solution c 1 , c2 for each k 1 , 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, d dt
t
t0
t
d at e − f (s) ds = e a(t s)
dt
e−as f (s) ds
t0
=
aeat
=
t
e−as f (s) ds + eat e−at f (t)
t0
t
a
ea(t−s) f (s) ds + f (t).
t0
Therefore, with u(t) = u 0 ea(t−t0 ) +
t
ea(t−s) f (s) ds,
t0
we have du (t) dt
=
au0 ea(t−t0 ) + a
t
ea(t−s) f (s) ds + f (t)
t0
=
au(t) + f (t).
Thus u satisfies the differential equation. We also have u(t0 ) = u 0 ea(t0 −t0 ) +
t0
t0
so the initial condition holds as well.
ea(t−s) f (s) ds = u 0 + 0 = u 0 ,
18
CHAPTER 4. ESSENTIAL ORDINARY DIFFERENTIAL EQUATIONS 9. u(t) = 12 sin 2 (t)
11. u(t) = 13.
1 2
e−t + cos (t) + sin (t)
2
(a) The solution is x(t) = e t
/2
.
(b) The solution is x(t) = (1 + cos (1)
3
− cos(t))/t . 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
d2 u du +b + cu = f (t), 2 dt dt
where u 1 , u 2 are solutions of the homogeneous version of this ODE, and assume that dc1 dc2 u1 + u2 = 0. dt dt We then have du du1 du2 = c 1 + c2 , dt dt dt d2 u d2 u1 d2 u2 dc 1 du1 dc 2 du2 = c + c + + . 1 2 dt2 dt2 dt2 dt dt dt dt Substituting into the ODE, we obtain d2 u1 d2 u2 dc 1 du1 dc 2 du2 a c1 2 + c2 2 + + dt dt dt dt dt dt
⇒ ⇒ c1
d2 u1 du1 + cu1 a 2 +b dt dt
du1 du2 + b c1 + c2 dt dt
d2 u2 du2 + cu2 a 2 +b dt dt
+ c2
dc1 du1 dc 2 du2 + = f (t) dt dt dt dt dc1 du1 dc2 du2 + = a −1 f (t), dt dt dt dt
a
⇒
+a
+ c(c1 u1 + c2 u2 ) = f (t)
dc1 du1 dc 2 du2 + = f (t) dt dt dt dt
as desired. Notice that we have used the fact that u 1 , u 2 solve the homogenous ODE: a
d2 u1 du1 d2 u2 du2 + + = 0, +b + cu2 = 0. b cu a 1 2 2 dt dt dt dt
17. Let us define v(s) = u(es ); then u(t) = v(ln(t)). It follows that du 1 dv d2 u 1 d2 v (t) = (ln(t)), (t) = (ln (t)) dt t ds dt2 t2 ds2
(ln(t)). − t1 dv ds 2
It follows that d2 u du (t) + at (t) + bu(t) = 0 dt2 dt 2 d v dv dv (ln (t)) (ln(t)) + a (ln(t)) + bv(ln(t)) = 0 2 ds ds ds d2 v dv (ln (t)) + (a 1) (ln(t)) + bv(ln(t)) = 0 ds2 ds d2 v dv (s) + (a 1) (s) + bv(s) = 0. ds2 ds
t2
⇔ ⇔ ⇔
−
−
−
Thus the change of variables t = e s transforms the Euler equation into a constant coefficient ODE.
19
4.3. LINEAR SYSTEMS WITH CONSTANT COEFFICIENTS
4.3
Linear systems with constant coefficients
1. The solution is x(t)
√ 13 u − √ 12 e− u + √ 16 e− − e− + e− − e− .
=
t
1
=
3. The general solution is
t 1 2 1 1 3 3 t 1 1 + + e 3 2 1 3
2t
1 6 2t
1 e 2t 6
−
x(t) =
2t
2
−
c1 et + c2 e−t c1 et c2 e−t
u3
.
−
5. x0 must be a multiple of the vector (1, 1). 7. 1 x(t) = 60 9.
(a) The solution is
− 55 − 3e− − 45e− + 20t − 7cos(t) − 11sin(t) 10 + 6e− + 20t − 16cos(t) − 8sin(t) −5 − 3e− + 45e− + 20t − 37cos(t) + 19 sin (t) 2t
t
2t
2t
t
3 −t 1 3t 3 1 3t e + e , y(t) = e−t e . 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 (r + s)e−t + (r 2
− s)e
3t
, y(t) =
1 (r + s)e−t + (s 2
3t
− r)e
.
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 + θ2 v = 0, v(0; s) = 0, (0; s) = 0 dt2 dt for all s, and define
t
u(t) =
v(t
0
Then du (t) = v(0; t) + dt
t
0
d2 u dv (t) = (0; t) + 2 dt dt
dv (0; t) dt
(using the fact that v(0; t) = 0 and
dv (t dt
t
0
t
− s; s) ds =
d2 v (t dt2
= f (t) +
t
d2 v (t dt2
0
0
= f (t)
Also,
0
0
dv (t dt
− s; s) ds, t
− s; s) ds = f (t) +
t
d2 v (t dt2
0
d2 v (t dt2
− s; s) ds
= 0 for all t). We then have
d2 u (t) + θ2 u(t) = f (t) + dt2
since
− s; s) ds.
−
d2 v (t dt2
s; s) ds + θ2
t
0
2
v(t
− s; s) ds
− s; s) + θ v(t − s; s)
ds
2
− s; s) + θ v(t − s; s) = 0 for all s.
0 du dv (t) = (t s; s) ds = 0, u(t) = v(t s; s) ds = 0, dt 0 0 dt and thus u satisfies both the differential equations and the initial conditions.
−
−
20
CHAPTER 4. ESSENTIAL ORDINARY DIFFERENTIAL EQUATIONS
13. Suppose that, for all s, v(t; s) satisfies the ODE dk u dk−1 u + a (t) + − k 1 dtk dtk−1
+ a (t)u = 0 ··· + a (t) du dt 1
0
and also the initial conditions du dk−2 u dk−1 u (s) = 0, . . . , (s) = 0, (s) = f (s). dt dtk−2 dtk−1
u(s) = 0, Define
t
u(t) =
v(t; s) ds.
0
Then
du (t) = v(t; t) + dt (since v(t; t) = 0 for all t). Similarly, d2 u dv (t) = (t; t) + 2 dt dt .. .
t
0
t
0
d2 v (t; s) ds = dt2
t
dk−1 u dk−2 v (t) = (t; t) + dtk−1 dtk−2 and finally dk u dk−1 v (t) = (t; t) + dtk dtk−1 We then see that
dv (t; s) ds = dt
0
t
0
t
0
t
0
dv (t; s) ds dt
d2 v (t; s) ds, dt2
dk−1 v (t; s) ds = dtk−1
t
dk−1 v (t; s) ds dtk−1
0
dk v (t; s) ds = f (t) + dtk
t
dk v (t; s) ds. dtk
0
dk u dk−1 u du + a (t) + + a1 (t) + a0 (t)u k−1 dtk dtk−1 dt t dk v dk−1 v dv =f (t) + (t; + (t) (t; s) + + a1 (t) (t; s) + a0 (t)v(t; s) ds s) a k 1 − k k 1 − dt dt dt 0
···
=f (t), since
·· ·
dk v dk−1 v (t; + (t) (t; s) + s) a k 1 − dtk dtk−1
by assumption. Also,
dj u (0) = dtj
15. The solution to
0
0
(t; s) + a (t)v(t; s) = 0 ··· + a (t) dv dt 1
0
dj v (t; s) ds = 0, j = 0, 1, . . . , k dtj
− 1.
d2 u du +3 + 2u = 0, 2 dt 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), 2 dt dt u(0) = 0, du (0) = 0 dt
is
t
u(t) =
0
t
v(t
− s; s) ds =
0
e−(t−s)
− e−
2(t s)
−
f (s) ds.
21
4.4. NUMERICAL METHODS FOR INITIAL VALUE PROBLEMS
4.4 1.
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 u 1 = x, u 2 = dx/dt, u 3 = y, u 4 = dy/dt. Then the system is du1 dt du2 dt du3 dt du4 dt
=
u2 ,
=
u1 + 2u4
+ µ ) µ (u − µ ) − µr (u − r (u , u ) , (u , u ) 2
1
1
1
3
1 3
1
1
2
1
3
2 3
=
u4 ,
=
−2u + u − r (uµ ,uu ) − r (uµ ,uu ) 2
2 3
3
1
1
1 3
3
3
2
1
3
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 0
y
−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. 2 ) t ( x
0
−2 0 1 ) t ( y
1
2
3
4
5
6
7
4
5
6
7
3 4 Step index
5
6
7
t
0
−1 0
1
2
3 t
0.15 p e 0.1 t s e m 0.05 i T
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.
− t ) ≤ t ≤ b + (t − t ) if and only if a ≤ t − (t − t ) ≤ b, so the function v is
(a) We first note that a + (t1 well-defined. In fact,
0
1
0
1
{v(t) : t ∈ [a + (t − t ), b + (t − t )]} {u(t − (t − t )) : t ∈ [a + (t − t ), b + (t − t )]} {u(t) : t ∈ [a, b]} 1
= = (just replace t We have
0
1
0
1
0
1
0
0
1
0
− (t − t ) by t in the last step). 1
0
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
− (t − t )) = u(t ) = u . 1
0
0
0
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) = 21 (t − 1)
2
7. The exact solution is x(t) =
2
t
− 2e− , 1 t e, 2
+ 1.
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:
23
4.5. STIFF SYSTEMS OF ODES ∆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 1 x(t) = 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 x(t) = 2
.
3e−t e−100t 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 x49 > x0 , while n 50 yields xn x0 .
≤
≥
(c)
With x i+1 = x i + ∆tAxi and y 1(i)
≤
= u 1 xi , we have
·
u1 xi+1 = u 1 xi + ∆tu1 Axi
⇒ ⇒ ⇒ ⇒ A similar calculation shows that
·
y1(i+1) y1(i+1) y1(i+1) y1(i+1)
= y 1(i) = y 1(i) = y 1(i) = (1
·
·
·
+ ∆t(Au1 ) xi
·
+ ∆tλ1 u1 xi
+ ∆tλ1 y1(i) + ∆tλ1 )y1(i) .
y2(i+1) = (1 + ∆tλ2 )y2(i) .
For stability, we need
|1 + ∆tλ | ≤ |1 + ∆tλ | ≤ 1
1,
2
1.
Since the eigenvalues of A are 1, 100, it is not hard to see that the upper bound for ∆t is ∆t as was determined by experiment. (d) For the backward Euler method, a similar calculation shows that
− −
y1(i+1)
=
y2(i+1)
=
Stability is guaranteed for any ∆t, since
1 2
1
|1 − ∆tλ |− i
for every ∆t > 0.
1 (i) 1 , 1 (i) 2 .
− ∆tλ )− y (1 − ∆tλ )− y (1
< 1
≤ 0.02, just
Chapter 5
Boundary Value Problems in Statics 5.1 1.
The analogy between BVPs and linear algebraic systems (a) Since M D is linear, it suffices to show that the null space of M D is trivial. If M D 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 (M D ) = 0 . (b) If u
1 D
N
{}
D u = f ,
then we have
∈ C [0, ] and M
du (x) = f (x), 0 < x < and u(0) = 0, dx which imply that
x
u(x) =
f (s) ds.
0
But we also must have u() = 0, which implies that
f (s) ds = 0.
0
If f 3.
(a)
(b)
∈ C [0, ] does not satisfy this condition, then it is impossible for M u = f to have a solution. 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 . We have already seen that Lu = f has a solution for any f ∈ C [0, ] (see the discussion immediately preceding D
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 2
⇒ a = 0.
Then u(x) = b, and
⇒ b = 0 ⇒ b = 0. Therefore u is the zero function, and so N (L) ∩ S = {0}. The uniqueness property then follows from u(x) dx = 0
0
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 (L) S = 0 . The uniqueness property then follows from the first part of this exercise.
⇒ N ∩ { }
5. The condition
2
du (x) dx
0
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 (Lm ˜ ) is trivial.
N
(b) For any f
∈ C [0, ], the function
u(x) =
x
2 belongs to C m ˜ [0, ] and satisfies
z
f (s) dsdz
0
2
− ddxu (x) = f (x), 0 < x < . u = f has a solution for any f ∈ C [0, ], and therefore R(L 2
This shows that L m ˜
∈
(c) Suppose u, v
m ˜)
2 C m ˜ [0, ].
= C [0, ].
Then
(Lm ˜ u, v)
= =
− − −
0
0
=
=
+
0
0
du dv (x) (x) dx dx dx
du dv (x) (x) dx (since v () = 0, du/dx(0) = 0) dx dx
−
dv u(x) (x) dx
=
du (x)v(x) dx
=
d2 u (x)v(x) dx dx2
0
u(x)
0
d2 v (x) dx dx2
d2 v u(x) 2 (x) dx (since u() = 0, dv/dx(0) = 0) dx 0 (u, Lm ˜ v).
Thus L m ˜ 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)
=
0
=
du (x)u(x) dx
=
0
>
d2 u (x)u(x) dx dx2
− − du (x) dx
+
0
0
du (x) dx
2
dx
2
dx (since u() = 0, du/dx(0) = 0)
0.
2 The last step follows because every nonzero function in C m ˜ [0, ] has a nonzero derivative.
27
5.2. INTRODUCTION TO THE SPECTRAL METHOD; EIGENFUNCTIONS 9. Define u(x) = x(1
2
− x), v(x) = x (1 − x). Then a direct calculation shows that (Mu,v) =
7 , 30
(u,Mv) =
4 . 15
but
11.
(a) Suppose u, v
2 R
∈ C [0, ]. Then
(LR u, v)
− − −
=
κ
0
=
κ
d2 u (x)v(x) dx dx2
du (x)v(x) dx
+κ
0
0
du dv (x) (x) dx dx dx
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
= =
−
− −
Applying the boundary conditions on u, v, we obtain du dv dv ()v() + κ (0)v(0) + u()κ () − u(0)κ (0) −κ du 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 L R 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
(α + κ)a + αb
=
0,
=
0.
The determinant is computed as follows:
−κ
α + κ
The only solution is a = b = 0, that is, u = 0, so
5.2
α α
=
−α(2κ + α) < 0.
N (L
R)
is trivial.
Introduction to the spectral method; eigenfunctions
1. If n = m, then
(un , um )
=
− − − − nπx mπx sin dx
sin
0
= =
1 2
cos
(n
0
m)πx
1 sin 2 (n m)π
(n
−
(n + m)πx
cos
m)πx
dx
sin (n + m)π
(n + m)πx
.
0
This last expression simplifies to four terms, each including the sine function evaluated at an integer multiple of π, and hence each equal to zero. Thus (un , um ) = 0 for n = m.
3. The eigenpairs are (2n
2 2
− 1) π 42
, cos
(2n
− 1)πx 2
, n = 1, 2, 3, . . . .
28
CHAPTER 5. BOUNDARY VALUE PROBLEMS IN STATICS 5. The characteristic polynomial is r 2
− r + (λ − 5), so the characteristic roots are √ √ 1 − 21 − 4λ 1 + 21 − 4λ r = , r = . 1
2
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) = c 1 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) = c 1 er1 x + c2 er2 x . The boundary conditions lead to the system
c1 r1 e
r1
c1 + c2
=
0,
r2
=
0.
+ c2 r2 e
This system has the unique solution c 1 = 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 r1 = 2
−
1 θi, r2 = + θi, θ = 2
√ 4λ − 21 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 c 1 = 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 v 1 , v2 , as given by the above formula. A direct calculation shows that . (v1 , v2 ) = 0.2092, so v 1 and v 2 are not orthogonal.
−
29
5.3. SOLVING THE BVP USING FOURIER SERIES 7.
(a)
∞
2( 1)n+1 nπ
sin (nπx)
(b)
∞
4 sin ( nπ 2 ) n2 π2
sin(nπx)
(c)
∞
12( 1)n+1 n3 π3
−
n=1 n=1
−
n=1
sin(nπx)
n+1
∞ 720(−1) sin(nπx) (d) n=1 n5 π5 The errors in approximating the original functions using 10 terms 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
−0.01 0
1
−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 terms 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 − x 3 ; 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
∞
an cos
n=1
where (a) an =
−
(2n
− 1)πx 2
,
4(2+( 1)n (2n 1)π) ; π2 (2n 1)2
−
−
−
2
((2n−1)π/8) (b) an = 32cos ((2n−1)π/4)sin ; π2 (2n−1)2
(c) an =
8(24+12( 1)n (2n 1)π+(2n 1)2 π2 ) ; π4 (2n 1)4
(d) an
8(5760+2880( 1)n (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 terms of the Fourier quarter-wave cosine series are graphed in Figure 5.2.
5.3 1.
Solving the BVP using Fourier series ∞
(a)
n=1
(b) x + 3.
(a)
−−
∞
n=1
sin(nπx)
2(1 ( 1)n ) n3 π3
−−
32( 1)n+1 n=1 (2n 1)4 π4
∞
(b) 1 + (c)
2(1 ( 1)n ) n3 π3
− −
sin
sin(nπx)
(2n 1)πx 2
−
16(( 1)n π+2( 1)n+1 nπ 2) n=1 (2n 1)4 π4
−
∞
−
−
2(1 ( 1)n ) n=1 nπ(2+n2 π2 ) 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
−0.04 0
1
0.01
0.05
0
0
−0.01
−0.05
−0.02
−0.1
−0.03 0
0.5 x
−0.15 0
1
0.5 x
1
0.5 x
1
Figure 5.2: Errors in approximating four functions by 10 terms of their Fourier quarter-wave cosine series (see Exercise 5.2.11). 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 .
(d) x + 5. We have
2( 1)n n=1 nπ(1+n2 π2 ) sin(nπx)
∞
−
d2 u (x) = dx2
−x, 0 < x < 1,
so du (x) dx
= = =
x
du (0) + dx du (0) dx du (0) dx
− 0
d2 u (s) dx dx2
x
s dx
0 2
− x2 .
Write C for du/dx(0); then another integration yields x
=
u(x)
du (s) ds dx
−
u(0) +
0
x
=
0+
s 2 2
C
0
ds
3
=
Cx
− x6 .
Finally, we choose C so that u(1) = 0, which yields C = 1/6. This yields u(x) =
1 x 6
x3 .
−
7. The Fourier quarter-wave sine coefficients of T d2 u/dx2 are
−
2T
2 an =
bn =
−
0
d2 u (x)sin dx2
(2n
− 1)πx 2
dx,
while those of u are
0
u(x)sin
(2n
− 1)πx 2
dx.
31
5.3. SOLVING THE BVP USING FOURIER SERIES Using integration by parts twice, almost exactly as in (5.19), we can express b n in terms of a n : 2T
− =
− 2T
0
d2 u (x)sin dx2
− − − − − − − − − (2n
(2n
du (x)sin dx
1)πx 2
dx
1)πx 2
x=0
(2n 1)πx − (2n 2 1)π du (x)cos dx 2 2T (2n − 1)π du (2n 1)πx (x)cos dx
dx
0
=
=
22
2 dx (since sin (0) = du/dx() = 0)
0
−
2T (2n 1)π 22
u(x)cos
+
=
=
T (2n
2 2
− 1) π
42 T (2n
2
(2n
(2n
x=0
1)π
2
1)πx 2
0
1)πx 2
dx
1)πx dx 2 0 (since u(0) = cos ((2n 1)π/2) = 0) u(x)sin
(2n
−
2 2
1) π an . 42
−
(2n
u(x)sin
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 a 1 , a2 , a3 , . . . be the Fourier quarter-wave sine coefficients of u; then
∞
d2 u κ(2n 1)2 π 2 κ 2 (x) = an sin dx 2002 n=1
−
−
(2n
− 1)πx
200
,
where κ = 3/2. We also have
∞
0.004 0.001 = sin (2n 1)π n=1
−
(2n
− 1)πx
200
.
Setting the two series equal and solving for a n , we find
∞
3.2 102 u(x) = sin 3(2n 1)3 π3 n=1
· −
(2n
− 1)πx
200
.
The temperature distribution is graphed in Figure 5.3. 3.5
3
2.5 e r u 2 t a r e p m 1.5 e t 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 D
1. Suppose f , g
∈ C [0, ], so that
f (0) = f () = g(0) = g() = 0.
Then
− − − − − − d dx
=
k(x)
d dx
0
=
k(x)
0
=
∈
3
3
k(x)
0
0
0
0
d dx
k(x)
d dx
k(x)
dg dx
f,
− c) (d − x)
+
f (x)
0
=
df dg (x) (x) dx (integration by parts) dx dx
df dg (x) (x) dx (using g(0) = g() = 0) dx dx
5. Let v
df (x) g(x) dx dx
k(x)
dg k(x)f (x) (x) dx
=
2 [0, ] C D
,g
df k(x) (x)g(x) dx
=
3. Take p(x) = (x
df dx
d f (x) dx
dg k(x) (x) dx (integration by parts) dx
dg (x) dx (using f (0) = f () = 0) dx .
and v [c,d] as suggested in the hint.
be a test function. We have d dx d dx
− ⇒−
⇒− 0
du k(x) (x) + p(x)u(x) = f (x), 0 < x < dx du k(x) (x) v(x) + p(x)u(x)v(x) = f (x)v(x), 0 < x < dx
d dx
du k(x) (x)v(x) dx
⇒−
⇒ ⇒
k(x)
0
k(x)
0
+
0
0
p(x)u(x)v(x) dx =
0
f (x)v(x) dx
p(x)u(x)v(x) dx =
0
f (x)v(x) dx
0
p(x)u(x)v(x) dx =
0
0
du dv k(x) (x) (x) dx + dx dx
du dv (x) (x) dx + dx dx
du k(x) (x) v(x) dx + dx
f (x)v(x) dx
0
du dv (x) (x) + p(x)u(x)v(x) dx = dx dx
f (x)v(x) dx.
0
The last step follows from the fact that v (0) = v() = 0. Thus the weak form of the given BVP is find u
∈
2 [a, b] C D
such that
0
du dv k(x) (x) (x) + p(x)u(x)v(x) dx = dx dx
f (x)v(x) dx for all v
0
2 D
∈ C [a, b].
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 C D [a, b] such that
∈
0
du dv du k(x) (x) (x) + c(x) (x)v(x) + p(x)u(x)v(x) dx = dx dx dx
f (x)v(x) dx for all v
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
− 1 2
∂ ∂t
=
c
0
Aρ(x)
0
∂u ∂t
∂u ∂t
2
1 dx + 2
2
dx, t > 0.
0
Ak(x)
∂u ∂x
2
dx
2 D [a, b].
∈ C
33
5.5. THE GALERKIN METHOD
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 = 0.
(c) No, a(f, f ) = 0 for f (x) = 1, but f = 0.
3. Let
{
}
F N = span sin(πx), sin(2πx), . . . , sin(N πx) . We will apply the Galerkin method to the BVP 2
−k ddxu = f (x), 0 < x < , 2
u(0) = 0, u() = 0, using F N as the approximating subspace. The weak form of the BVP is find u
∈
2 [a, b] C D
such that
0
du dv k (x) (x) dx = dx dx
f (x)v(x) dx for all v
0
2 D [a, b],
∈ C
and the Galerkin method takes the form N
find v N =
uj sin( jπx) such that
0
j=1
dvn dv (x) (x) dx = k dx dx
f (x)v(x) dx for all v
0
∈ F
N .
We find the coefficients u j , j = 1, 2, . . . , N , by solving Ku = f , where
K ij =
0
dφj dφi (x) (x) dx, F i = k dx dx
f (x)φi (x) dx
0
and φ j (x) = sin( jπx). We have
0
dφj dφi (x) (x) dx = k dx dx
2
kijπ cos( jπx)cos(iπx) dx =
0
kj 2 π2 , 2
0,
i = j, i = 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 L 2 (0, 1) inner product. Also, notice that
f (x)φj (x) dx =
0
f (x)sin( jπx) dx =
0
cj , 2
where c j is the Fourier sine coefficient of f . Since K is diagonal, it is easy to solve Ku = f : cj /2 cj = 2 2. kj 2 π 2 /2 kj π
uj = Thus
N
vN (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
2
− 3e)x
+ (4e
− 10)x.
34
CHAPTER 5. BOUNDARY VALUE PROBLEMS IN STATICS
P
(d) The bilinear form is not an inner product on 2 , since a(1, 1) = 0 but 1 = 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) = e x in the energy norm .
−
7. Write p(x) = x(1 x) and q (x) = x(1/2 Ku = f . The stiffness matrix K is
K
= =
−
−
∈
− x)(1 − x). The approximation will be v (x) = u p(x) + u q (x), where 1
1 dp dp (1 + x) dx (x) dx (x) dx 0 1 dp dq (1 + x) dx (x) dx (x) dx 0 1 1 2 30 , 1 3 30 40
−
1 (1 + 0 1 (1 + 0
−
dq dp (x) dx (x) dx x) dx dq dq (x) dx (x) dx x) dx
2
and the load vector f is f =
Therefore,
1 0 1 0
xp(x) dx xq (x) dx
. u=
=
1 12 1 120
−
0.16412 0.038168
−
.
.
The exact solution is u(x) =
− 14 x
2
1 + x 2
− 4 l1n 2 ln (1 + x).
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 and 1003 times as long, or 106 seconds (about 11.5 days) to solve a 10000 10000 system.
×
× 1000 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 as long, or 1 second, to solve a 10000 10000 tridiagonal system.
×
× 1000 tridiagonal system, and 100 times
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 2
− ddxu = x, 2
0 < x < 1,
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 4 8 4 K = 0 4 8
−
− − − (as computed in Example 5.18 in the text). The load vector F ∈ R
3
is defined by
1
F i =
f (x)φi (x) dx.
0
A direct calculation shows that F 1 = 1/16, F 2 = 1/8, F 3 = 3/16. Solving Ku = f yields
u =
The exact solution to the BVP is u(x) = (x approximation.
3
5 128 1 16 7 128
.
− x )/6. Figure 5.5 shows the exact solution and the piecewise linear
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 error = O 5. The exact solution is u(x) = x(1
3
1 n2
.
− x )/12. Here are the errors for n = 10, 20, 40, 80:
36
CHAPTER 5. BOUNDARY VALUE PROBLEMS IN STATICS maximum error 1.1286 10−3 2.9710 10−4 7.6186 10−5 1.9288 10−5
n 10 20 40 80
· · · ·
We see that, when n is doubled, the error descreases by a factor of approximately four. Thus error = O
1 n2
.
7. The weak form is
find u
∈ V such that
du dv k(x) (x) (x) + p(x)u(x)v(x) dx dx dx
0
=
f (x)v(x) dx for all v
0
The bilinear form is
∈ V.
du dv a(u, v) = k(x) (x) (x) + p(x)u(x)v(x) dx. dx dx 0 9. We wish to apply the piecewise linear finite element method to the BVP 2
− ddxu + 2u = 21 − x, 0 < x < 1, 2
u(0) = 0, u(1) = 0. We use a sequence of increasingly fine uniform meshes. The weak form of the BVP is u
∈
2 C D [0, 1],
1
0
du dv (x) (x) dx + dx dx
1
1
2u(x)v(x) dx =
0
− 1 2
0
x v(x) dx for all v
2 D
∈ C [a, b],
and the Galerkin method requires the solution of 1
∈ V ,
vn
n
0
dvn dv (x) (x) dx + 2 dx dx
1
1
vn (x)v(x) dx =
0
− 0
1 2
x v(x) dx for all v
The second integral on the left side of the variational equation leads to the matrix M 1
M ij =
∈ V . n
(n 1) (n 1)
∈R
− × − , where
φ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
− − M ij =
2h , 3
j = i, frach6, 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 1 1 h f i = 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 2
−D ddxu 2
=
10−7 , 0 < x < 100,
u(0)
=
0,
u(100)
=
0.
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 ∞ 2 10−3 (1 ( 1)n ) nπx us (x) = sin . 3 3 Dn π 100 n=1 3. The solution is
·
u(x, t) = A graph of u( , 0.1) is given in Figure 6.1.
·
→ ∞ −−
∞
2( 1)n+1 −n2 π2 t e sin(nπx). nπ n=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) + and v (x, t) = u(x, t)
x (h(t)
− g(t))
− p(x, t), we have ρc
∂v ∂t
2
∂ v − κ ∂x
2
= =
ρc
∂u ∂t
2
− κ ∂ ∂xu −
f (x, t)
2
− ρc 39
ρc
∂p ∂t
2
∂ p − κ ∂x
dg x (t) + dt
2
dh (t) dt
−
dg (t) dt
.
40
CHAPTER 6. HEAT FLOW AND DIFFUSION We also have v(x, t0 ) = u(x, t0 )
− p(x, t ) = ψ(x) − g(t ) − x (h(t ) − g(t )) 0
0
0
0
and v(0, t)
=
v(, t)
=
− p(0, t) = g(t) − g(t) = 0, u(, t) − p(, t) = h(t) − h(t) = 0. u(0, t)
We define g(x, t) = f (x, t)
− ρc
and φ(x) = ψ(x)
dg x (t) + dt
dh (t) dt
−
dg (t) dt
− g(t ) − x (h(t ) − g(t )). 0
0
0
Then v satisfies ρc
∂v ∂t
2
∂ v − κ ∂x
2
=
g(x, t), 0 < x < , t > 0,
v(x, t0 )
=
φ(x), 0 < x < ,
v(0, t)
=
0, t > t0 ,
v(, t)
=
0, t > t0 .
−
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 ∂t
2
∂ v − ∂x
2
=
x sin(t), 0 < x < 1, t > 0,
v(x, 0)
=
0, 0 < x < 1,
v(0, t)
=
0, t > 0,
v(1, t)
=
0, t > 0.
This IBVP has solution
∞
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 + 1) nπ π 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 ∂t
2
− κ ∂ ∂xu 2
=
0, 0 < x < 100, t > 0,
u(x, 0)
=
5, 0 < x < 100,
u(0, t)
=
0, t > 0,
u(100, t)
=
0, t > 0.
41
6.2. PURE NEUMANN CONDITIONS AND THE FOURIER COSINE SERIES The solution is u(x, t) =
∞
10(1
n
− (−1) nπ
n=1
κn2 π 2 ) − 100 nπx t 2 ρc e sin . 100
The temperature at the midpoint after 20 minutes is
. u(50, 1200) = 1.58 degrees Celsius. 11.
(a) The steady-state temperature is u s (x) = x/20. (b) The temperature u(x, t) satisfies the IBVP ρc
∂u ∂t
2
− κ ∂ ∂xu 2
=
0, 0 < x < 100, t > 0,
u(x, 0)
=
5, 0 < x < 100,
u(0, t)
=
0, t > 0,
u(100, t)
=
5, t > 0.
The solution is u(x, t) =
∞
κn2 π 2 10 − 100 x nπx t 2 ρc + e 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: κπ2 πx t . x 10 − 100 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
A graph shows that
πx 100
200sin πx
for all x, so we need
e
200sin πx
κπ 2 t − 100 2 ρc
e
κπ2 − 100 2 ρc t
πx 100
≤
≤ 0.01. 2
≤ 0.005.
This yields 2
t
ln 0.005 . ≥ −100 ρc = 4520. κπ 2
About 75 minutes and 20 seconds are required.
6.2
Pure Neumann conditions and the Fourier cosine series
1. The solution is u(x, t) = d 0 +
∞
2
dn e−n
π2 t
cos(nπx),
n=1
where
1
d0 =
0
x(1
−
1 x) dx = , dn = 2 6
The graphs are given in Figure 6.3.
1
x(1
0
−
2 ( 1)n+1 x)cos(nπx) dx = n2 π2
−
−1
.
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 terms in the Fourier series. 3.
(a) If u, v
2 N [0, ],
∈ C
then
(LN u, v)
= =
−
0
dv (x) dx
u(x)
=
u(x)
0
=
+
0
0
du dv (x) (x) dx dx dx
du dv (x) (x) dx (since dx dx
0
=
− − − du (x)v(x) dx
=
d2 u (x)v(x) dx dx2
u(x)
0
0 2
du (0) dx
du () dx
=
= 0)
d2 v (x) dx dx2
d v (x) dx (since dx2
dv (0) dx
=
dv () dx
= 0)
(u, LN v).
This shows that L N is symmetric. (b) Suppose λ is an eigenvalue of L N and u is a corresponding eigenvector, normalized so that (u, u) = 1. Then λ = λ(u, u)
=
(λu,u)
=
(LN u, u)
=
−
=
0
du (x) dx
0
≥
− du (x)u(x) dx
=
d2 u (x)u(x) dx dx2 +
0
0 2
du (x) dx
dx (since
2
du (0) dx
dx =
du () dx
= 0)
0.
Thus L N cannot have any negative eigenvalues. 5.
(a) Suppose u is a solution to the BVP. Then
f (x) dx =
0
−κ
0
d2 u (x) dx = κ dx2
du (0) dx
This is the compatibility condition:
f (x) dx = κ(a
0
2 (b) The operator K : C N [0, ]
− b).
→ C [0, ] defined by 2
Ku =
−κ ddxu + u 2
du () − dx
= κ(a
− b).
43
6.2. PURE NEUMANN CONDITIONS AND THE FOURIER COSINE SERIES has eigenpairs κn2 π 2 nπx , γ n (x) = cos , n = 1, 2, 3, . . . . 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 key is that 0 is not an eigenvalue of K , as it is of L N .)
∈ C [0, ]. (The
7. As shown in this section, the solution to the IBVP is u(x, t) = d 0 +
∞
2
dn e−n
π2 t/2
cos
n=1
nπx ,
where d 0 , d1 , d2 , . . . are the Fourier cosine coefficients of ψ . We have
∞
2
dn e−n
2
nπx cos
2
π t/
n=1
∞
| |
≤
n=1
∞
≤
2 2 2 nπx dn e−n π t/ cos
|
2 2 2 dn e−n π t/
|
n=1
∞
2 2 e−π t/
≤
|
2 2 2 dn e−(n −1)π t/ .
n=1
|
This last series certainly converges (as can be proved, for example, using the comparison test), and e−π This shows that d0 +
∞
2
dn e−n
2
t/2
→ 0 as t → ∞.
π2 t/2
cos
n=1
The limit d 0 is
1 the average of the initial temperature distribution.
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
∞
nπx 2 u(x, t) = an (t)sin , an (t) = n=1
u(x, t)sin
0
nπx dx
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:
−
= = = =
−
2
−
2
2nπ 2
0
∂ 2 u nπx (x, t)sin dx 2 ∂x
−
∂u nπx (x, t)sin ∂x
0
0
2nπ (( 1)n u(, t) 2
0
+
nπ
u(x, t)sin
0 2 2
− u(0, t)) + nπ 2
∂u nπx (x, t)cos dx ∂x
∂u nπx (x, t)cos dx ∂x
2nπ nπx u(x, t)cos 2
−
0
nπ
nπx dx
an (t).
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 terms of a 1 (t), a2 (t), . . ..
−
44
CHAPTER 6. HEAT FLOW AND DIFFUSION (b) The Fourier sine series of u(x, t) = t is
∞
2 (1
n
− (−1) ) t sin(nπx), nπ
n=1
and the formal calculation of the sine series of ∂ 2 u/∂x2 yields
−
∞
2 (1
n=1
n
− (−1)
However,
) nπt sin(nπx).
2
− ∂ ∂xu (x, t) = 0, 2
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
∂u ∂t
2
− κ ∂ ∂xu 2
=
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) = c 0 +
∞
κn 2
cn e− 25ρc t cos
n=1
where
nx , 5
π 4 , 9 10( 1)n+1 , n = 1, 2, 3, . . . . n4
c0
=
25 +
cn
=
−
(c) The steady-state temperature is the constant u s = 25 + π4 /9 degrees Celsius. (d) We must choose t so that
|u − u(x, t)| ≤ 0.01 |u | s
s
holds for every x 5.
∈ [−5π, 5π]. By trial and error, we find that about 360 seconds (6 minutes) are required.
(a) To show that L p 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 du (x)v(x) = ( )v( ) dx − dx
du − − − dx ()v().
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.
45
6.3. PERIODIC BOUNDARY CONDITIONS AND THE FULL FOURIER SERIES (b) Suppose L p u = λu, where u has been normalized: (u, u) = 1. Then
λ = λ(u, u) = (λu,u) = (Lp u, u)
=
d2 u (x)u(x) dx − dx2
− −
= =
≥
du (x)u(x) + dx − − 2 du (x) dx − dx 0.
du (x) dx
2
dx
The boundary terms 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
f (x) dx =
−
∂ 2 u (x, t) dx = κ − ∂x 2
−
−
∂u (x, t) = κ κ ∂x −
∂u ( , t) ∂x
−
−
∂u (, t) = 0. ∂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
−
cf (x) dx = 0 for all c
−
∈ R,
or simply
f (x) dx = 0.
−
9. Consider the IBVP ρc
∂u ∂t
2
− κ ∂ ∂xu + pu = f (x, t), − < x < , t > t , u(x, t ) = ψ(x), − < x < , u(−, t) = u(, t), t > t , ∂u ∂u (−, t) = (, t), t > t . ∂x ∂x 0
2
0
0
0
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 u(x) = a 0 +
∞
nπx nπx + bn sin
an cos
n=1
.
Substituting this series into the left side of the PDE yields ∂u ∂ 2 u (x, t) κ 2 (x, t) + pu(x, t) ∂t ∂x ∞ da0 dan = ρc + pa0 (t) + ρc (t) + dt dt n=1 ρc
−
ρc
dbn (t) + dt
κn2 π2 nπx + p an (t) cos + 2
κn2 π 2 nπx + p bn (t) sin 2
46
CHAPTER 6. HEAT FLOW AND DIFFUSION Following the derivation in Section 6.3.3, we find a 0 , a1 , . . ., b 1 , b2 , . . . by solving the IVPs da0 + pa0 = c 0 (t), a0 (t0 ) = p 0 , dt κn2 π 2 + p an = c n (t), an (t0 ) = p n , n = 1, 2, . . . , 2 ρc
dan + ρc dt ρc
dbn + dt
κn2 π2 + p bn = d n (t), bn (t0 ) = q n , n = 1, 2, . . . . 2
(b) The eigenvalues of the spatial operator are λ 0 = p and λn = Therefore, λ n
κn2 π 2 + p, n = 1, 2, . . . . 2
≥ λ for all n, and all the eigevalues are positive if and only if p > 0. 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) = a 0 +
∞
nπx nπx + bn sin
an cos
n=1
,
where 1 f (x) dx, 2 − 1 nπx an = f (x)cos dx, n = 1, 2, . . . , −
a0 =
bn =
−
For any function g : ( , )
1 nπx f (x)sin dx, n = 1, 2, . . . . −
→ R, if g is even, then
g(x) dx = 2
−
g(x) dx,
0
while if g is odd, then
g(x) dx = 0.
−
Since f is even,
1 1 a0 = f (x) dx = f (x) dx. 2 − 0 Also, f (x)cos(nπx/) is an even function of x (the product of even functions is even), and therefore
1 nπx 2 an = f (x)cos dx = −
f (x)cos
0
nπx dx, n = 1, 2, . . . .
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 1 nπx bn = f (x)sin dx = 0, n = 1, 2, . . . . −
Thus the full Fourier series reduces to
∞
f (x) = a 0 +
n=1
where
an cos
nπx ,
1 2 nπx a0 = 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).
47
6.4. FINITE ELEMENT METHODS FOR THE HEAT EQUATION
6.4
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 must hold, and n
(x, Gx) =
n
(Gx)i xi =
i=1
n
n
Gij xj xi =
i=1 j=1
n
(uj , ui )xj xi
=
i=1 j=1
n
n
xj uj , ui
i=1
n
xj uj ,
j=1
xi ui
i=1
2
n
=
ui
j=1
n
=
. Then (x, Gx) = 0
n
∈ R
xj uj
.
j=1
Thus Gx = 0 implies that the vector n j=1 xj uj is the zero vector. But, since u1 , u2 , . . . , un is linearly independent, this in turn implies that x 1 = x 2 = = x n = 0, that is, that x = 0. Since the only vector x R n satisfying Gx = 0 is the zero vector, G is nonsingular.
3.
{
···
(a) M =
2/9 1/18
(b) The system of ODEs is M dα = dt
1/18 2/9
, K =
−3
6 3
−
6
1 11cos(t) , f (t) = 162 11cos(t)
−Kα + f (t), that is,
2 dα1 1 dα2 + 9 dt 18 dt 1 dα1 2 dα2 + 18 dt 9 dt 5.
}
= =
−6α + 3α + 11cos(t) , 162 11cos (t) 3α − 6α + . 162 1
2
1
2
(a) Write ρ 1 = 8.97, ρ 2 = 7.88,
ρ(x) =
0 < x < 50, 50 < x < 100,
ρ1 , ρ2 ,
and similarly for c(x) and κ(x). Then the IBVP is ∂u ρ(x)c(x) ∂t
−
∂ ∂x
=
0, 0 < x < 100, t > 0,
=
5, 0 < x < 100,
u(0, t)
=
0, t > 0,
u(100, t)
=
0, t > 0.
∂u κ(x) ∂x u(x, 0)
(b) The mass matrix M is tridiagonal and symmetric, and its nonzero entries are M i,i+1 = and M ii =
ρ1 c1 h , 6 ρ2 c2 h 6 ,
2ρ1 c1 h , 3 (ρ1 c1 +ρ2 c2 )h , 3 2ρ2 c2 h , 3
i = 1, 2, . . . , i =
n n 2, 2
n 2
− 1,
+ 1, . . . , n
i = 1, 2, . . . , n2 i = i =
n , 2 n + 1, n2 2
− 1,
− 1,
+ 2, . . . , n
− 1.
The stiffness matrix K is tridiagonal and symmetric, and its nonzero entries are K i,i+1 = and K ii =
− −
κ1 , h κ2 , h
2κ1 , h κ1 +κ2 , h 2κ2 , h
i = 1, 2, . . . , n2 i =
n n , 2 2
+ 1, . . . , n
i = 1, 2, . . . , i = i =
− 1,
n , 2 n + 1, n2 2
− 1,
n 2
− 1,
+ 2, . . . , n
− 1.
.
∈
48
CHAPTER 6. HEAT FLOW AND DIFFUSION 7. The solution is u(x, t) =
∞
nπx , 100
− an (t)sin
n=1
where 400(2 + ( 1)n ) an (t) = κ2 π 7 n7
−
4
10 ρc e
2 π2 − κn t 104 ρc
2 2
1 + κn π 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).
9.
(a) The IBVP is ∂u ρ(x)c(x) ∂t
−
∂ ∂x
=
0, 0 < x < , t > 0,
=
5, 0 < x < ,
u(0, t)
=
0, t > 0,
u(, t)
=
0, t > 0.
∂u κ(x) ∂x u(x, 0)
(b) The weak form is to find u satisfying
0
∂u ∂u dv ρ(x)c(x) (x, t)v(x) + κ(x) (x, t) (x) dx = 0, t > 0, ∂t ∂x dx
∀ v ∈ V.
(c) The temperature distribution after 120 seconds is shown in Figure 6.5. t = 120 (seconds) 5 4.5 4 3.5 e r 3 u t a r e2.5 p m e 2 t
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.
49
6.5. FINITE ELEMENTS AND NEUMANN CONDITIONS
6.5
Finite elements and Neumann conditions
1. The BVP to be solved is 2
−k ddxu 2
du (0) dx du (100) dx
=
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
4 e r u 3 t a r e p m 2 e t 1
0
−1 0
20
40
60
80
100
x
Figure 6.6: The computed steady-state temperature distribution from Exercise 6.5.1. 3.
(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 2
−κ ddxu 2
du (0) dx du k (100) dx k
=
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 e r u4 t a r e p m 3 e t 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 n
˜ u = u K
·
n
n
˜ ij uj ui = K
n
i=0 j=0
n
a (φj , φi ) uj ui
=
i=0 j=0
n
a
i=0
uj φj , φi
j=0
n
=
a
n
uj φj ,
j=0
=
ui
ui φi
i=0
a(v, v).
Therefore, ˜ u = 0 K
⇒ u · K ˜ u = 0 ⇒ a(v, v) = 0. ˜ and a(v, v) = 0. By definition of a(·, ·), this is equivalent to (b) Suppose v ∈ V
dv (x) dx
κ(x)
0
2
dx = 0.
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 is the zero function. Therefore, v is a constant function. ˜ u = 0, it follows that v(x) = n ui φi (x) is a constant function, where u0 , u2 , . . . , un are the (c) Thus, if K 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 u 0 , u1 , . . . , un , so we see that u = C uc , which is what we wanted to prove.
7.
ˆ = v (a) Define 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: ˆ and then integrate by parts. multiply the differential equation by an arbitrary test function v V ,
∈
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 du k() ()v() dx
−
d du k(x) (x) + f (x) v(x) dx = 0 for all v dx dx
0
∈ Vˆ .
ˆ this implies that the differential equation ⊂ V ,
Since V
d du k(x) (x) + f (x) = 0, 0 < x < dx dx holds, and we then have k() Choosing any v
du ()v() = 0 for all v dx
∈ Vˆ .
∈ V ˆ with v() = 0 shows that the Neumann condition holds at x = .
9. The temperature distribution after 300 seconds is shown in Figure 6.8.
51
6.5. FINITE ELEMENTS AND NEUMANN CONDITIONS t=300 (seconds) 8 7 6 5
e r u t a r e4 p m e t
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 u Ku
·
⇒ a(v, v) = 0,
−1 where v = n 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 u.
− ct) + ψ(x + ct)). Figure 7.1 shows three snapshots of
−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 join 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 u ∂t 2
−c
2 2 ∂ u ∂x 2
−∞ < x < ∞, t > 0, u(x, 0) = ψ(x), −∞ < x < ∞, ∂u = 0, −∞ < x < ∞, ∂t where the support 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 ψ(x − ct) = 0 if any of the following conditions hold: • x ≤ a; • a < x < b and t ≥ − (since ψ (x − ct) = 0 after the trailing edge of the wave passes x ); • x ≥ b and (t ≤ − or t ≥ − ) (since ψ(x − ct) = 0 before the leading edge of the wave reaches x and = 0,
1
1
x1 a c
1
1
x1 b c
1
x1 a c
1
1
after the trailing edge of the wave passes x 1 ).
Similarly, ψ(x + ct) is a left-moving wave, and therefore ψ (x1 + ct) = 0 if any of the following is true:
• x ≥ b; • a < x < b and t ≥ • x ≤ a and (t ≤ − 1
1
1
b x1 ; c a x1 or t c
−
≥
b x1 ). c
−
Since u(x1 , t) is guaranteed to be zero if both ψ (x1 of the following is true:
− ct) and ψ (x + ct) are zero, we see that u(x , t) = 0 if any
1
1
1
a x1 b x1 or t ); c c x1 a b x1 and t ; c c x1 b x1 a or t ). c c
• x ≤ a and (t ≤ − • a < x < b and t ≥ • x ≥ b and (t ≤ − 1
1
−
≥
−
≥
≥
−
−
7. Consider the IBVP ∂ 2 u ∂t 2
2 2 ∂ u ∂x 2
−c
= 0, 0 < x < , t > 0,
u(x, 0) = ψ(x), 0 < x < , ∂u = 0, 0 < x < , ∂t u(0, t) = 0, t > 0, u(, t) = 0, t > 0, where the support of ψ is [a, b] and 0 < a < b < . If we define 1 u(x, t) = (ψ(x 2
− ct) + ψ(x + ct)) ,
then we know that u satisfies the same PDE and initial condition as does the solution of the above IBVP. Also, since the support 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 < t f = 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 c n = 0 in Example 7.4: u(x, t) =
∞
n=1
bn cos
cnπt 25
sin
A graph analogous to Figure 7.6 from the text is given in Figure 7.3.
nπx . 25
{
−
}
55
7.2. FOURIER SERIES METHODS FOR THE WAVE EQUATION
t n e m e c a l p s i d
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
∞
an (t)sin
bn
− c (2n50−c1) π
√
− √ −
u(x, t) =
n=1
(2n
− 1)πx 50
,
where an (t) = and bn =
8
2
n
2
2 2
cos
c(2n
− 1)πt
50
+
502 cn c2 (2n 1)2 π2
−
2sin(nπ/2) 2cos(nπ/2) + ( 1)n 4000 , cn = . 2 2 5(2n 1) π (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 1 t n 0.5 e m e c 0 a l p s i d−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
∞
+
cnπt 25
bn cos
n=1
cos
nπx , 25
where bn =
2 (2cos(nπ/2) 1 5n2 π2
− − (−1)
n
)
, n = 1, 2, 3, . . . .
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) =
∞
n=1
an (t)sin
(2n
− 1)πx 2
,
56
CHAPTER 7. WAVES 0.6
0.4 t n e
m e c a l p s i d 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
=
4cn 2 c (2n 1)2 π 2
− − 2(−1) , 125(2n − 1) π bn
n+1
cos
c(2n
− 1)πt 2
+
4cn , 1)2 π2
c2 (2n
−
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
2
x 10
t n e m e c1 a l p s i d
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 account 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 f 0 = 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 f 0 = 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
·
c 1 = f 0 = 2(0 + ∆) 2(0 + ∆)
k∆ 0
m 0 +∆
=
√ k∆
2
m0 (0 + ∆)
Some straightforward algebra yields ∆ =
4m20 f 02 . = 8.0586. k 4m0 f 02
Thus the string must be stretched to 48.0586 cm.
−
⇒
k∆ = f 02 . 4m0 (0 + ∆)
57
7.3. FINITE ELEMENT METHODS FOR THE WAVE EQUATION
7.3
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. 0.2
t n e m e c a l p s i d
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 24cm 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 ∂ 2 u ∂t 2
2 2 ∂ u ∂x 2
−c
=
0, 0 < x < 50, t > 0,
u(x, 0) ∂u (x, 0) ∂t u(0, t)
=
0, 0 < x < 50,
=
γ (x), 0 < x < 50,
=
0, t > 0,
u(50, t)
=
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.)
·
≤ ≤ 0.1
t n e m e c a l p s i d
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
∂ 2 u ∂u dv ρ(x) 2 (x, t)v(x) + k(x) (x, t) (x) dx = ∂t ∂x dx
0
·
0
f (x, t)v(x) dx
58
CHAPTER 7. WAVES 1.2
1
0.8 t n e
0.6 m e c a l p s 0.4 i d 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
≥t
0
and v
˜ where ∈ V ,
˜ = 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), 2
dt
except that the entries in the mass matrix are now ˜ ij = M
ρ(x)φj (x)φi (x) dx, i,j = 1, 2, . . . , n ,
0
and the entries in the stiffness matrix are now ˜ ij = K
k(x)
0
dφj dφi (x) (x) dx, i,j = 1, 2, . . . , n . dx dx
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) =
∞
an (t)sin
n=1
Following the usual procedure, a n must solve the IVP d2 an cnπ + 2 dt
2
nπx .
dan (0) = 0, dt
an = c n (t), an (0) = 0,
where c n (t), n = 1, 2, . . ., are the Fourier sine coefficients of the right-hand-side function f (t) = sin(2πωt): 2 cn (t) =
sin (2πωt)sin
0
nπx 2 (1 dx = nπ
=
− (−1)
n
)sin(2πωt)
0,
if n is odd, if n is even.
4 sin(2πωt), nπ
Therefore, a n (t) = 0 if n is even, while, if n is odd, an (t) = cnπ
t
sin
0
cnπ (t
− s)
cn (s) ds.
If n is odd and ω = cn/(2), we obtain
42 (cn sin(2πωt) 2ω sin an (t) = cn2 π2 (c2 n2 42 ω2 )
− −
cnπt
)
.
59
7.4. RESONANCE If n is odd and ω = cn/(2), we obtain an (t) =
cnπt
−
2( sin
cnπt cos 2 3 c n π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 a n (t) is larger (due to the factor of c 2 n2 42 ω 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 a n (t) grows without bound as t , and resonance is observed.
−
→∞
3. The experiment described in this problem is modeled by the IBVP ∂ 2 u ∂t 2
−c
2 2 ∂ u ∂x 2
= 0, 0 < x < 1, t > 0,
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 c 2 = 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
−
∂ 2 w ∂t 2
−
2 2 ∂ w ∂x 2
4π 2 ω 2 B sin(2πωt)(x 1), 0 < x < 1, t > 0, k w(x, 0) = 0, 0 < x < 1, 2πωB ∂w (x, 0) = (x 1), 0 < x < 1, ∂t k ∂w k (0, t) = 0, t > 0, ∂x w(1, t) = 0, t > 0.
−c
=
−
−
−
We solve this IBVP by the usual method. We write u(x, t) =
∞
an (t)sin(nπx).
n=1
The Fourier coefficients a n are found by solving the IVP d2 an dan + c2 n2 π2 an = c n (t), an (0) = 0, (0) = b n , 2 dt dt where
1
cn (t) = 2
0
and
4π2 ω2 B sin(2πωt)(x k bn =
1
− 4πωB k
0
It follows (see Section 4.2.3 from the text) that bn 1 an (t) = sin(cnπt) + cnπ cnπ
t
0
sin(cnπ(t
−
an (t) =
− 1)sin(nπx) dx = 4ωB . kn
4ωB 1 s))cn (s) ds = sin (cnπt) + 2 kcn π cnπ
If ω = cn/2, then
(x
2
− 1)sin(nπx) dx = − 8πωkn B sin(2πωt)
4ωB sin (cnπt) kcn2 π
2
t
sin(cnπ(t
0
sin(2πωt) − 2ω sin(cnπt)) − 8ω B(cnckn , π(c n − 4ω ) 2
2 2
2
− s))c
n (s) ds.
60
CHAPTER 7. WAVES while if ω = cn/2, an (t) =
4ωB sin(cnπt) kcn2 π
2
4ω B − ckn π 2
sin(cnπt) cnπ
− t cos(cnπt)
.
(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 terms 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
1
0.5
1
−3
x 10
5
0
−5 0
0.5
x 10
0
0.5
1
−5 0
Figure 7.11: Absence of resonance in Exercise 7.3.3(b).
61
7.5. FINITE DIFFERENCE METHODS FOR THE WAVE EQUATION 5. We solve the IBVP by the usual method, taking
∞
u(x, t) =
an (t)sin(nπx),
n=1
where the Fourier coefficient a n solves the IVP d2 an dan + c2 n2 π2 an = c n (t), an (0) = 0, (0) = 0. dt2 dt Here c n is the Fourier coefficient of the right-hand-side function f (x, t): 1
cn (t) = 2
13/5
f (x, t)sin(nπx) dx = 2 sin (2πωt)
0
− sin(nπx) dx
3/5
=
2 sin (2πωt) nπ
3nπ 5
cos
cos
13nπ 20
.
It follows that t 1 sin(cnπ(t cnπ 0 2 cos 3nπ cos 5 = cn2 π2
an (t) =
− −
s))cn (s) ds
13nπ 20
(cn sin(2πωt) 2ω sin(cnπt)) . π(c2 n2 4ω 2 )
− −
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 dφ φ(x + ∆x) (x) = dx ∆x or
dφ (x) dx
2
− φ(x) − 1 d φ (c)∆x 2 dx2 2
− φ(x) = 1 d φ (c)∆x. − φ(x + ∆x) ∆x 2 dx 2
Since
d2 φ d2 φ (c) (x) dx2 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, φ(x + ∆x) = φ(x) + φ(x
− ∆x) = φ(x) −
dφ (x)∆x + dx dφ (x)∆x + dx
where c 1 is between x and x + ∆x and c 2 φ(x + ∆x) + φ(x
−
1 d2 φ 1 d3 φ 2 (x)∆x + (x)∆x3 + 2 dx2 6 dx3 1 d2 φ 1 d3 φ 2 (x)∆x (x)∆x3 + 2 3 2 dx 6 dx
1 d4 φ (c1 )∆x4 , 24 dx4 1 d4 φ (c2 )∆x4 , 4 24 dx
− is between x and x − ∆x. Adding these two equations yields
1 d 2 φ ∆x) = 2φ(x) + 2 (x)∆x2 + dx 24
d4 φ d4 φ (c ) + (c2 ) ∆x4 . 1 dx4 dx4
Solving this equation for the second derivative yields d2 φ φ(x + ∆x) (x) = 2 dx
− 2φ(x) + φ(x − ∆x) − ∆x2
1 24
d4 φ d4 φ (c ) + (c2 ) ∆x2 . 1 dx4 dx4
Since
d4 φ d4 φ d4 φ (c ) + (c ) 2 (x) 1 2 dx4 dx4 dx4 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 . u(j) = u(xi , tj ), i = 1, 2, . . . , k, j = 1, 2, . . . , n i (using the same notation as in the text). We apply the usual 2-2 scheme at each (xi , tj ), i = 1, 2, . . . , k 1, resulting in equation (7.41) (page 301 in the text) for computing u(j+1) , i = 1, 2, . . . , k 1. To compute u(j+1) , i k we approximate the Neumann condition ∂u (xk , tj ) = 0 ∂x by u(j) u(j) (j) (j) k+1 k−1 =0 uk+1 = u k−1 . ∆x
−
−
−
⇒
We then obtain equation (7.48) for u (j+1) . Finally, we use (7.44) to compute u (1) 1, and (7.41) i , i = 1, 2, . . . , k k (1) with i = k and ψ (xi+1 ) replaced by ψ(xi−1 ) to compute u k . 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 verify 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) = u 0 (x), 0 < x < 1, u(0, t) = φ(t), t > 0. Characteristics passing through (s, 0), 0 < s < 1, have the form
⇔
x = 4τ + s, t = 3τ
This characteristic can also be written as t = (3/4)(x
τ =
t , s = x 3
− 34 t.
− s). On such characteristics,
−
4 t . 3
u(x, t) = v(s, τ ) = u 0 (s) = u 0 x Characteristics through (0, s), s > 0, have the form
⇔
x = 4τ, t = 3τ + s
τ =
x , s = t 4
− 34 x.
This characteristic can also be written as t = s + (3/4)x. On such characteristics,
−
u(x, t) = v(s, τ ) = φ(s) = φ t
3 x . 4
The characteristic dividing the domain into two is t = (3/4)x. Thus u(x, t) =
4 t 3 3 x 4
− − u0 x φ t
, ,
t 34 x, t > 34 x.
≤
(b) Now consider the initial/boundary conditions
u(x, 0) = u 0 (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 passes through a point (1, t0 ), 0 < t0 < 3/4. These means that the solution is over-determined on part of the domain, and u 0 , ψ 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) = u 0 (x), 0 < x < 1, u(0, t) = φ(t), t > 0. The solution is u(x, t) =
t < c−1 x, t > c−1 x.
u0 (x ct), φ(t c−1 x),
−
−
We assume that u 0 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 t 0 = c−1 x0 . We have u0 (x
− ct) → u (x − ct ) = u (0) as (x, t) → (x , t ), φ(t − c− x) → φ(t − c− x ) = φ(0) as (x, t) → (x , t ), 0
1
0
0
0
1
0
0
which shows that u is continuous if and only if u 0 (0) = φ(0).
0
0
0
0
65
8.2. FIRST-ORDER QUASILINEAR PDES (b) Next, we have ∂u (x, t) = ∂x and du0 (x ct) dx dφ c−1 (t c−1 x) dt
−
du0 (x dx 1 dφ (t dt
−c−
t < c−1 x, t > c−1 x,
− ct),− − c x), 1
du (x − ct ) = (0) as (x, t) → (x , t ), − → du dx dx dφ (t − c− x ) = −c− (0) as (x, t) → (x , t ), − → −c− dφ dt dt 0
0
0
0
1
1
0
0
0
1
0
0
0
Therefore, ∂u/∂x is continuous if and only if du0 (0) = dx
−c−
1 dφ
dt
(0).
Similar reasoning shows that this condition also ensures that ∂u/∂t is continuous. 7. Consider the IVP ∂u ∂ u + = x + y, (x, y) ∂y ∂x u(x, x) = x 2 , x
−
2
∈R ,
∈ 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) = s 2 . ∂t
−
The first two equations yield x = s + t, y = and therefore
−s + t ⇔
∂v = 2t, v(s, 0) = s 2 ∂t
Thus, u(x, y) =
8.2
2
s =
⇒
x
2
y
2
+
2
v(s, t) = s 2 + t2 .
x+y 2
2
− x
− y , t = x + y ,
=
1 2 (x + y 2 ). 2
First-order quasilinear PDEs
1. Consider the following IVP: x
∂u ∂u +y ∂x ∂y
− (x + y)u = 0, 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 = y 0 x. Thus the characteristic through (1, y0 ) is the line with slope y0 passing 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: t v(s, t) = φ(s)e(1+s)(e −1) .
Changing variables, we obtain
y x+y−1−y/x e . 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) ∂s ∂y (s, 0) ∂s
∂x (s, 0) ∂t ∂y (s, 0) ∂t
=
1 0
φ(s) φ(s)
= φ(s),
we see that the initial curve is noncharacteristic if and only if φ(s) = 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 s =
x+
− x2 2
4y
, t =
x
−
− x2 2
4y
.
67
8.3. BURGERS’S EQUATION Since v (s, t) = s, it follows that u(x, y) = This solution is defined for x > 0, y
x+
2
≤ x /4.
− x2 2
4y
.
(b) The characteristic indexed by s passes 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 = x 0 (x x0 ).
−
−
−
(c) Since the slopes of the characteristics increase as x 0 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 = x 0 (x x0 ) and y = x 1 (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 = x 2 /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) = u 0 (x) =
1 . 1 + x2
We have seen that the characteristics are the lines x = x 0 + u0 (x0 )t
⇔
t = u 0 (x0 )−1 (x
− x ). 0
(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 passing through any such point. Since the slope (in the x, t plane) of every characteristic is positive, the characteristic through (x0 , 0) passes through (x, t) only if x 0 < x < 0 and t = (1 + x20 )(x
3 0
2 0
− x ) ⇔ x − xx + x + (t − x) = 0. If we define ψ(x ) = x − xx + x + (t − x), then it is easy to verify that ψ(0) > 0, ψ(x ) < 0 for x 3 0
0
2 0
0
0
0
0
0
sufficiently large and negative, and ψ (x0 ) > 0 for all x 0 < 0. It follows that ψ(x0 ) = 0 has a unique solution x0 < 0, and hence there is a unique characteristic passing 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 + 3x20 (1 + x20 )2 , t = , x0 > 0. 2x0 2x0
The corner of the caustic can be found by setting the derivatives of x and t (with respect to x 0 ) to zero and solving, which yields 1 8 x0 = , x = 3, t = . 3 3 3 Moreover, we can solve 1 + 3x20 x = 2x0 for x 0 in terms of x to get x x2 3 x0 = . 3 Substituting into the equation for t yields the two branches of the caustic:
√
√
√
± √ −
3 1+ t1 =
x
√
−
x2 3
−3
√ 2(x − x − 3) 2
2
2
3 1+ , t2 =
√
x+
x2 3
−3
√ 2(x + x − 3) 2
2
2
, x
≥ √ 3
68
CHAPTER 8. FIRST-ORDER PDES AND THE METHOD OF CHARACTERISTICS (where t 2 > 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
3. Consider the IVP ∂u + (1 ∂t
≥ √
3 and (t < t 1 or t > t 2 .
x
− u) ∂u = 0, −∞ < x < 0, t > 0, ∂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, ∂τ 1 ∂v = 0, v(s, 0) = . ∂τ s
−
Solving yields x(s, t) = s +
s
Setting s = x 0 , we obtain
− 1 τ, s
x0 1 t x0 These are the equations of the characteristics. x = x 0 +
−
t(s, t) = τ, v(s, t) =
⇔
(b) We can solve x = s +
t =
s
x0 (x x0 1
−
1 . s
− x ). 0
− 1t s
for s as follows: x = s +
s
− 1t ⇒ s
2
− x)s − t = 0 ⇒ s = x − t ± (t − x)
s
+ (t
⇒ s = x
− − − 2 (t 2
t
2
+ 4t
x)2 + 4t
,
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) =
x
2 (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 x0
−1
(x
− x ), t = x x− 1 (x − x ) 1
0
1
1
for (x, t) yields (x, t) = (x0 Notice that t < 0, as expected.
− x x + x , −x x ). 0 1
1
0 1
5. Consider the IVP ∂u ∂u + a(u) = 0, 0.
−∞
∞, t > 0,
The characteristics are defined by x(s, τ ) = s + a(φ(s))τ, t(s, τ ) = τ , v(s, τ ) = φ(s).
69
8.3. BURGERS’S EQUATION We can write the characteristic passing through (x0 , 0) as x = x 0 + a(φ(x0 ))t
⇔
t =
−
x x0 . a(φ(x0 ))
Assuming that a(φ(x)) is an decreasing function of x, we know that characteristic passing through (x0 , 0) will intersect with the characteristic passing through (x1 , 0) if x 1 > x0 . We can compute the p oint of intersection as follows: x = x 0 + a(φ(x0 ))t, x = x1 + a(φ(x1 ))t
⇒ x + a(φ(x ))t = x + a(φ(x ))t ⇒ t = a(φ(x x)) −− xa(φ(x )) = − 0
0
1
0
1
1
1
1
a(φ(x1 )) a(φ(x0 )) x1 x0
0
− −
.
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 da (φ(x0 )) dφ (x0 ) du dx
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 (s, τ ) ∂s ∂t (s, τ ) ∂s
∂x (s, τ ) ∂τ ∂t (s, τ ) ∂τ
=
1+
da (φ(s)) dφ (s)τ du dx
0
a(φ(s)) 1
= 1+
da dφ (φ(s)) (s)τ. du dx
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) = −k x
u(x) =1 x
−
(the uncompressed length of the part of the bar between x = s and x = is u(x) =
The Green’s function is G(x; s) =
s , k x , k
− −
0 < x < s, s < x < .
s , k x , k
− −
0 < x < s, s < x < .
(b) The Green’s function is
G(x; s) =
dz . k(z ) max{x,s}
(c) If k is constant, then G(x; s) =
− max{x, s} = k
s , k x , k
− −
0 < x < s, s < x < .
(d) ln 2 1 x 2 x + 2 4 4 2 5. By direct integration, we obtain the following solution to (9.11): u(x) =
− ln (12+ x) .
− −
x
u(x) = p
0
{ }
Since x = min x, for x
ds . k(s)
∈ [0, ], we can write min x,s
u(x) = p
0
{ } ds k(s)
that is, u(x) = g(x; ) p.
71
,
− x). Therefore we obtain
72
CHAPTER 9. GREEN’S FUNCTIONS 2 7. Let the operators K : C m [0, ]
2 m [0, ]
→ C [0, ] and M : C [0, ] → C Ku =
−
d dx
du k(x) dx
be defined by
and
(M f )(x) =
G(x; y)f (y) dy,
0
where
min x,y
G(x; y) =
{
} 1 k(z )
0
We wish to prove directly that (M K )u = u for all u
2 m [0, ].
∈ C
x
(M f )(x) =
G(x; y)f (y) dy +
x
0
=
0
0
G(x; y)f (y) dy
x
y
x
Given the definition of G, we will write
0
=
dz.
x
z
1 f (y) dzdy + k(z ) 1 f (y) dydz + k(z )
y
1 f (y) dzdy k(z )
1 f (y) dydz, k(z )
x
0
x
0
x
where the last step follows from changing the order of integration. Therefore, x
(M Ku)(x) =
x
x
− − − − − − − 0
z
1 d k(z ) dx
du k(y) (y) dydz dx
0
x
1 du du = k(x) (x) k(z ) (z ) dz dx dx 0 k(z ) x 1 du = k(z ) (z ) dz dx 0 k(z ) x du = (z ) dz = u(x) 0 dx
1 d du k(y) (y) dydz dx x k(z ) dx x 1 du du k() () k(x) (x) dz dx dx 0 k(z )
−
(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 x
≤ y ≤ x, ≤ y ≤ 1.
3. The Green’s function is G(x; y) = where v 1 (x) = eθx
θx
− e−
, u 2 (x) = eθx
−e
−
kv1 (y)v2 (x),
−kv (x)v (y), 1
2
0 x
≤ y ≤ x, ≤ y ≤ 1,
θ(2 x)
− , and k =
1 2(e2θ
− 1)θ .
5. Let G = G(x, y) be any continuous function of two variables, a by
≤ x ≤ b, a ≤ y ≤ b, and define M : C [a, b] → C [a, b]
b
(M f )(x) =
a
G(x, y)f (y) dy.
73
9.2. GREEN’S FUNCTIONS FOR BVPS IN ODES: THE SYMMETRIC CASE
∈ [a, b]. We then have
Suppose G is symmetric in the sense that G(y, x) = G(x, y) for all x, y b
(M f , g)L2 =
(M f )(x)g(x) dx
a b
=
b
G(x, y)f (y) dy g(x) dx
a b
= = =
a b
G(x, y)f (y)g(x) dydx
a b
a b
a b
a
G(x, y)f (y)g(x) dxdy (changing the order of integration) b
f (y)
a b
=
G(x, y)g(x) dx dy
a b
f (y)
a b
=
G(y, x)g(x) dx dy (since G is symmetric)
a
f (y)(M g)(y) dy = (f , M g)L2 .
a
Thus M is a symmetric operator with respect to the L 2 (a, b) inner product. 7. Define H : R
→ R by H (x) =
We wish to prove that
holds for all φ
∞
δ (x)φ(x) dx =
−∞
0, 1.
x < 0, 0 < x
∞
−
H (x)
−∞
dφ (x) dx dx
∈ D(R). That is, we wish to prove that if φ ∈ D(R), then ∞
H (x)
−∞
dφ (x) dx = dx
−φ(0).
This is a direct calculation. Choose b > 0 sufficiently large that φ(x) = 0 for all x
∞
H (x)
−∞
∞ dφ
dφ (x) dx = dx
0
dx
b
(x) dx =
0
dφ (x) dx = φ(b) dx
≥ b. Then
− φ(0) = −φ(0),
as desired. 9. Let G be the Green’s function for the BVP d dx
−
du P (x) 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
Then G(x; y) =
− −
v1 (y)v2 (x) , P (x)W (x) v1 (x)v2 (y) , P (x)W (x)
a x
where v 1 and v 2 are certain solutions of the homogeneous ODE
−
d dx
du P (x) dx
≤ y ≤ x, ≤ y ≤ b,
+ R(x)u = 0, a < x < b.
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
−
du P (x) (x)v(x) dx
−
b
∈ b
+
a
a
b
du dv P (x) (x) (x) + R(x)u(x)v(x) dx = dx dx
a
δ (x
− y)v(x) dx = v(y).
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
du du (y ) = lim (x) = dx x→y − dx
(y)v (y) − vP (y)W (y) ,
du du (y+) = lim (x) = dx x→y + dx
(y)v (y) − vP (y)W (y) .
−
2
1
1
2
We now substitute u into the left-hand side of the variational equation to obtain
−
du = P (a) (a)v(a) dx b
+
−
y
b
b
du P (x) (x)v(x) dx
+
a
a
du dv P (x) (x) (x) + R(x)u(x)v(x) dx dx dx
du P (b) (b)v(b) + dx
y
a
du dv P (x) (x) (x) + R(x)u(x)v(x) dx dx dx
du dv P (x) (x) (x) + R(x)u(x)v(x) dx. dx dx
The integrands of both integrals on smooth on the intervals of integration, and hence integration by parts is valid: b
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 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 y b d du d du + 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 a y
− − −
du = P (y) (y− )v(y) dx du = P (y) (y− )v(y) dx
− −
− − − b
du P (y) (y+ )v(y) + dx du P (y) (y+ )v(y), dx
d dx
a
−
du P (x) (x) + R(x)u(x) v(x) dx dx
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
−
= P (y)v(y) = P (y)v(y) 2 since W (y) = v 1 (y) dv (y) dx
−
b
− −
du P (x) (x)v(x) dx
+
a
du − (y ) dx
du dv P (x) (x) (x) + R(x)u(x)v(x) dx dx dx
a
du + (y ) dx
2 1 (y) dv (y)v2 (y) v1 (y) dv dx dx = v(y) P (y)W (y)
dv1 (y)v2 (y). dx
This completes the proof.
11. We wish to find a formula for the solution of
−
d dx
du P (x) dx
+ R(x)u = 0, a < x < b,
du (a) = v a , dx du β 1 u(b) + β 2 (b) = v b , dx
α1 u(a) + α2
75
9.3. GREEN’S FUNCTIONS FOR BVPS IN ODES: THE GENERAL CASE in terms of the Green’s function G for the BVP
−
d dx
du P (x) 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
We proceed as in Section 9.2.2; in particular, we apply (u,Lv)
−
dv (Lu,v) = P (a) u(a) (a) dx
−
du (a)v(a) dx
−
dv P (b) u(b) (b) dx
−
du (b)v(b) 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 u(y) = P (a) u(a) (a) dx
−
−
du (a)v(a) dx
dv P (b) u(b) (b) dx
−
du (b)v(b) . dx
We now apply the boundary condition satisfied by u and v: du 1 α (a) = ⇒ dx v − u(a), α α du 1 β ⇒ dx (b) = v − u(b), β β
du (a) = v a dx du β 1 u(b) + β 2 (b) = v b dx dv α1 v(a) + α2 (a) = 0 dx dv β 1 v(b) + β 2 (b) = 0 dx
α1 u(a) + α2
2
2
⇒ ⇒
dv (a) = dx dv (b) = dx
− −
1
a
b
2
1 2
α1 v(a), α2 β 1 v(b). β 2
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 α1 1 α 1 β 1 = P (a) u(a)v(a) va v(a) + u(a)v(a) P (b) u(b)v(b) α2 α2 α2 β 2
−
−
= β 2−1 P (b)v(b)vb
1
− α− P (a)v(a)v = β − P (b)G(b; y)v − α− P (a)G(a; y)v = β − P (b)G(y; b)v − α− P (a)G(y; a)v . 2 2
1
2
b
1
b
−
− −
−
1 β 1 vb v(b) + u(b)v(b) β 2 β 2
a
2 2
1
a
1
a
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
9.3
1
− α− P (a)G(y; a)v . 2
a
Green’s functions for BVPs in ODEs: the general case
1. The Green’s function is G(x; y) = where v 1 (x) = e x
−e
2x
, v 2 (x) = e x
−e
− −
v1 (y)v2 (x) , P (x)W (x) v1 (x)v2 (y) , P (x)W (x)
0 x
≤ y ≤ x, ≤ y ≤ 1,
2x 1
− , and P (x)W (x) is the constant 1 − e−1 . The solution to the BVP is 1
u(x) =
0
where w(x) = e −3x .
G(x; y)f (y)w(y) dy,
76
CHAPTER CHAPTER 9. GREEN’S GREEN’S FUNCTIO FUNCTIONS NS 3. The Green’s Green’s function function is G(x; y ) = where v1 (x) = e x sin(2x sin(2x), v2 (x) = e x (sin(2 sin(2x x) solution to the BVP is
− −
v1 (y )v2 (x) , P ( P (x)W ( W (x) v1 (x)v2 (y) , P ( P (x)W ( W (x)
≤ y ≤ x, ≤ y ≤ 1,
0 x
(tan (2)) (2)) cos cos (2x (2x)), and P ( constan antt 2 tan (2). The − (tan P (x)W ( W (x) is the const 1
u(x) =
G(x; y )f ( f (y)w(y) dy,
0
where w where w((x) = e −2x . 5. The Green’s Green’s function function is G(x; y ) =
ey 6 + 6x 6 x + 3x 3x2 + x3 , x e 6 + 6y 6 y + 3y 3 y2 + y3 ,
The solution to the BVP when f when f ((x) = x is
u(x) =
1 2
23 23 + x + 6x 6x2 + 2x 2 x3 2 2
≤ y ≤ x, ≤ y ≤ 1.
x
x 1
− 7e − . b
7. Let Let G be G be the Green’s function derived in this section, and let M be M be defined by (M (M f )(x )(x) = a G( G (x; y )f ( f (y ) dy. dy. (Here 2 maps complex C complex [ ] into complex C complex [ ].) We wish to prove that ( (M ) = ( ) for all f all M C [a, b C a, b M f , g) g w f , M gg)w f , g C [a, b]. In the following calculation, we use the facts that G that G is symmetric (G (G(x; y) = G( all x, y [a, b]) and G and G is G (y ; x) for all x, real-valued, so that G that G((x; y ) = G( G(x; y). We have b
(M f , g) g )w =
b
(M f )( f )(x x)g (x)w(x) dx = dx =
a
b
a b
b
G(x; y )f ( f (y )w(y )g (x)w(x) dydx
a
a
b
=
b
G(x; y )f ( f (y )w(y )g (x)w(x) dxdy
a
a
b
=
b
f ( f (y )
a
G(x; y )g (x)w(x) dx w(y ) dy
b
f ( f (y )
a
G(x; y )g (x)w(x) dx w(y ) dy
a
b
=
a
b
=
G(x; y )f ( f (y )w(y ) dy g(x)w(x) dx
a
=
b
f ( f (y )
a
G(y ; x)g (x)w(x) dx w(y ) dy
a
b
=
∈
∈
)(M g)(y )(y )w(y ) dy = f ( f (y )(M dy = (f , M gg))w .
a
(Notice how the key step was the interchanging of the order of integration.) This proves the desired result.
9.4
Introd Introduc ucti tion on to Green Green’s ’s func functio tions ns for IVPs IVPs
1. We wish to show that the solution of d2 u du + θ2 u = f = f ((t), u(0) = u = u 0 , (0) = v = v 0 2 dt dt is
∂G u(t) = (t; 0)u 0)u0 + G(t; 0)v 0)v0 + ∂t
where G(t; s) = We can write the proposed solution explicitly as cos (θt) u(t) = cos θt)u0 +
∞
G(t; s)f ( f (s) ds,
0
sin sin (θ(t s)) , θ
t > s, t < s.
t
sin(θ sin(θ(t θ
−
0,
sin(θt sin(θt)) v0 + θ
0
− s)) f (s) ds.
77
9.4. INTRODUCTION INTRODUCTION TO GREEN’S FUNCTIONS FUNCTIONS FOR FOR IVPS Then du (t (t) = dt
sin sin (θ(t θ sin(θt sin(θt))u0 + cos cos (θt) θt)v0 + θ
−
− t)) f (t) +
0
d2 u (t ( t) = dt2 =
−θ
2
−θ
2
0
cos(θt cos(θt))u0
cos(θ cos(θ(t
0
))f ((s) ds − s))f
− s))f ))f ((s) ds t
sin(θt))v + cos cos (θ(t − t))f ))f ((t) − − θ sin(θt 0
−
− θ sin(θt sin(θt))v + f ( f (t) 0
θ sin(θ sin(θ(t
0
θ sin(θ sin(θ(t
0
t
cos(θt cos(θt))u0
cos(θ cos(θ(t
0
t
= −θ sin(θt sin(θt))u + cos cos (θt) θt)v +
t
))f ((s) ds − s))f
− s))f ))f ((s) ds.
From these formulas, we can easily verify that u that u satisfies the ODE and the initial conditions: u(0) = cos cos (0)u (0)u0 + du (0) = dt
0
sin sin (0) v0 + θ
0
sin(0)u + cos(0)v cos(0)v + −θ sin(0)u 0
0
sin(θ sin(θ(t θ
− s)) f (s) ds = ds = u u , 0
0
cos(θ cos(θ(t
0
))f ((s) ds = − s))f ds = v v , 0
and d2 u (t ( t) + θ2 u(t) = dt2
−θ
2
t
cos(θt cos(θt))u0
sin(θt))v + f ( − θ sin(θt f (t) 0
t
2
cos(θt))u0 + θ sin(θt sin(θt))v0 + θ cos(θt
−
sin(θ(t θ sin(θ
0
= f ( f (t). 3. The Green’s Green’s function function is G(t; s) = 5. The Green’s Green’s function function is G(t; s) =
G(t; s) =
0)P 0 + P ( P (t) = G( G(t; 0)P
∞
−
))f ((s) ds − s))f
2(t 2(t s)
− e−
− , t > s,
0,
t < s.
02(t−s) e0.02(t , 0,
t > s, t < s,
02t 55.5e0.02t + 450 + 10t 10t G(t; s)f ( f (s) ds = ds = 55.
0
0.02t 02t
450e − 450e
= 450 + 10t 10t
9. Let the solution solution of d2 u du + a1 (t) + a0 (t)u = 0, 2 dt dt u(s) = 0, du (s (s) = f ( f (s) dt be u( u (t) = S (t; s)f ( f (s), and define the Green’s function G( G (t; s) by G(t; s) =
))f ((s) ds+ − s))f ds+
t > s, t < s.
0,
e−(t−s)
7. The Green’s Green’s function function is
and the solution to the IVP is
sin(2(t sin(2(t s)) , 2
sin(θ(t θ sin(θ
0
S (t; s), 0,
t > s, t < s,
We wish to show that G that G satisfies ∂ 2 G ∂G (t; s) + a1 (t) (t ( t; s) + a0 (t)G(t; s) = δ (t 2 ∂t ∂t G(t0 ; s) = 0, ∂G (t0 ; s) = 0, ∂t
− s),
0.02t 02t
394.5e − 394.
.
78
CHAPTER CHAPTER 9. GREEN’S GREEN’S FUNCTIO FUNCTIONS NS provided t provided t 0 < s. We notice that G that G((t, s) = H ( H (t
− s)S (t, s), where H where H is is the Heaviside function. We will show that 2
∂G (t, ( t, s) = H ( H (t ∂t
2
∂ G ∂ S − s) ∂S (t, ( t, s) and (t, s) = δ (t − s) + H (t − s) (t, ( t, s) ∂t ∂t ∂t (in the weak sense). sense). Let φ be any smooth function having compact support in [t [t , ∞), that is, any function for which φ which φ((t) = 0 for all t all t ∈ where t < a < b < ∞. We first show that (a, b), where t ∞ ∞ ∂S − G(t; s) dφ (t (t) dt = dt = H (t − s) (t, ( t, s)φ(t) dt, dt ∂t 2
2
0
0
t0
t0
which implies that
∂G (t, ( t, s) = H ( H (t ∂t
− s) ∂S (t, ( t, s) ∂t
in the weak sense. We have
∞
−
t0
dφ G(t; s) (t ( t) dt = dt = dt
∞
− −
H (t
t0
(t (t) dt − s)S (t, s) dφ dt
b
=
S (t, s)
s
=
dφ (t (t) dt dt b
b s
− S ( S (t, s)φ(t)| + b
∂S (t, ( t, s)φ(t) dt (by dt (by integration by parts) ∂t
s
∂S (t, ( t, s)φ(t) dt (since and φ((b) = 0) dt (since S S ((s, s) = 0 and φ s ∂t ∞ ∂S = ( t, s)φ(t) dt. H (t s) (t, ∂t t0 =
−
This proves the desired formula for ∂G/∂t. ∂G/∂t. Next, we will show that
∞ ∂G
−
t0
dφ (t ( t; s) (t ( t) dt = dt = φ φ((s) + ∂t dt
which implies that
∂ 2 G (t, s) = δ (t ∂t 2
∞
t0
2
H (t
( t, s)φ(t) dt, − s) ∂ ∂tS (t, 2
2
− s) + H (t − s) ∂ ∂tS (t, ( t, s) 2
in the weak sense. We have
−
∞ ∂G
t0
∂t
(t; s)
dφ (t (t) dt = dt = dt
∞
− −
t0 b
=
s
=
−
H (t
dφ (t, ( t, s) (t ( t) dt − s) ∂S ∂t dt
∂S dφ (t, ( t, s) (t ( t) dt ∂t dt
b
= φ( φ(s) +
s
= φ( φ(s) +
b
b
−
∂S (t, ( t, s)φ(t) ∂t
+
s
s
2
∂ 2 S (t, ( t, s)φ(t) dt (by dt (by integration by parts) ∂t 2
∂ S (t, ( t, s)φ(t) dt (since dt (since ∂t 2
∞
H (t
t0
s)
∂S (s, (s, s) ∂t
= 1 and φ and φ((b) = 0)
∂ 2 S (t, ( t, s)φ(t) dt. ∂t 2
This proves the desired formula for ∂ G2 /∂t 2 . Since S Since S ((t, s), as a function of t, t , satisfies the homogeneous version of the ODE, we see that ∂ 2 G ∂G (t; s) + a1 (t) (t; s) + a0 (t)G(t; s) = δ (t ∂t 2 ∂t as desired. Also,
∂G ∂S (t0 , s) = H ( H (t0 s) (t ( t0 , s) = 0 ∂t ∂t s < 0 by assumption. This completes the proof.
G(t0 , s) = H ( H (t0 since H since H ((t0
because t − − s) = 0 because t 0
− s),
− s)S (t , s) = 0,0 , 0
−
79
9.5. GREEN’S FUNCTIONS FOR THE HEAT EQUATION
9.5
Green’s functions for the heat equation
1. Given that
1 √ e− 2 πkt
S (x, t) = we have
∞
1 √ 2 πkt
S (x, t) dx =
−∞
√
x2 /(4kt)
∞
2
e−x
, /(4kt)
dx.
−∞
Making the change of variables u = x/(2 kt), we obtain
∞
∞ 2 1 1 S (x, t) dx = e−u du = π = 1. π −∞ π −∞
∈ R,
√ √
√
Also, for any x, y
∞
S (x
− y, t) dx =
−∞
∞
S (v, t) dv = 1.
−∞
3. For t = 660, six terms of the Fourier series are sufficient to give an accurate graph. For t = 61, even 100 terms are not sufficient to give a correct graph. 5. Following the example in Section 9.5.2, we obtain the following Green’s function: G(x, t; y, s) =
2 ρc
∞ e−k(2n−1)2 π2 (t−s)/(42 ) sin n=1 0,
(2n 1)πy 2
−
sin
(2n 1)πx 2
−
,
≤ s ≤ t,
0
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
0.02 t=120 t=240 t=360 t=480 t=600
0.015
e r u t a r e p
m e t
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 terms of the Fourier series were used to create these graphs. 7. We have
u(x, t) =
G(x, t; y, 0)φ(y) dy
0
=
0
=
∞
n=1
=
∞
n=1
2
∞
2
e−kn
π2 t/2
n=1
2
sin
0
2
bn e−kn
sin
nπy nπx sin φ(y) dy
2 2 2 nπy nπx φ(y) dy e−kn π t/ sin
π2 t/2
sin
nπx ,
80
CHAPTER 9. GREEN’S FUNCTIONS where 2
bn =
nπy φ(y) dy, n = 1, 2, 3, . . .
sin
0
are the Fourier sine coefficients of the function φ. We recognize u as the solution of the IBVP ∂u ∂t
2
− k ∂ ∂xu = 0, 0 < x < , t > 0, 2
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 u(x, t) = 2c
t
x+c(t s)
−
0
f (y, s) dy ds.
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 ∂u 1 (x, t) = 2c ∂t = =
1 2c 1 2c
x+c(t t)
−
{
x c(t t) x
− −
1 f (y, t) dy + 2c
f (y, t) dy +
x t
cf (x + c(t
0
1 ∂ 2 u (x, t) = (cf (x + c(t 2 ∂t 2c
1 2c
t
1 ∂u (x, t) = ∂x 2c ∂ 2 u 1 (x, t) = 2 2c ∂x
t
{
t
cf (x + c(t
0
df (x + c(t dx
f (x + c(t
0
t
0
− s), s) + cf (x − c(t − s), s)} ds
− s), s) + cf (x − c(t − s), s)} ds
− s), s) + cf (x − c(t − s), s)} ds,
− t), t) + cf (x − c(t −
0
cf (x + c(t
0
{
1 1 = (cf (x, t) + cf (x, t)) + 2c 2c c = f (x, t) + 2
t
{
t
2 df
c
0
t
1 t), t)) + 2c
dx
− s), s) −
(x + c(t
df (x dx
c
0
2 df
dx
(x + c(t
2 df
− s), s) − c
dx
ds,
− c(t − s), s)
(x
− s), s) − c
2 df
− c(t − s), s)
dx
(x
− c(t − s), s)
ds
− s), s) − f (x − c(t − s), s)} ds,
df (x + c(t dx
− s), s) −
df (x dx
− c(t − s), s)
ds.
Thus ∂ 2 u (x, t) ∂t 2
−
∂ 2 u c2 2 (x, t) ∂x
c = f (x, t) + 2
− 2c = f (x, t),
t
0 t
0
df (x dx
df (x + c(t dx
− s), s) −
− c(t − s), s)
ds
df (x + c(t dx
df (x − c(t − s), s) − s), s) − dx
ds
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)
ds
81
9.6. GREEN’S FUNCTIONS FOR THE WAVE EQUATION
∈
− | − |
for all x,y,t, s R such that c(t s) = x y . (The two expressions differ when c(t when H is defined as in the text: H (t) = 1 if t > 0 and H (t) = 0 if t 0.) We have
≤
− s) = |x − y| > 0, at least
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. |x − y|. This completes the proof. The reasoning can be reversed if we add the condition c(t − s) = 5. Let u be the solution to ∂ 2 u ∂t 2
2 2 ∂ u ∂x 2
−c
˜ t), = f (x,
−∞ < x < ∞, t > 0, u(x, 0) = 0, −∞ < x < ∞, ∂u (x, 0) = 0, −∞ < x < ∞, ∂t
˜ 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) = where f 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 ∂ 2 v ∂ 2 u (x, = ( x, t), t) ∂t 2 ∂t 2 ∂v ∂u (x, t) = ( x, t), ∂x ∂x ∂ 2 v ∂ 2 u (x, = ( x, t). t) ∂x 2 ∂x 2
− −
− − − − −
Therefore, ∂ 2 v (x, t) ∂t 2
2
−c
∂ 2 v (x, t) = ∂x 2
2
2
− ∂ ∂tu (−x, t) + c ∂ ∂xu (−x, t) ∂ u ∂ u =− ( −x, t) − c (−x, t) ∂t ∂x ˜ −x, t) = f (x, ˜ t) = −f ( 2
2
2
2
2
2
2
2
˜ odd). Thus v satisfies the same wave equation as does u. Moreover, it is obvious that v also satisfies (since f is the same initial conditions; hence v = u and the proof is complete. 7. We outline the derivation of the Green’s function for ∂ 2 u ∂t 2
2 2 ∂ u ∂x 2
−c
= f (x, t), 0 < x < , t > 0,
u(x, 0) = 0, ∂u (x, 0) = 0, ∂t ∂u (0, t) = 0, ∂x ∂u (, t) = 0, ∂x
0 < x < , 0 < x < , t > 0, t > 0.
82
CHAPTER 9. GREEN’S FUNCTIONS ˜ be the even, periodic extension of f , and define u We define f to ˜ to be the solution of ∂ 2 u ˜ ∂t 2
−c
2 ˜ 2 ∂ u ∂x 2
˜ t), = f (x,
−∞ < x < ∞, t > 0, ˜(x, 0) = 0, −∞ < x < ∞, u ∂ ˜ u (x, 0) = 0, −∞ < x < ∞. ∂t
From earlier results, we know that t
˜(x, t) = u
0
Then
∞ 1 ˜ s) dyds = 1 H (c(t − s) − |x − y |)f (y, 2c 2c −∞
t
∂ ˜ u (x, t) = ∂x
∂ ˜ u (x, 0) = ∂x
It follows that
∂ ˜ u (x, ) = ∂x ˜ even, we have f (c(t ˜ Since f is
˜ + c(t f (x
0
t
˜ f (c(t
0
t
˜ c(t f (
0
x+c(t s)
− −
− c(t − s), s)
˜ s) dy ds. f (y,
x c(t s)
ds.
ds,
˜ − c(t − s), s) − s), s) − f (
− s), s) = f (˜ −c(t − s), s), and hence
−
− s), s) − − − s), s)
˜ + c(t f (
0
− s), s) −
˜ f (x
t
ds.
∂ ˜ u (x, 0) = 0. ∂x Also, ˜ f (
˜ even) − c(t − s), s) = f (˜ − + c(t − s), s) (since f is ˜ + c(t − s), s) (since f ˜ is 2-periodic). = f (
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) =
∞
1 H (c(t 2c n=−∞
{
− s) − |x − y − 2n|) + H (c(t − s) − |x + y − 2n|)} .
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 du P (x) (x)v(x) dx
−
b
=
du −P (b) du (b)v(b) + P (a) (a)v(a) du dx
=
(b)v(b) (since v (a) = 0) −P (b) du du
a
β 1 = P (b)u(b)v(b) β 2
du since (b) = dx
−
β 1 u(b) . β 2
Therefore, b β 1 du dv (Lu,v)w = P (b)u(b)v(b) + P (x) (x) (x) dx + β 2 dx dx a The expression on the right is symmetric in u and v, which shows that
b
R(x)u(x)v(x) dx.
a
(u, L(v))w = (L(v), u)w = (L(u), v)w , verifying 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 (a)(−1) = κ (a). −κ du 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
−
κ
(a) = −γT , −γu(a) + κ du dx
83
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.
(a) The eigenpairs are
nπ ln 2
2
, un (x) = sin
nπ λn = 4 + β + ln 2
2
nπ ln x , n = 1, 2, 3, . . . . ln 2
2
, un (x) = x sin
(b) 0 is an eigenvalue if and only if β =
−4 −
for some positive integer n.
nπ ln 2
nπ ln x , n = 1, 2, 3, . . . . ln 2
2
(c) There is exactly one negative eigenvalue if and only if 2
2
4π π −4 − (ln2) ≤ β < −4 − (ln2) 2
2
,
while there are exactly k negative eigenvalues if and only if 2 2
2 2
1) π k π −4 − (k + ≤ β < −4 − (ln2) (ln2) 2
2
.
9. The eigenpairs are nπ λn = 2 + ln 2
10.3
2
nπ nπ ln x nπ ln x , un (x) = x −1 cos + sin ln 2
ln 2
ln 2
, n = 1, 2, 3, . . . .
Numerical Methods for Sturm-Liouville problems
1. We wish to apply the finite element method to
−
d dx
du P (x) dx
+ R(x)u = λw(x)u, a < x < b, 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
{
}
×
b
Aij =
a b
W ij =
dφj dφi (x) (x) + R(x)φj (x)φi (x) dx, P (x) dx dx
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.)
85
10.4. EXAMPLES OF STURM-LIOUVILLE PROBLEMS
{
}
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 A ij , W ij . 5. Let λ C and x We then obtain
∈
∈ C
n
×
satisfy Ax = λWx, and let (x, y)W = (Wx) y. We can assume x satisfies (x, x)W = 1.
·
·
λ = λ(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
n
u(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 1
0
1 u(x) dx =
Since the graph of each φ i , i = 1, . . . , n
−
{
n
1 αi φi (x) dx = i=0
0
n
}
αi
φi (x) dx.
0
i=0
1, forms a triangle of height 1 and base 2h (h = /n), it follows that
φi (x) dx = h, i = 1, 2, . . . , n
0
− 1;
similarly,
φ0 (x) dx =
0
φn (x) dx =
0
Thus the average value of u on [0, ] is h
1 α0 + α1 + 2
h , i = 1, 2, . . . , n 2
··· + α − + 21 α n 1
n
− 1.
.
3. The fundamental frequency is about 3254 Hz.
10.5
Robin boundary conditions
1. The solution to the IBVP is u(x, t) =
∞
an (t)cos
n=1
where 2
an (t) = c n e−κ(2n−1)
π2 t/(42 ρc)
(2n
, cn =
− 1)πx 2
,
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: 2
− ddxu = λu, 2
0 < x < ,
u(0) = 0, du () + αu() = 0 dx where κ > 0, α < 0. We can write the second boundary condition as κ
du () dx
− αu() = 0,
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, α . λ. There are solutions s k = (2k + 1)π/2, k = 1, 2, 3, . . .. The corresponding eigenpairs are tan(s) =
where s =
√
s2k sk x , ψk (x) = sin , k = 1, 2, 3, . . . . 2 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 α = 1, then 0 is not an eigenvalue. Case 3 (λ < 0) Negative eigenvalues must satisfy the equation
λk =
∈
√ −λ √ tanh( −λ) = α
or
s , s > 0, α 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/). tanh(s) =
√ −
≤
−
Therefore, when α < 0, there is a sequence of eigenpairs λk =
s2k sk x , ψk (x) = sin , k = 1, 2, 3, . . . , 2
where
. (2k + 1) 2 π 2 λk = . 42 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) dx
+ R(x)u = F (x), a < x < b, u(a) = 0,
du β (b) + αu(b) = 0, dx
87
10.6. FINITE ELEMENT METHODS FOR ROBIN BOUNDARY CONDITIONS
∈
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
du P (x) dx
+ R(x)u = F (x), a < x < b,
du (a) = 0, dx du β 1 u(b) + β 2 (b) = 0 dx
−α u(a) + α 1
is du P (a) (a)v(a) dx
−
du P (b) (b)v(b) + dx
2
b
du dv P (x) (x) (x) dx + dx dx
a
for all test functions v. The boundary conditions imply du α1 du (a) = (b) = u(a), dx α2 dx
b
b
R(x)u(x)v(x) dx =
a
F (x)v(x) dx
a
− β β u(b), 1 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
b
a
du dv P (x) (x) (x) dx + dx dx
b
b
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 P (a), i = j = 0, α2 β1 Gij = P (b)β, i = j = n, β2 0, otherwise.
{
}
∈
5. Suppose u satisfies
= F (x), 0 < x < , −κ du dx du (0) = 0, dx κ
du () + αu() = 0. dx
Then
0
F (x) dx =
− κ
0
d2 u (x) dx = dx2 =
− κ
du () dx
du κ () dx
−
−
du (0) dx
du since (0) = 0 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 un U = 1 for all n and Lun V as n . (Here we are using the alternate definition of the norm of L: L = sup Lu V : u U, u U = 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 Lun V implies that Lunk V , 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 −
d dx
du P (x) dx
+ R(x)u = λw(x)u, a < x < b,
du (a) = 0, dx du β 1 u(b) + β 2 (b) = 0, dx
α1 u(a) + α2
where P and w are assumed to be positive on the interval [a, b]. Let W (x) =
φ(x) dφ (x) dx
ψ(x) dψ (x) dx
= φ(x)
dψ (x) dx
(x)ψ(x) − dφ dx
be the Wronskian of φ, ψ. Since φ, ψ satisfy the linear homogeneous ODE
−
d dx
du P (x) dx
+ (R(x)
− λw(x)u = 0
on [a, b], we can show that φ, ψ is linearly independent by proving that W (x) = 0 for any given x have dφ dφ α1 α1 φ(a) + α2 (a) = 0 (a) = φ(a), dx dx α2 and similarly dψ α1 (a) = ψ(a). dx α2 Therefore, dψ dφ α1 α1 W (a) = φ(a) (a) (a)ψ(a) = φ(a)ψ(a) + φ(a)ψ(a) = 0. dx dx α2 α2 This shows that φ, ψ is linearly dependent, as desired.
{
}
⇒
−
−
{
}
−
−
∈ [a, b]. We
Chapter 11
Problems in Multiple Spatial Dimensions 11.1
Physical models in two or three spatial dimensions
1. Let Ω
2
⊂R
be the rectangular domain
Ω= x and let F : R 2
2
→R
2
: a < x1 < b, c < x2 < d ,
∈R
be a smooth vector field. Then b
d
∇ · { − } { − − · · · · F =
Ω
a b
c
a d
c
c
a
∂F 1 ∂ F 2 (x1 , x2 ) + (x1 , x2 ) dx 2 dx1 ∂x 1 ∂x 2
d
∂F 1 (x1 , x2 ) dx2 dx1 + ∂x 1
b
∂F 1 (x1 , x2 ) dx1 dx2 + ∂x 1
= =
b
a
∂F 2 (x1 , x2 ) dx2 dx1 ∂x 2
d
∂F 2 (x1 , x2 ) dx2 dx1 ∂x 2
c
b
a
d
=
d
c
b
F 1 (b, x2 )
F 1 (a, x2 ) dx 2 +
c
F 2 (x1 , d)
}
F 2 (x2 , c) dx 1
a
d
=
d
F 1 (b, x2 ) dx2
F 1 (a, x2 ) dx2 +
c
=
b
c
F n+
Γ2
b
F 2 (x1 , d) dx1
a
F n+
Γ4
F n+
Γ3
F n =
Γ1
−
F 2 (x2 , c) dx1
a
F n,
·
∂ Ω
as desired. Here Γ1 , Γ2 , Γ3 , Γ4 denote the four sides of ∂ Ω, labeled as in Figure 11.3 on page 452 of the text.
∇ · F(x) = 1, and therefore
3. We have
∇ ·
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
·
∂ Ω
2π
F · n =
·
sin(2θ) + sin2 (θ) dθ = π.
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,
− ∇ · (k∇u)v = Ω
∇ · ∇ − ∇ · k u
Ω
89
v
Ω
(k( u)v).
∇
90
CHAPTER 11. PROBLEMS IN MULTIPLE SPATIAL DIMENSIONS By the divergence theorem,
∇ ·
(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, (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
∇ · F(u(x)) = ∂x∂ (F (u(x))) + ∂x∂ (F (u(x))) + ∂x∂ (F (u(x))) 1
1
2
2
3
3
∂u ∂ u ∂u (x) + F 2 (u(x)) (x) + F 3 (u(x)) (x) ∂x 1 ∂x 2 ∂x 3 = F (u(x)) u(x). = F 1 (u(x))
·∇
9. Using the product rule for scalar functions, we obtain 3
∇ · (v∇u) = 11. Suppose u, v
2 m (Ω).
∈ C
3
i=1
∂ ∂ u = v ∂x i ∂x i
i=1
∂ v ∂u ∂ 2 u +v 2 ∂x i ∂x i ∂x i
3
=
i=1
=
3
∂v ∂u ∂ 2 u +v ∂x i ∂x i ∂x 2i i=1
∇v · ∇u + v∆u.
Then
(Lmu, v) =
−
v∆u
∇ · ∇ − ∇ · ∇ − −
=
u
Ω
v
Ω
=
v
∂ Ω
u
∂u (Green’s first identity) ∂ n
v
Ω
=
u
∂ Ω
=
∂v ∂ n
u∆v (Green’s first identity)
Ω
u∆v = (u, Lmv).
Ω
The boundary terms vanish because the product that forms the integrand is zero over the entire boundary. For example, ∂u =0 v ∂ Ω ∂ n
since v = 0 on Γ1 and ∂u/∂ n = 0 on Γ2 .
11.2 1.
Fourier series on a rectangular domain (a) m
cmn =
− 4 (5 + 7(−1)m n) (π−1 + ( −1)
n
)
3 3 6
(b) The graphs of the error are given in Figure 11.1. 3. The solution is given by u(x) =
∞ ∞
cmn sin(mπx1 )sin(nπx2 ), λmn m=1 n=1
where the c mn are the coefficients from Exercise 1 and λmn = (m2 + n2 )π 2 .
91
11.2. FOURIER SERIES ON A RECTANGULAR DOMAIN x 10
−3
2 0 −2 1 0.5 0
0
0.2
x
0.8
1
x
2
x 10
0.6
0.4 1
−3
2 0 −2 1 0.5 0
0
0.2
x
0.6
0.4
0.8
1
x
2
1
Figure 11.1: The error in approximating f (x, y) by the first 4 terms (top) and the first 25 terms (bottom) of the double Fourier sine series. (See Exercise 11.2.1.) 5.
(a) The IBVP is ρc
∈
The domain Ω is the rectangle x (b) The solution is
∂u ∂t
− κ∆u
=
0.02, x
u(x, 0)
=
5, x
u(x, t)
=
∈ Ω, t > 0,
∈ Ω, 0, x ∈ ∂ Ω, t > 0.
R2 : 0 < x1 < 50, 0 < x2 < 50 .
∞ ∞
u(x, t) =
amn (t)sin
m=1 n=1
where amn (t)
=
bmn
=
cmn
=
λmn
=
mπx1 nπx2 sin , 50 50
cmn cmn bmn 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
−
− − − −
− − − −
(c) The steady-state temperature u s (x) satisfies the BVP
−κ∆u u(x)
The solution is
∞ ∞
=
0.02, x
=
0, x
∈ Ω,
∈ ∂ Ω.
cmn mπx1 nπx2 sin sin us (x) = . 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 b oundary, and the wave travels at a speed of 261 2 . 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 x
20 0
0
2
x
1
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) =
bmn sin (mπx1 )sin(nπx2 ),
m=1 n=1
then we find that
0,
bmn = We then have that u(x, t) =
sin2 (mπ/2) , 1250m2 π2
∞ ∞
m = n, m = n.
amn (t)sin(mπx1 )sin(nπx2 ),
m=1 n=1
where a mn (t) satisfies
d2 amn + λmnc2 amn 2 dt amn(0) damn (0) dt The result is amn (t) = Thus the solution is u(x, t) =
∞
=
0,
=
bmn ,
=
0.
0, bmm cos (c λmmt),
√
m = n, m = n.
√
bmm cos (c λmmt)sin(mπx1 )sin(mπx2 ).
m=1
Four snapshots of u are shown in Figure 11.3. 13. The solution is u(x) =
∞ ∞
4sin(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 ∂t
− κ∆u
=
0.02, x
u(x, 0)
=
u(x, t) ∂u (x, t) ∂ n
=
∈ Ω, 0, x ∈ Γ ∪ Γ , t > 0, 0, x ∈ Γ ∪ Γ , t > 0.
=
∈ Ω, t > 0,
5, x
1
4
2
3
The domain Ω is the rectangle
∈ x
R2 : 0 < x1 < 50, 0 < x2 < 50 .
93
11.2. FOURIER SERIES ON A RECTANGULAR DOMAIN −4
−4
x 10
x 10
1
1
0
0
−1 1
−1 1
1
0.5 0 0
2
1
0.5
0.5
x
0.5
x
x
2
1
−4
0 0
x
0 0
x
1
−4
x 10
x 10
1
1
0
0
−1 1
−1 1
1
0.5 0 0
2
1
0.5
0.5
x
0.5
x
x
2
1
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
x
0
2
x
1
Figure 11.4: The solution to Exercise 11.2.13. (b) The solution is u(x, t) =
∞ ∞
amn (t)sin
m=1 n=1
where
amn(t)
=
bmn
=
cmn
=
bmn
(2m
=
1
100
− κλc
mn mn
sin
(2n
e−κλmn t/(ρc) +
− 1)πx
2
100
cmn , κλmn
80 (2m 1)(2n 1)π2 8 25(2m 1)(2n 1)π 2
−
−
− ((2m − 1)
−
2
λmn
− 1)πx
+ (2n 10000
2
− 1) )π
2
.
(c) The steady-state temperature u s (x) satisfies the BVP
−κ∆u u(x) ∂u (x) ∂ n
The solution is
∞ ∞
cmn sin us (x) = κλmn m=1 n=1
= = =
∈ Ω, 0, x ∈ Γ ∪ Γ , 0, x ∈ Γ ∪ Γ .
0.02, x
(2m
1
4
2
3
− 1)πx 100
1
sin
(2n
− 1)πx 100
2
.
,
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) u s (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 x
20 0
0
2
x
1
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, c
50
. . = 0.0116362, c60 =
3. The solution is u(r, θ) =
∞
−0.00722147.
am0 J 0 (s0m r),
m=1
where
am0 = the c m0 are the coefficients of f (r, θ) = 1 that
− r, and the s a10 a20 a30 a40 a50 a60
. = . = . = . = . = . =
cm0 , s20m
0m are
the positive roots of J 0 . A direct calculation shows
0.226324,
−0.0157747, 0.00386078,
−0.00148156, 0.000716858,
−0.000399554.
The solution (approximated by six terms of the series) is graphed in Figure 11.6. 5. The solution is u(r, θ) =
∞
am0 J 0 (s0m r),
m=1
where
am0 =
cm0 , s20m
95
11.3. FOURIER SERIES ON A DISK
0.25 0.2 0.15 0.1 0.05 0 1 0.5
1 0.5
0
0
−0.5 x
−0.5 −1
2
−1
x
1
Figure 11.6: The approximation to solution u, computed with six terms 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 J 0 . 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 terms of the series) is graphed in Figure 11.7.
0.15
0.1
0.05
0 1 0.5
1 0.5
0
0
−0.5 x
−0.5 −1
2
−1
x
1
Figure 11.7: The approximation to solution u, computed with six terms of the (generalized) Fourier series (see Exercise 11.3.5). 7. If we write φ(r, θ) =
∞
cm0 J 0 (αm0 r) +
m=1
∞
(cmn cos (nθ) + dmn sin (nθ)) J n (αmn r),
m=1 n=1
where the c mn and d mn are known, and u(r,θ,t) =
∞ ∞
am0 (t)J 0 (αm0 r) +
m=1
∞ ∞
(amn(t)cos(nθ) + bmn (t)sin(nθ)) J n (α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 c m1 , m = 1, 2, 3, . . .. Therefore,
∞
u(r,θ,t) =
am1 (t)cos(θ)J 1 (αm1 r),
m=1
with a m1 (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 5
10 5
0 0
−5
−5 −10
x
2
−10
x
1
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 c
π2 ( πA)2
√
2 + (√ ππA)2
2π
=
√ πA. The fundamental
c 1 . c = 0.398942 . A 2π A
√
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 .
K=
while the load vector is
−−
4.4815 1.1759 1.1759 0.0000
−1.1759 −1.1759 4.9259 0.00000 1.3426
− .
F=
0.0000 1.3426 1.3426 5.8148
− −
0.0000 4.9259 1.3426
−
0.11111 0.11111 0.11111 0.11111
.
The resulting weights for the finite element approximation are given by
u = K
−1 F =.
0.048398 0.044979 0.044979 0.039879
.
,
97
11.4. FINITE ELEMENTS IN TWO DIMENSIONS 1
0.8 3
4
1
2
0.6 2
x
0.4
0.2
0 0
0.2
0.4
0.6 x
0.8
1
1
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
0.6 2
x
0.4
0.2
0 0
1
2
0.2
3
0.4
0.6
4
0.8
1
x
1
Figure 11.10: The mesh for Exercise 11.4.3.
5. Let V n be the space of continuous piecewise linear functions relative to a given triangulation of the domain Ω, and let φ1 , . . . , φn be the standard basis for V n . Suppose that u : Ω R belongs to L 2 (Ω). We wish to find v V n that minimizes w u L2 (Ω) over all w V n . We can write v = n i=1 αi φi and use the projection theorem: v must satisfy (u v, w)L2 (Ω) = 0 for all w V n , 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
αj φj , φi = (u, φi ), i = 1, 2, . . . , n
j=1
n
α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, M ij =
φj φi , i , j = 1, 2, . . . , n ,
Ω
and f is the vector defined by f i =
vφi , i = 1, 2, . . . , n .
Ω
7.
(a) A direct calculation shows that the inhomogeneous Dirichlet problem
−∇ · (k(x)∇u)
=
u(x)
=
f (x), x
∈ Ω, g(x), x ∈ ∂ Ω,
has weak form
∈ G + V,
2 D
∈ V = C (Ω), 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. u
a(u, v) = (f, v) for all v
In the Galerkin method, we replace V by a finite-dimensional subspace V n and solve
∈ V , a(w, v) = (f, v) − a(G, v) for all v ∈ V .
w
n
n
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 f i = (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 f i = The load vector is .
f =
while the solution to Ku = f is .
u=
−a(G, φ ), i = 1, 2, 3, 4.
i
0.22222222222222 1.55555555555556 0.22222222222222 1.55555555555556 0.27777777777778 0.61111111111111 0.27777777777778 0.61111111111111
,
.
99
11.4. FINITE ELEMENTS IN TWO DIMENSIONS 9.
(a) Consider the BVP
−∇ · (k(x)∇u) = f (x), in Ω, ∂u (x) = h(x), x ∈ ∂ Ω. ∂ n ˜ = C (Ω) and applying Green’s first identity yields Multiplying the PDE by a test function v ∈ V k∇u · ∇v = f v + kvh for all v ∈ V˜ .
(11.1)
2
Ω
Ω
∂ Ω
˜m equal to the space of all continuous piecewise linear function We will apply Galerkin’s method with V ˜m as φ1 , φ2 , . . . , φn , and the nodes as relative to a given triangulation . We label the standard basis of V ˜ z1 , z2 , . . . , zn . (Every node in the mesh is now free.) Thus V n = span φ1 , φ2 , . . . , φn . We then define the Galerkin approximation w n by
T
{
∈ V ˜ , a(w
wn Writing w n =
m j=1
n , φi ) = (f, φi ) +
n
}
khφi , i = 1, 2, . . . , m .
∂ Ω
β j φj and
w =
β 1 β 2 .. . β n
˜ = F ˜ . The stiffness matrix K ˜ we obtain the system Kw
,
m m
∈R
× is given by
˜ ij = a(φj , φi ), i , j = 1, 2, . . . , n , K ˜ and the load vector F
m
∈R
by ˜i = (f, φi ) + F
khφi , i = 1, 2, . . . , n .
∂ Ω
˜i is zero unless zi is a boundary node. The boundary integral in the expression for F (b) The compatibility condition for (11.1) is determined by the following calculation:
Ω
f =
− ∇· Ω
(k u) =
∇
∂u = k ∂ Ω ∂ n
−
−
kh.
∂ Ω
(c) Now consider (11.1) with 3 3x21 , h(x) = , 4 5 where Ω is the unit square. We will produce the finite element solution using the regular grid with 18 ˜ triangles (16 nodes). The stiffness matrix K R16×16 is singular, as should be expected, and there are infinitely many solutions. The unique solution u with last component equal to zero is k(x) = 1, f (x) = x 1 x2 +
−
∈
u =
0.3622 0.2984 0.1344 0.0843 0.3728 0.3302 0.2081 0.0255 0.3803 0.3447 0.2334 0.0591 0.3665 0.3257 0.1878 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(α) = −x sin (α) + x cos (α), 1
2
as desired. 3. We have 2
2
d ψ dψ = φ(r) ⇒ −r − ddrψ − 1r dψ − dr = rφ(r) dr dr (r) = rφ(r). ⇒ − drd r dψ dr 2
2
This last ODE can be solved by integrating twice to yield the general solution r
ψ(r) = c1 + c2 ln (r)
−
r φ(s)sds. s
ln
0
5. Let V be the space of all smooth functions defined on R2 and having compact support. We wish to show that the variational problem
∇
u(x)
R2
is equivalent to
2π
0
in polar coordinates.
· ∇v(x) dx = v(y) for all v ∈ V
∞ ∂u ∂v
1 ∂u ∂v + 2 ∂r ∂r r ∂θ ∂θ
0
v rdrdθ = v(y) for all v
∈ V
(a) We first convert
∇
u(x)
R2
· ∇v(x) dx
to p olar coordinates, where y R2 is fixed and represents the origin for the polar coordinates. By the chain rule, we have ∂u ∂u sin(θ) ∂u ∂u ∂u cos(θ) ∂u = cos (θ) , = sin (θ) + ∂x 1 ∂r r ∂θ ∂x 2 ∂r r ∂θ (cf. the derivation in Section 11.3.1) and similarly for v. Therefore,
∈
−
∇u · ∇v =
sin (θ) ∂u r ∂θ ∂u cos(θ) ∂u sin(θ) + ∂r r ∂θ ∂u cos(θ) ∂r
−
sin (θ) ∂v r ∂θ ∂v cos(θ) ∂v sin(θ) + ∂r r ∂θ ∂v cos(θ) ∂r
−
+
∂u ∂v sin(θ)cos(θ) ∂u ∂v sin(θ)cos(θ) ∂u ∂v sin2 (θ) ∂u ∂v + + ∂r ∂r r ∂θ ∂r r ∂r ∂θ r2 ∂θ ∂θ cos2 (θ) ∂u ∂v ∂u ∂v sin (θ)cos(θ) ∂u ∂v sin (θ)cos(θ) ∂u ∂v sin2 (θ) + + + ∂r ∂r r ∂θ ∂r r ∂r ∂θ r2 ∂θ ∂θ 2 2 ∂u ∂v sin (θ) cos (θ) ∂ u ∂v = (sin2 (θ) + cos2 (θ)) + + ∂r ∂r r2 r2 ∂θ ∂θ 1 ∂u ∂v ∂u ∂v = + 2 . ∂r ∂r r ∂θ ∂θ = cos2 (θ)
−
−
101
11.5. THE FREE-SPACE GREEN’S FUNCTION FOR THE LAPLACIAN It follows that
2π
∇
u(x)
R2
· ∇v(x) dx =
∞ ∂u ∂v
0
0
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 by a test function, and integrating. The right-hand side becomes
−
δ (x
R2
− y) in polar coordinates, multiplying
− y)v(x) dx = v(y).
On the left side, we have 2π
0
=
∞
0 2π
0 2π
=
∂ 2 u ∂r 2
∂r 2
∂u r v ∂r
0
0
∞
+
r2
∂u v ∂θ
0
2π
0
1 ∂u ∂v ∂u ∂v + 2 ∂r ∂r r ∂θ ∂θ
0
2π
0
∂u ∂v dθ ∂θ ∂θ
∞ 1 ∂u ∂v
0
∂r
0
2π
r2 ∂θ ∂θ
r2
0
∞ ∂u
2π
∞ 1
vdrdθ
r dr +
∂r ∂r
0
0
∂r
0
∞ ∂u ∂v
vrdrdθ
∞ ∂u
0
∂u ∂v rdrdθ + ∂r ∂r
2π
=
2π
vrdr dθ
∞ 1
2π
1 ∂u r ∂r
∞ ∂ 2 u
0
0
=
1 ∂ 2 u r2 ∂θ 2
− − − − − − − − −
∂ 2 u v dθ r dr ∂θ 2
0 2π
∞ ∂u
v dr dθ
0
0
∂r
vdrdθ
−
r dr
rdrdθ
rdrdθ.
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 ∂g (ρ) = ∂ρ
− ρc
1 2
and therefore 2π
π
∞ ∂g ∂v
0
0
0
∂ρ ∂ρ
2π
2
ρ sin(φ) dρdφdθ =
−c
1
0
π
0
2π
=
−c
1
0
−c
1
0
−c
1
0
= c 1 v(y)
∞ ∂v
0
∂ρ
0
ρ2 sin(φ) dρdφdθ
sin (φ) dρdφdθ
π
∞ ∂v
sin(φ)
0
2π
=
ρ2 ∂ρ
0
π
2π
=
∞ 1 ∂v
− 0
∂ρ
π
dρ dφdθ
sin(φ) ( v(y)) dφdθ
0 2π
π
sin(φ) dφdθ
0
0
2π
= 2c1 v(y)
dθ
0
= 4πc1 v(y).
From this we see that g is a solution of (11.83) if and only if c 1 = 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
−
(∆u)v =
Ω
⇒−
fv
Ω
∂u v + ∂n
∂ Ω
The b oundary condition yields
∂u α = − ∂n α
1
∇ · ∇ u
v =
Ω
f v.
Ω
u on ∂ Ω,
2
and hence we obtain
α1 α2
as desired.
∇ · ∇
uv +
u
∂ Ω
v =
Ω
f v,
Ω
3. Let G be the free-space Green’s function for the Laplacian and define ∂ G (x; y) ∂ nx
w(x; y) = α1 G(x; y) + α 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 u2 (x) =
∂G (x; y). ∂ nx
vΩ (x; y)f (y) dy.
Ω
Then, for all v
2
∈ C (Ω), − (∆
⇒− ⇒ ⇒
v v = 0
x Ω
Ω
∂ Ω
α1 α2 α1 α2
∇ − ∇
∂v Ω (x; y)v(x) dσx + ∂ nx vΩ (x; y)v(x) dσx
∂ Ω
1 α2
vΩ (x; y)v(x) dσx +
∂ Ω
v (x; y)
x Ω
Ω
· ∇v(x) dx = 0
w(x; y)v(x) dσx +
∂ Ω
∇
v (x; y)
x Ω
Ω
·∇
v (x; y)
x Ω
Ω
1 v(x) dx = α2
· ∇v(x) dx = 0
w(x; y)v(x) dσx .
∂ Ω
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 α1 α2
Ω
∂ Ω
vΩ (x; y)v(x) dσx f (y) dy+
∇
v (x; y)
x Ω
Ω
=
1 α2
Ω
Ω
· ∇v(x) dx
f (y) dy
w(x; y)v(x) dσx f (y) dy.
∂ Ω
Changing the order of integration in the two integrals on the left yields α1 α2
∂ Ω
Ω
vΩ (x; y)f (y) dy v(x) dσx +
∇
v (x; y)f (y) dy
x Ω
Ω
=
1 α2
Ω
Ω
∂ Ω
·∇
v(x) dx
w(x; y)v(x) dσx f (y) dy.
11.6. THE GREEN’S FUNCTION FOR THE LAPLACIAN ON A BOUNDED DOMAIN Since
∇
∇
v (x; y)f (y) dy =
x Ω
Ω
we obtain
α1 α2
vΩ (x; y)f (y) dy =
Ω
u2 (x)v(x) dσx +
∂ Ω
∇
u2 (x)
Ω
Ω
∇u (x), 2
· ∇v(x) dx
1 α2
=
103
w(x; y)v(x) dσx f (y) dy,
∂ Ω
as desired. 5. Let G be the free-space Green’s function for the Laplacian, and suppose v Ω (x; y) satisfies the BVP
−∆ v
x Ω =
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 u
∈
2 C D (Ω),
∇
u(x)
Ω
· ∇v(x) dx =
f (x)v(x) dx for all v
Ω
2 D (Ω).
∈ C
We prove that G Ω is the desired Green’s function by proving that u(x) =
GΩ (x; y)f (y) dy
Ω
satisfies the ab ove variational problem. We have u(x) =
G(x; y)f (y) dy
Ω
where u1 (x) =
−
vΩ (x; y)f (y) dy = u 1 (x)
Ω
G(x; y)f (y) dy, u2 (x) =
Ω
− u (x), 2
vΩ (x; y)f (y) dy.
Ω
Letting V be the space of all smooth functions with compact support, we have
∇
u1 (x)
R2
· ∇v(x) dx =
f (x)v(x) dx for all v
R2
∈ V,
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 C D (Ω), we can extend v to be identically zero outside of Ω and thereby obtain an element v of V . We then have
∈
∇
u1 (x)
Ω
· ∇v(x) dx =
f (x)v(x) dx for all v
Ω
On the other hand, we have
−∆u (x) = 0 for all x ∈ Ω; 2
multiplying by v
∈
2 C D (Ω)
and integrating yields
−
∆u2 (x)v(x) dx = 0 for all v
Ω
2 D
∈ C (Ω).
Integrating by parts on the left, we obtain
∇ Ω
u2 (x)
2 D (Ω).
· ∇v(x) dx = 0 for all v ∈ C
2 D
∈ C (Ω).
104
CHAPTER 11. PROBLEMS PROBLEMS IN MULTIPLE MULTIPLE SPA SPATIAL DIMENSIONS DIMENSIONS (Notice that the boundary term vanishes because v = 0 on ∂ Ω.) Ω.) Putting Putting together the variational variational equations equations satisfies by u by u 1 and u and u 2 , we see that u that u = u = u 1 u2 satisfies
−
∇
u(x)
Ω
Also, for x
· ∇v(x) dx =
all v f ( f (x)v (x) dx for all v
Ω
2 D
(Ω). ∈ C (Ω).
Ω, we have ∈ ∂ Ω,
u(x) =
G(x; y)f ( f (y) dy
Ω
This shows that u that u
−
vΩ (x; y)f ( f (y) dy =
Ω
G(x; y)f ( f (y) dy
Ω
2 hence u satisfies satisfies D (Ω); hence u
∈ C
−
G(x; y)f ( f (y) dy = 0.
Ω
the given BVP, which shows that G that GΩ is the desired Green’s function.
7. Green’s Green’s second identity identity is
( v ∆u
− u∆v) =
Ω
Let G Let G Ω be the Green’s function for the BVP
∂u v ∂ n
∂ Ω
−
∂v u ∂ n
.
= f ((x) in Ω −∆u = f α1 u + α2 and, for y, z
define u((x) = G ∈ Ω, define u
), v((x) Ω (x; y), v
∂u = 0 on ∂ on ∂ Ω Ω, ∂n
= GΩ (x; z). Since G Since G Ω (x; y) satisfies
−∆
= δ (x y) in Ω GΩ = δ ∂G Ω = 0, on ∂ on ∂ Ω Ω, α1 GΩ + α2 ∂ nx
−
x
we see that
v ∆u =
−v(y) = −G(y; z),
u∆v =
−u(z) = −G(z; y).
Ω
Ω
Therefore,
(v ∆u
Ω
We also have
G (z; y) − G(y; z). − u∆v) = G( ∂u = ∂ n
or
− αα
1
u on ∂ Ω ∂ Ω,
2
∂G Ω (x; y) = ∂ nx
− αα
GΩ (x; y).
∂G Ω (x; z) = ∂ nx
− αα
GΩ (x; z).
1 2
Similarly, 1 2
Therefore,
∂ Ω
∂u v ∂ n
−
−
∂v = u ∂ n
∂ Ω
α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) that is, G is, G((z; y) = G( G(y; z), as desired. 9. Let Let Ω
2
⊂R
− G(y; z) = 0,
be given, and let u let u be the solution of Ω, −∆u = 0 in Ω,
u = φ = φ((x) on ∂ on ∂ Ω Ω.
105
11.6. 11.6. THE GREEN’S GREEN’S FUNCTION FUNCTION FOR THE LAPLACIAN LAPLACIAN ON A BOUNDED BOUNDED DOMAIN DOMAIN Let v Let v((x) = GΩ (x; y), where G where G Ω is the Green’s function for the BVP = f ((x) in Ω, Ω, −∆u = f u = 0 on ∂ on ∂ Ω Ω.
{ ∈
−
}
∈
⊂
Finally, define Ω = x Ω : x y > , where y Ω is given and > 0 is small enough that B (y) Ω. Notice that ∆u ∆u = 0 in Ω, ∆v ∆v = 0 in Ω , and v and v((x) = 0 for all x ∂ Ω. Ω. Applying Green’s second identity to u to u and v yields
− u∆v) =
(v ∆u
Ω
⇒ 0 = ⇒ 0 =
∂ Ω
∂ Ω
∂u v(x) (x) ∂ n
∂ Ω
∂u 0 (x) ∂ n
−
∂u v ∂ n
−
Now, recall that G that G Ω = G = G
−v
Ω,
∂v u dσ x ∂ n
∂v u(x) (x) dσ x + ∂ n
S (y)
∂G Ω φ(x) (x; y) dσ x + ∂ nx
∂G Ω (x; y) dσx = φ(x) ∂ nx ∂ Ω
⇒
−
∈
∂u v (x) (x) ∂ n
∂u GΩ (x; y) (x) ∂ n
S (y)
∂u GΩ (x; y) (x) dσx ∂ n S (y)
−
−
u(x)
S (y)
∂v u(x) (x) dσ x ∂ n
−
∂G Ω u(x) (x; y) dσ x ∂ nx
∂G Ω (x; y) dσx . ∂ nx
where G where G is the free-space Green’s Green’s function function for the Laplacian and v and v Ω (x; y) satisfies satisfies
−∆
Ω, v = 0 in Ω,
x Ω
vΩ (x; y) = G( G(x; y) on ∂ on ∂ Ω Ω. We therefore have ∂GΩ (x; y) dσx = φ(x) ∂ nx ∂ Ω
∂u G(x; y) (x) dσx ∂ n S (y)
vΩ (x; y)
S (y)
−
u(x)
S (y)
∂u (x) dσx + ∂ n
S (y)
∂ G (x; y) dσx ∂ nx
u(x)
∂v Ω (x; y) dσx . ∂ nx
The integrands for the last two integrals are continuous, with no singularities as integrals tend to zero as as 0+ . We have
+
→ 0
→
G(x; y) = and therefore G therefore G((x; y) is constant for x
−
. Therefore, Therefore, these two two
− 21π ln (x − y) ,
∈ S (y). It follows that
G(x; y)
S (y)
∂u (x) dσx ∂ n
is just a multiple of
S (y)
∂u (x) dσx = ∂ n
∆u(x) dx = 0
B (y)
(applying the divergence theorem and using the fact that ∆u ∆u = 0 in Ω). It remains only to calculate
u(x)
S (y)
∂ G (x; y) dσx . ∂ nx
Now,
∇
x
G(x; y) =
and, on S on S (y), nx =
− 2πxx−−yy
2
− xx −− yy
(notice that this is the outward-pointing normal to Ω ). Therefore, ∂G Ω (x; y) = ∂ nx
∇
x
·
G(x; y) nx =
1 . 2π
106
CHAPTER 11. PROBLEMS PROBLEMS IN MULTIPLE MULTIPLE SPA SPATIAL DIMENSIONS DIMENSIONS Thus
1 ∂ G u(x) (x; y) dσx = ∂ nx 2π S (y)
which is the average value of u u on S (y), and therefore
u(x)
S (y)
as as
+
→0
u(x) dσx ,
S (y)
∂ G (x; y) dσx ∂ nx
→ u(y) +
→ 0 in Green’s second identity, that = −u(y),
. Therefore, it finally follows, by taking the limit as
φ(x)
∂ Ω
∂G Ω (x; y) dσx ∂ nx
which, given the symmetry of G G Ω , yields the desired result. 11. In the text it is shown shown that, that, if u u solves in B −∆u = 0 in B
(0), R (0),
u = φ = φ((x) on S on S R (0), (0),
then u(x) =
SR (0)
R2 x 2 φ(y) dσy . 2πR x y 2
− −
Now suppose that u that u is is to represented as a function of polar coordinates (r, (r, θ). Let the circle S circle S R (0) b e represent represented ed by y = (R cos(ψ cos(ψ ), R sin(ψ sin(ψ)), 0 ψ 2π . Then
≤ ≤
dσy =
∂y 1 ∂ψ
2
+
2
∂y 2 ∂ψ
dψ = dψ =
R2 cos2 (ψ) + R2 sin2 (ψ ) dψ = dψ = Rdψ. Rdψ.
With x represented by the polar coordinates (r, (r, θ), we have x 2
2
cos(ψ)) − R cos(ψ
2
2
2
+ R2
cos(ψ − θ) − 2rR cos(ψ (where we used the trignometric trignometric identity identity cos (θ)cos(ψ )cos(ψ) + sin sin (θ)sin(ψ )sin(ψ) = cos cos (ψ − θ)). Changing Changing variables variables from
x − y
= (r cos(θ cos(θ)
2
= r and + (r (r sin(θ sin(θ) − R sin(ψ sin(ψ )) = r
y to (r, θ) thus yields
SR (0)
x 2 1 R2 φ(y) dσy = 2πR x y 2 2π
− −
2π
0
R2 + r 2
R2 r 2 2Rr cos(ψ cos(ψ
−
−
− θ) φ(ψ) dψ,
as desired. 13. In Chapter 9, it was shown shown that the solution to
−
d dx
du P ( P (x) dx
+ R(x)u = 0, a < x < b, u(a) = u a , u(b) = u b
is
∂G ∂G (y ( y ; a) P ( P (b)ub (y ( y ; b), ∂y ∂y where G where G is the Green’s function for the inhomogeneous differential equation with homogeneous boundary conditions. If we take P take P ((x) = 1, we obtain u(y) = P ( P (a)ua
u(y ) =
−
∂G − − ∂G (y ( y ; a)u + (y ( y ; b)u ∂y ∂y a
b
.
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, ( a, b) (that is, the set a, b ). Notice Notice that the unit unit normal normal to Ω at x at x = b = b is n = 1, while at x at x = a = a it is n is n = 1, so that
{ }
−
(y ( y; a) − ∂G ∂y
107
11.7. GREEN’S FUNCTION FOR THE WAVE EQUATION 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 ∈ B (z), u(x) = 0, x ∈ S (z) R
R
is u(x) = 2
Substituting x = z , we see that z
− y
u(z) =
= R
x z 2 R2 φ(y) dσy . 2πR x y 2
− − − (since y ∈ S (z), and of course z − z
SR (z) 2
2
R
R2 1 φ(y) dσy = 3 2πR 2πR
SR (z)
= 0; therefore,
φ(y) dσy .
SR (z)
Thus u(z) is the average value of u over the circle S R (z) (recall that u(y) = φ(y) on S R (z)).
11.7
Green’s function for the wave equation
1. Let ψ : R 3
3
→ R, c ∈ R, and y ∈ R
be given. Then, by the chain rule,
∂ (ψ(x + cty)) = ∂t
3
∂ψ ∂ (x + cty) (xi + ctyi ) ∂x i ∂t
i=1 3
=
∂ψ (x + cty)(cyi ) ∂x i
i=1
3
= c
i=1
Similarly,
∇ 3
∂ 2 ∂ (ψ(x + cty)) = 2 ∂t ∂t
c
i=1
3
= c
∂ψ (x + cty)yi = ∂x i
i=1
∂ ∂t
3
= c
3
i=1 j=1 3
= c
3
i=1 j=1 3
= c
2
∂ψ (x + cty)yi ∂x i
∂ψ (x + cty)yi ∂x i
∂ 2 ψ ∂ (x + cty)yi (xj + ctyj ) ∂x j ∂x i ∂t ∂ 2 ψ (x + cty)yi cyj ∂x j ∂x i 3
yi
i=1
∇ψ(x + cty) · y.
j=1
∂ 2 ψ (x + cty)yj ∂x j ∂x i
3
= c2
yi
2
ψ(x + cty)y
i=1
= c2 y
2
· ∇ ψ(x + cty)y.
i
3. We wish to derive the solution of ∂ 2 u ∂t 2
2
2
− c ∆u = 0, x ∈ R , t > 0, u(x, 0) = 0, x ∈ R , ∂u (x, 0) = ψ(x), x ∈ R . ∂t 2
2
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 t u(x, t) = 4π
ψ(x + cty) dσy .
∂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 R 2 , is the solution of the IVP given above. Therefore, it remains only to write u as an integral over a two-dimensional region, namely, B 1 (0). If S represents the upper hemisphere of ∂ B1 (0), then
ψ(x + cty) dσy = 2
∂B1 (0)
√ √ − − ψ(x + cty) dσy
S
1 x2
1
−
ψ(x1 + cty1 , x2 + cty2 ) dy2 dy1 1 y12 y22 −1 − 1− ψ(x + cty) =2 dy. 1 y 2 B1 (0) =2
−
x2
It follows that
ψ(x + cty) dy. 1 y 2 B1 (0) Performing the change of variables z = x + cty yields the alternate formula t 2π
u(x, t) =
1 u(x, t) = 2πc 5. The solution of ∂ 2 u ∂t 2
−
Bct (x)
c 2 t2
ψ(z) 2
− z − x
2
dy.
3
− c ∆u = 0, x ∈ R , t > 0, u(x, 0) = 0, x ∈ R , ∂u (x, 0) = ψ(x), x ∈ R . ∂t 3
3
is
t u(x, t) = 4π(ct)2
Let ψ(x) =
ψ(z) dσz .
∂Bct (x)
1, x < 1, 0, otherwise,
and let u be the corresponding solution of the above IVP. Suppose first that x R3 , x > R, and t < ( v This last condition implies that ct < x R, and therefore, if z S ct (x), then
∈ ∈ z − x = ct < x − R. It is easy to see that this condition implies that z > R (since z is closer to x than x is to S z ∈ S (x) ⇒ ψ(z) = 0, −
R (0)).
ct
and it follows from the formula for u(x, t) that u(x, t) = 0. Now suppose t > ( x + R)/c, that is, ct > x + R. Then
∈ S (x) ⇒ z − x = ct > x + R, which implies that z > R (since z − x ≤ z + x). Therefore, in this case also, z ∈ S (x) ⇒ ψ(z) = 0, z
ct
ct
and we see that u(x, t) = 0. 7. The causal Green’s function is G(x, t; y, s) =
∞
n=1
√ − √ c λ
sin (c λn (t s)) n
0,
ψn (x)ψn (y),
t > s, t < s.
− R)/c.
Therefore,
109
11.8. GREEN’S FUNCTIONS FOR THE HEAT EQUATION
11.8
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−x
/(4kt)
satisfies the homogeneous heat equation in R n . We have 2 ∂S ( x, t) = (4kπt)−n/2 e−x /(4kt) ∂t 2 ∂S (x, t) = (4kπt)−n/2 e−x /(4kt) ∂x i 2 ∂ 2 S (x, t) = (4kπt)−n/2 e−x /(4kt) 2 ∂x i 2
∆S (x, t) = (4kπt)−n/2 e−x
/(4kt)
x 2 n , 4kt2 2t xi , 2kt x2i 1 , 4k2 t2 2kt
− − − − x 2 4k2 t2
n 2kt
.
From these formulas, it follows immediately that ∂S ( x, t) ∂t
− k∆S (x, t) = 0.
3. We have seen that the IBVP ∂v ∂t
− k∆v = f (x, T − t), x ∈ Ω, 0 < 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 ∂v (x, t) − k∆u(x, t) = (x, T − t) − k∆v(x, T − t) = f (x, T − (T − t)) = f (x, t). − ∂u ∂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 T
u
0
=
u
Ω
=
Ω T
0
uv
Ω
∂v ∂t
k∆v
T
v
0
∂u ∂t
k∆u
d x dt T
u∆v dx dt +
k
0
T . We have
t
T
∂v dt dx ∂t
T 0
v
Ω and 0
Ω
∂u dt d x ∂t
v
Ω
T
k
u
0 T
k
0
∂ Ω
0
∂u dt dx + k ∂t
(u(x, T )v(x, T )
u
Ω
∂ Ω
v dx dt +
Ω
u
v dx dt
Ω
v
0
Ω
Ω
∂u v dσx ∂ Ω ∂ n
u(x, 0)v(x, 0)) d x + k
v∆u dx dt
0
T
∂v dσ x ∂ n
T
=
T
∂u ∂ n
u
∂v dσ x . ∂ n
0
v
∂u dx dt+ ∂t
Chapter 12
More about Fourier Series 12.1
The complex Fourier series
1. The complex Fourier series of f is
∞
cn eiπnx ,
n=
−∞
where c 0 = 2/3 and cn =
− 2(n−π1)
n
, n =
2 2
±1, ±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.5
0 x
0.5
1
−0.5
0 x
0.5
1
−0.5
0 x
0.5
1
0.04 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) = a 0 +
∞
nπx nπx + bn sin
an cos
n=1
111
,
112
CHAPTER 12. MORE ABOUT FOURIER SERIES 1
0.5
0 −1
−0.5
0 x
0.5
1
−0.5
0 x
0.5
1
−0.5
0 x
0.5
1
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
=
1 g(x) dx, 2 − 1 nπx g(x)cos dx, n = 1, 2, 3, . . . , −
1 nπx g(x)sin dx, n = 1, 2, 3, . . . , −
(see (6.25) in the text), and similarly let h(x) = p 0 +
∞
nπx nπx + q n sin
pn cos
n=1
.
The complex Fourier coefficient of f is given by
1 cn = f (x)e−iπnx/ dx. 2 −
For n > 0, we have cn
= =
=
1 nπx nπx (g(x) + ih(x)) cos i sin dx 2 − 1 1 1 nπx nπx g(x)cos dx i g(x)sin dx 2 − − 1 1 nπx nπx +i h(x)cos dx + h(x)sin dx − − 1 (an + q n + i( pn bn )) . 2
− − −
Similarly, if n > 0, then c−n =
1 (an 2
− q + i( p
n +
n
bn )) .
Finally, c0 = a 0 + ip0 . 7. We have N 1
−
n=0
e
πijn/N
2
N 1
=
−
n=0 2πij
=
e
−1 , −1
e2πij/N
n
e2πij/N
113
12.2. 12.2. FOURIER FOURIER SERIES SERIES AND AND THE FFT FFT using the formula for the sum of a finite geometric series: m
rm+1 1 . r 1
rn =
−
n=0
−
Moreover, since e since e iθ is a 2π 2π-periodic -periodic function function of θ, θ , we see that e that e 2πij = 1 and hence that N 1
−
eπijn/N
n=0
2
= 0.
On the other hand, by Euler’s formula, N 1
−
n=0
2
eiπjn/N
N 1
−
2
− −
=
jnπ N
cos
n=0
N 1
−
=
jnπ N
cos2
n=0
N 1
−
=
cos
n=0
N 1
−
+2i +2i
N 1
−
n=0
sin2
n=0
jnπ N
cos
jnπ N
sin2
jnπ N
2
jnπ N
+ i sin
sin
+ 2i 2i cos
jnπ N
sin
jnπ N
jnπ N
jnπ . N
Since this expression equals zero, we must have N 1
−
N 1
−
− 2
cos
n=0
jnπ N
jnπ N
=
0,
sin
jnπ N
=
0.
n=0
N 1
−
sin2
cos
n=0
jnπ N
The second second equation equation is one of the result resultss we set out to prove prove.. The first equatio equation n can be written written,, using using the trigonometric identity sin2 (θ) = 1 cos2 (θ), as
−
N 1
−
N 1
−
− − cos2
n=0
jnπ N
1
cos2
n=0
jnπ N
which simplifies to N 1
−
2
cos2
n=0
or
N 1
−
cos2
n=0
which is the desired result. The result
N 1
−
sin2
n=0
then follows.
12.2 12.2
jnπ = N = N,, N
jnπ N = , 2 N jnπ N = N 2
Fourie ourier r seri series es and and the the FFT FFT
{| |}
1. The graph graph of F n is shown in Figure 12.3. 3. As shown shown in the text,
. cn = F n , n =
16, −15, 15, . . . , 15, 15, −16,
= 0,
114
CHAPTER CHAPTER 12. MORE ABOUT ABOUT FOURIER FOURIER SERIES SERIES 0.5
0.4
0.3 |
n
F |
0.2
0.1
0 0
5
10 n
15
20
Figure 12.3: The magnitude of the sequence { F n } 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 function f ((x) = x(1 x (1 − x2 ) (the curve) and the estimates of f of f ((xj ) (the circles) computed from the complex Fourier coefficients of f of f (see Exercise 12.2.3.
where F n
15 n= 16 is
{ }
−
the DFT of f j
15 j = 16 ,
{ { }
−
f j = f ( f (xj ), xj =
16, −15, 15, . . . , 15 −1 + 16j , j = −16,
(note that f that f (( 1) = f (1)). Therefore, we can estimate f ( f (1)). Therefore, f (xj ) by taking the inverse DFT (using the inverse FFT) of cn . The results are shown in Figure 12.4.
{ { }
−
{
}
5. We have have N 1
−
An e
2πijn/N
=
n=0
=
1 N
N 1 N 1
1 N
N 1 N 1
−
−
n=0 m=0
−
−
πi (j −m)n/N am e2πi(
n=0 m=0 N 1
=
am e−2πimn/N e2πijn/N
−
1 am N m=0
N 1
−
e
n=0
If j j = m, m, then N 1
−
n=0
N 1 πi(j −m)n/N = e2πi(
−
n=0
2πi( πi (j m)n/N
1 = N, N ,
−
.
12.3. 12.3. RELATIO RELATIONSHI NSHIP P OF SINE AND COSINE COSINE SERIES TO THE FULL FOURIER FOURIER SERIES
115
while if j j = m, m , then this sum is a finite geometric series: N 1
−
N 1
e
2πi( πi(j m)n/N
−
−
=
n=0
n=0
πi(j −m)/N e2πi(
=
πi(j −m)/N e2πi( πi (j −m) e2πi( 1
N
n
πi(j −m)/N e2πi(
−1
−1
− −1
=
πi(j −m)/N e2πi( 0.
= Therefore, N 1
am
−
e
2πi( πi(j m)n/N
−
=
n=0
= j, m = j, m = j,
N aj , 0,
and so we obtain N 1
−
N 1
An e
2πijn/N
n=0
N 1
−
1 = am N m=0
−
e
2πi( πi(j m)n/N
−
n=0
=
1 N aj = a j , N
as desired. 7. Below we show the exact Fourier Fourier sine coefficients of f of f , the coefficients estimated by the DST, and the relative error. Since both a both a n and the estimated a estimated a n are zero for n for n even, we show only a only a n for n for n odd. n
an
1 2.5801 · 10−1 3 9.5560 · 10−3 5 2.0641 · 10−3 7 7.5222 · 10−4 9 3.5393 · 10−4 11 1.9385 · 10−4 13 1.1744 · 10−4 15 7.6448 · 10−5
estimated a estimated a n
relative error
2.5801 · 10−1 9.5511 · 10−3 2.0555 · 10−3 7.3918 · 10−4 3.3531 · 10−4 1.6778 · 10−4 8.0874 · 10−5 2.4279 · 10−5
6.2121 10−6
· 5.1571 · 10− 4.1824 · 10− 1.7340 · 10− 5.2608 · 10− 1.3448 · 10− 3.1135 · 10− 6.8241 · 10−
4 3 2 2 1 1 1
9. We have have N 1
−
jnπ 2 F n sin N n=1
N 1 N 1
=
4
−
−
f m sin
n=1 m=1 N 1
=
4
−
N 1
f m
m=1
−
sin
n=1
mnπ jnπ sin N N
mnπ jnπ sin . N N
By the hint, N 1
−
sin
n=1
mnπ jnπ sin N N
is either N either N/ /2 or 0, depending on whether m whether m = j = j or not. It then follows immediately that N 1
2
−
n=1
which is what we wanted to prove.
F n sin
jnπ N = 4 f j = 2N f j , N 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 terms) 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 terms, 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 f odd is continuous on [ , ] only if f (0) = 0. Otherwise, f odd has a jump discontinuity at x = 0. (b) The even extension f even 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 ˜ defined by odd extension of the result. That is, f is
→
˜ f (x) =
−−
f (x), f (2 x), f ( x), f (x + 2),
− −
−
→
0 x , < x 2, x < 0, 2 x < .
≤ ≤ ≤ − ≤ − ≤ −
Since this function is odd, its full Fourier series has only sine terms (all of the cosine coefficients are zero), and because of the other symmetry, the even sine terms also drop out.
12.4
Pointwise convergence of Fourier series
1. The Fourier series of f converges pointwise to the function F defined by F (x) =
0, (x
3
− 2kπ) ,
x = (2k
±(2k − 1)π, k = 1, 2, 3, . . . , − 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 f (x) =
1 2
0
≤ x ≤ , < x ≤ 1.
x, x + 1,
1 2
5. If f N converges uniformly to f on [a, b], then it obviously converges pointwise (after all, the maximum difference, over x [a, b], between f N (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
{ } ∈
M N = max f (x)
{|
N (x)
− f
| : a ≤ x ≤ b}.
Then
|
≤
=
M N b
b
a
By hypothesis, M N
b
f (x)
2 N (x) dx
− f
|
a
→ 0 as N → ∞, which gives the desired result.
2 dx M N
√ − a.
117
12.5. UNIFORM CONVERGENCE OF FOURIER SERIES
→ R is continuous. Then f , the even, periodic extension of f to R, is defined by 2k ≤ x < (2k + 1), k = 0, ±1, ±2, . . . , f (x − 2k), f (x) = f (−x + 2k), (2k − 1) ≤ x < 2k, k = 0, ±1, ±2, . . . . Obviously, then, f 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 f (x) = lim f (−x + 2k) = lim f (x) = f (0) → → →
7. Suppose f : [0, ]
even
even
even
x
2k−
even
x
2k−
x
0+
and lim
→2k+
x
f even (x) =
lim f (x
→2k+
x
− 2k) =
lim f (x) = f (0).
→0+
x
Since f even (2k) = f (0) by definition, this shows that f even is continuous at x = 2κ. A similar calculation shows that x
lim
→(2k−1)−
f even (x) =
lim
→(2k−1)+
x
f even (x) = f even ((2k
− 1)) = f (),
and hence that f even is continuous at x = (2k 1). ˜ : [ 2, 2] 9. Given f : [0, ] R, define f R by
→
−
→
˜ = f (x)
−
−−
f (x), f (2 x), f ( x), f (x + 2),
− −
≤ ≤ ≤ − ≤ − ≤ −
0 x , < x 2, x < 0, 2 x < .
˜ to R . Assuming f is piecewise smooth, the quarter-wave Then define F : R R to be the periodic extension of f 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 f per also has discontinuities at x = 1). Therefore, the Fourier coefficients of f decay like 1/n2 . Finally, g per 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 terms, 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
∞
1 2 n n=1 converges to a finite value.1 Therefore, 1/n function in L 2 (0, 1).
2
{ } ∈
1A
standard result in the theory of infinite series is that
∞
n=1
1 nk
and so, by Theorem 12.36, the sine series converges to a
118
CHAPTER 12. MORE ABOUT FOURIER SERIES 1.5
1
y
0.5
0 21 terms 41 terms 81 terms 161 terms −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 terms, and the line y = 1.09. 1.105 21 terms 41 terms 81 terms 161 terms
1.1
1.095
y
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 terms 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
−
1
cn = 2
m(1
0
, − x)sin(nπx) = 2m nπ
n = 1, 2, 3, . . . .
Therefore, in order that c n = 1/n, we must have m = π/2. Therefore, f (x) =
π (1 2
− x).
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 terms, from Exercise 12.6.1.
∞
converges if k is greater than 1. This can be proved, for example, by comparison with the improper integral
1
dx .
xk
119
12.7. A NOTE ABOUT GENERAL EIGENVALUE PROBLEMS 3. The infinite series
∞
2
∞
√ 1 n
n=1
=
1 n n=1
{ √ } ∈
is the harmonic series , a standard example of a divergent series. Therefore, 1/ n 2 , and so the series does not converge to an L 2 (0, 1) function. 1 5. (a) The function f (x) = x k belongs to L 2 (0, 1) if and only if the integral 0 x2k dx is finite. Provided k = 1/2, we have 1 1 1 1 r2k+1 , k > 12 , 1+2k x2k dx = lim x2k dx = lim = , k < 12 . 1 + 2k r→0+ r r→0+ 0 For k = 1/2, 1 1 1 dx dx = lim = lim ln (r) = . x2k dx = + x x r→0 r→0+ 0 0 r Thus we see that f L2 (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 terms, and the difference between the two, are shown in Figure 12.9. 6
y=x−1/4 sine series (63 terms)
4 y
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 terms) (top), and the difference between the two (bottom). See Exercise 12.6.5.
→ v, there exists a positive integer N such that n ≥ N ⇒ v − v < . 2 Then, if n, m ≥ N , we have v − v ≤ v − v + v − v < 2 + 2 = . Therefore, {v } is Cauchy.
7. Since vn
n
n
m
n
m
n
12.7
A note about general eigenvalue problems (K D u, v) =
− ∇· Ω
(k(x) u) v
∇
=
∇ · ∇ − ∇ · ∇ − ∇· ∇ − ∇· ∇ k(x) u
Ω
=
k(x)v
v
∂ Ω
k(x) u
∂u ∂ n
v (since v = 0 on ∂ Ω)
Ω
=
(k(x) v) u +
Ω
=
∂ Ω
∂v ∂ n
(k(x) v) u (since u = 0 on ∂ Ω)
Ω
=
k(x)u
(u, K D v).
2 D (Ω).
∈ C
1. The proof is a direct calculation, using integration by parts (Green’s identity): Suppose u, v
Then
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 c 1 is the first generalized Fourier coefficient of the initial temperature 5: c1 =
5ψ1 . ψ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 terms 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) = x 2 . With f (x) = 1, we have 1 f = (area of T R ) 1 = 2 T R
T R
f when f (x) = 1, f (x) = x 1 , or
·
and
1 1 1 1 f , = . 2 3 3 2
With f (x) = x 1 , we have 1
=
f
T R
1 x1
−
− · 0
x1 dx2 dx1
0
1
=
x21 dx 1 =
x1
0
On the other hand,
1 2
− 13 = 61 .
1 1 1 1 1 1 = = , f , 2 3 3 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 , 1 1−x1 1 1 1 1 f = x21 dx2 dx1 = x21 x31 dx 1 = = , 3 4 12 T R 0 0 0
while
·
−
−
1 1 1 1 1 1 f , = = . 2 3 3 2 9 18 3. We have
3/2
2x1
f =
T
1
−2
2
x1 x22 dx2 dx1 +
0
On the other hand, the mapping from T
3/2
R
4 2x1
−
x1 x22 dx2 dx1 =
0
to T is x1
=
1 + z 1 +
x2
=
z 2 .
The Jacobian is J =
121
1 0
1 2
1
1 z 2 , 2
,
7 1 1 + = . 120 15 8
122
CHAPTER 13. MORE ABOUT FINITE ELEMENT METHODS and its determinant is 1. Therefore, with
1 g(z 1 , z 2 ) = f 1 + z 1 + z 2 , z 2 = 2 we have
1
1 z1
−
1 1 + z 1 + z 2 z 22 , 2
f =
T
g =
T R
0
1 1 + z 1 + z 2 z 22 dz 2 dz 1 . 2
0
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 ci ; let us define j dj by
×
→
dj =
→
if j = c i , if j = c i for all i = 1, 2, . . . , K .
i 0
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 g i 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 b oundary). 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
−
F i =
f φi
Ω
by F i =
f φi
Ω
− a(G, φ ). i
The quantity a(G, φi ) is nonzero only if the free node nf i 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
−
n 1
i=1
2(n
2
− i)
+ (n
−
n 1
−
− i)
= 2
2
j +
j=1
j =
j=1
2n3 3
2
− n2 − n6 .
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
−
−
n 1
2
−
n 1
(n
i=1
− 1) = 2
−
j = n 2
j=1
The final step of back substitution requires n steps, with 2(n n
− n.
− i) + 1 operations per step, for a total of
n 1
{
2(n
i=1
The total for the two triangular solves is 2n2
− i) + 1 } = 2
−
j + n = n 2 .
j=1
− n, which also agrees with the count given in the text.
123
13.2. SOLVING SPARSE LINEAR SYSTEMS 3. With 1 x Ax 2
·
φ(x) =
− b · x,
we have φ(x + y) =
1 (x + y) A(x + y) 2
·
− b · (x + y)
1 1 x Ax + y Ax + y Ay b x 2 2 1 φ(x) + ( Ax b) y + y Ay. 2
·
= =
· − ·
·
− · −b·y
·
The fact that A is symmetric allows the simplification 1 1 x Ay + y Ax = y Ax, 2 2
·
·
·
2
∇φ(x) · y + O y
which was used above. We thus see that φ(x + y) = φ(x) + 5. An inner product has three properties: (a) (x, x)
∇φ(x) = Ax − b.
, with
≥ 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,
−∇φ
φ x(0)
α φ x(0)
− −1), and so we must minimize
= 7α2
− 5α − 18.
The minimizer is α = 5/14, and so (1)
x
=
23 23 , 7 14
.
(b) At x (1) = (23/7, 23/14), the negative gradient is ( 3/14, 3/7), so we must minimize
−
φ x(1)
− α∇φ
(1)
x
=
27 2 α 196
45 529 − 196 α− . 28
The minimizer is α = 5/6, and so (2)
x
=
87 ,2 . 28
The desired solution is (3, 2). (c) With x (0) = 0 , the CG algorithm produces (1)
x
. = (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
We can therefore replace b by b of x.
(0)
− Ax
(0)
⇒ Ay + Ax
= b
(0)
⇒ Ay = b − Ax
.
, apply the CG algorithm to estimate y , and add x (0) to get the estimate
