• Author / Uploaded
• gustavo leiziaga

Citation preview

Undergraduate Texts in Mathematics Editors

S. Axler F. W. Gehring K.A. Ribet

Springer New York Berlin Heidelberg Barcelona Budapest Hong Kong London Milan Paris Singapore Tokyo

Undergraduate Texts in Mathematics Anglin: Mathematics: A Concise History and Philosophy. Readings in Mathematics. Anglin/Lambek: The Heritage of Thales. Readings in Mathematics. Apostol: Introduction to Analytic Number Theory. Second edition. Armstrong: Basic Topology. Armstrong: Groups and Symmetry. Axler: Linear Algebra Done Right , Second edition. Beardon: Limits: A New Approach to Real Analysis. BakINewman: Complex Analysis. Second edition. BanchofflWermer: Linear Algebra Through Geometry. Second edition. Berberian: A First Course in Real Analysis. Bix: Conics and Cubics: A Concrete Introduction to Algebraic Curves. Bremaud: An Introduction to Probabilistic Modeling. Bressoud: Factorization and Primality Testing. Bressoud: Second Year Calculus. Readings in Mathematics. Brickman: Mathematical Introduction to Linear Programming and Game Theory. Browder: Mathematical Analysis: An Introduction. Buskes/van Rooij: Topological Spaces: From Distance to Neighborhood. Cederberg: A Course in Modem Geometries. Childs: A Concrete Introduction to Higher Algebra. Second edition. Chung: Elementary Probability Theory with Stochastic Processes. Third edition. Cox/Little/O'Shea: Ideals, Varieties, and Algorithms. Second edition. Croom: Basic Concepts of Algebraic Topology. Curtis: Linear Algebra: An Introductory Approach. Fourth edition.

Devlin: The Joy of Sets: Fundamentals of Contemporary Set Theory. Second edition. Dixmier: General Topology. Driver: Why Math? Ebbinghaus/Flum/Thomas: Mathematical Logic. Second edition. Edgar: Measure, Topology, and Fractal Geometry. Elaydi: Introduction to Difference Equations. Exner: An Accompaniment to Higher Mathematics. Fine/Rosenberger: The Fundamental Theory of Algebra. Fischer: Intermediate Real Analysis. FJanigan/Kazdan: Calculus Two: Linear and Nonlinear Functions. Second edition. Fleming: Functions of Several Variables. Second edition. Foulds: Combinatorial Optimization for Undergraduates. Foulds: Optimization Techniques: An Introduction. Franklin: Methods of Mathematical Economics. Gordon: Discrete Probability. HairerlWanner: Analysis by Its History. Readings in Mathematics. Halmos: Finite-Dimensional Vector Spaces. Second edition. Halmos: Naive Set Theory. Hlimmerlin/Hoffmann: Numerical Mathematics. Readings in Mathematics. Hijab: Introduction to Calculus and Classical Analysis. HiltonlHolton/Pedersen: Mathematical Reflections: In a Room with Many Mirrors. Iooss/Joseph: Elementary Stability and Bifurcation Theory. Second edition. Isaac: The Pleasures of Probability. Readings in Mathematics. (continued ofter index)

J. David Logan

A pplied Partial

Differential Equations With 35 Illustrations

Springer

J. David Logan Department of Mathematics and Statistics University of Nebraska at Lincoln Lincoln, NE 68588-0323 USA Editorial Board

S. Axler Mathematics Department San Francisco State University San Francisco, CA 94132 USA

F.W. Gehring Mathematics Department East Hall University of Michigan Ann Arbor, MI 48109 USA

K.A. Ribet

Department of Mathematics University of California at Berkeley Berkeley, CA 94720-3840 USA

Mathematics Subject Classification (1991): 35-01

Library of Congress Cataloging-in-Publication Data Logan, J. David (John David) Applied partial differential equations / J. David Logan. p. cm.-(Undergraduate texts in mathematics) Includes bibliographical references and index. ISBN -13: 978-0-387-98439-1 e-ISBN -13: 978-1-4684-0533-0 DOl: 10 ·1007/978-1-4684-0533-0 1. Differential equations, Partial. I. Title. II. Series. QA377.L578 1998 515'.353-dc21 97-48861 Printed on acid-free paper.

© 1998 Springer-Verlag New York, Inc. All rights reserved. This work may not be translated or copied in whole or in part without the written permission of the publisher (Springer-Verlag New York, Inc., 175 Fifth Avenue, New York, NY 10010, USA), except for brief excerpts in connection with reviews or scholarly analysis. Use in connection with any form of information storage and retrieval, electronic adaptation, computer software, or by similar or dissimilar methodology now known or hereafter developed is forbidden. The use of general descriptive names, trade names, trademarks, etc., in this publication, even if the former are not especially identified, is not to be taken as a sign that such names, . as understood by the 1l:ade Marks and Merchandise Marks Act, may accordingly be used freely by anyone. Production managed by Victoria Evarretta; manufacturing supervised by Jacqui Ashri. Photocomposed copy prepared from the author's files by The Bartlett Press, Inc.

9 8 7 6 5 4 3 2 1 ISBN-13: 978-0-387-98439-1

Preface

This textbook is for the standard, one-semester, junior-senior course that often goes by the title "Elementary Partial Differential Equations" or "Boundary Value Problems;' The audience usually consists of students in mathematics, engineering, and the physical sciences. The topics include derivations of some of the standard equations of mathematical physics (including the heat equation, the· wave equation, and the Laplace's equation) and methods for solving those equations on bounded and unbounded domains. Methods include eigenfunction expansions or separation of variables, and methods based on Fourier and Laplace transforms. Prerequisites include calculus and a post-calculus differential equations course. There are several excellent texts for this course, so one can legitimately ask why one would wish to write another. A survey of the content of the existing titles shows that their scope is broad and the analysis detailed; and they often exceed five hundred pages in length. These books generally have enough material for two, three, or even four semesters. Yet, many undergraduate courses are one-semester courses. The author has often felt that students become a little uncomfortable when an instructor jumps around in a long volume searching for the right topics, or only partially covers some topics; but they are secure in completely mastering a short, well-defined introduction. This text was written to proVide a brief, one-semester introduction to partial differential equations. It is limited in both scope and depth compared with existing books, yet it covers the main topics usually studied in the standard course and also provides an introduction to using computer algebra packages to solve and understand partial differential equations. The frontiers of mathematics and science v

vi

Preface

are receding rapidly, and a one-semester course must try to advance the students to a level where they can reach these boundaries more quickly than in the past. Not every traditional topic can be covered, and not every topic can be examined in great detail. An example is the method of separation of variables, which plays a dominant role in most texts. It is this author's view that a few well-chosen illustrations of the method of separation of variables should suffice. The level of exposition in this text is slightly higher than one usually encounters in the post-calculus differential equations course. The philosophy here is that a student should progress in his or her ability to read mathematics. Elementary calculus texts, and even ordinary differential equations texts, contain lots of examples and detailed calculations, but advanced mathematics and science books leave a lot to the reader. This text leaves some of the details that are easy to supply to the reader. The student is encouraged as part of the learning process to fill in these missing details (see I/Tb the Student"). The writing has more of an engineering and science style to it than a traditional, mathematical, theorem-proof format. Consequently, the arguments given are I/derivations" in lieu of carefully constructed proofs. Students who come out of reform calculus courses often have different skills from those of students who have had a traditional calculus course. For example, they rely more on calculators and computers in solving some of their problems. This text is a little less traditional than other elementary texts on partial differential equations. It encourages the student to use software packages, where applicable, in solving problems. Herein we use Maple, version V4, for the illustrations, but the text is in no way dependent on Maple; Mathematica, or any other computer algebra program, can be used. The exercises encourage students to think about the concepts and derivations rather than just grind out lots of routine solutions. The student who reads this book carefully and who solves most of the exercises will have a sound enough knowledge base to continue with a second-year partial differential equations course where careful proofs are constructed or with upper-division courses in science and engineering where detailed applications of partial differential equations are introduced. Both the exposition and exercises will build analytical skills that some students did not develop in the reform calculus courses. In Chapter 1 we introduce some basic partial differential equations of applied mathematics. Many of the basic equations come from a conservation law, or balance law, and describe physical processes like convection, diffuSion, and reaction. There are a variety of applications to show the central role that partial differential equations play in science, engineering, and mathematics. Applications include contaminant transport in aquifers, convection and diffusion processes in biology and chemical engineering, quantum mechanics, mechanical vibrations, heat flow, electromagnetic phenomena, and acoustics. The chapter ends with a dis-

Preface

vii

cussion of classification of partial differential equations; the emphasis is on the physical meaning of the equations (diffusion-like, wave-like, etc.) as well as on the analytical structure of the equations. The goal of the chapter is to give students a sense of the origins of partial differential equations and how their solutions differ. At the same time the exercises force the students to revisit the chain rule, the divergence theorem, and other concepts from multivariable calculus. Chapter 2 examines equations on unbounded domains, both infinite and semi-infinite. The author's view is that these problems are simpler than their counterparts on bounded domains with boundaries present. Easy formulae are obtained for solutions to the initial value problems for the heat and wave equations, and well-posedness is discussed. Problems with sources are handled with Duhamel's principle. Most students have studied Laplace transforms in an elementary ordinary differential equations course, so it is a natural transition to study transform methods for partial differential equations; the Fourier transform is a straightforward extension from the Laplace transform. We also include material on using computer algebra packages to find general solutions to equations and solutions based on transform methods. One of the fundamental ideas in applied mathematics is orthogonality. In Chapter 3, rather than adopt a strict focus on Fourier series, a general strategy is taken. Calculus courses have always included- Thylor series, and many calculus courses, especially reform courses, now include some material on Fourier series. Therefore, students are ready to be introduced to general expansions of functions in series, especially orthogonal series. These expansions are motivated by the separation of variables method, and classical Fourier series are studied as a special case. The chapter ends with a study of Sturm-Liouville problems. Chapter 4 contains traditional material on the separation of variables method for solving partial differential equations on bounded domains. Here we solve various equations with various boundary conditions in rectangular, cylindrical, and spherical geometry. Students are urged to use software packages to perform some of the calculations. There is a section on inverse problems and a long section on the finite difference method. I want to thank my wife, Tess, and my son, David, for their encouragement and support. I also thank my many students for enduring preliminary versions of this text and especially for their comments on the clarity of the exposition and suitability of some of the exercises.

viii

Preface

Suggestions for Using the 'Text: The author has taught this material on numerous occasions and uses approximately the following schedule: Chapter 1 (12 classes), Chapter 2 (11 classes), Chapter 3 (7 classes), Chapter 4 (14 classes). Under this schedule, the sections marked with an asterisk (*) in the Table of Contents are sometimes not covered in lectures. Often those sections are assigned as extra reading material to graduate students taking the course.

Lincoln, Nebraska

J. David Logan

To the Student

Partial differential equations (PDEs) is a subject about differential equations for unknown functions of several variables; the derivatives involved are partial derivatives. As such, it is a subject that is intimately connected with multivariable or third-semester calculus. Th be successful you should have, first, a good command of the concepts in the calculus of several variables. So keep a calculus text nearby and review concepts when they are needed. The same comments apply to elementary ordinary differential equations (ODEs). There is an appendix at the end of the ·book that reviews some ofthe basic solution techniques for ODEs. Second, a mathematics book should be read with a pencil and paper at hand. Elementary books fill in most ofthe steps in the exposition, but advanced books leave many details to the reader. This book has enough detail so that you can follow the discussion, but pencil and paper work is required in some portions. VerifYing all the statements in a text is a worthwhile endeavor and will help you learn the material. Many students find that studying PDEs provides an opportunity to reinforce many calculus concepts and calculations. Finally, the exercises are an important part of this text (perhaps the most important part!), and you should try to solve most or all of them. Some will require routine analytical or computer calculations, but others will require careful thought. We learn mathematics by doing mathematics, even when we are stymied by a problem. The effort put into a failed attempt will help you sort out the concepts and reinforce the learning process. View the exercises as a challenge and resist the temptation to give up. Lincoln, Nebraska

J. David Logan

ix

Contents

Preface .............................................................. v 'Ib the Student ..................................................... ix

Chapter 1: The Physical Origins of Partial 'Differential Equations ........................................................... 1 1.1 1.2 1.3 1.4 1.5 l.6 1.7 1.8 l.9 l.10

Mathematical Models ........................................ 1 Conservation Laws ........................................... 9 Diffusion ................................................... 15 Contaminant 1tansport in Aquifers * ........................ 20 Vibrations of a String " ...................................... 24 Quantum Mechanics* ....................................... 28 Heat Flow in Three Dimensions ............................. 31 Laplace's Equation .......................................... 36 Acoustics* .................................................. 42 Classification of PDEs ....................................... 44

Chapter 2: Partial Differential Equations on Unbounded Domains ......... ............................................. ..... 50 2.1 2.2 2.3 2.4

Cauchy Problem for the Heat Equation ...................... 50 Cauchy Problem for the Wave Equation ..................... 56 Ill-Posed Problems .......................................... 61 Semi-Infinite Domains ...................................... 63 xi

xii

2.5 2.6 2.7 2.8

Contents

Sources and Duhamel's Principle ............................ 68 Laplace 1tansforms ......................................... 72 Fourier 1tansforms ......................................... 78 Solving PDEs Using Computer Algebra Packages ............ 84

Chapter 3: Orthogonal Expansions ............................... 91 3.1 3.2 3.3 3.4

The Fourier Method ......................................... 91 Orthogonal Expansions ..................................... 94 Classical Fourier Series .................................... 102 Sturm-Liouville Problems .................................. 108

Chapter 4: Partial Differential Equations on Bounded Domains ................................................... '....... 116 4.1 4.2 4.3 4.4 4.5 4.6 4.7 4.8

Separation of Variables ..................................... 116 Flux and Radiation Conditions ............................. 124 Laplace's Equation ......................................... 131 Cooling of a Sphere ......................................... l39 Diffusion in a Disk ......................................... 144 Sources on Bounded Domains .............................. 149 Parameter Identification Problems* ........................ 152 Finite Difference Methods* ................................ 157

Appendix: Ordinary Differential Equa,tions .................... 169 Thble of Laplace Transforms ..................................... 175 References ........................................................ 177 Index .............................................................. 179

CHAPTER

1.1

The Physical Origins of Partial Differential Equations

Mathematical Models

Many important ideas of mathematics were developed within the framework of physical science. For example, calculus has its origins in efforts to accurately describe the motion of bodies. Mathematical equations have always provided a context in which to formulate concepts in physicsMaxwell's equations describe electrodynamical phenomena, Newton's equations describe mechanical systems, Schrbdinger's equation describes aspects of quantum mechanics, and so on. In recent years, however, mathematicians and scientists have extended these types of connections to include nearly all areas of science and technology, and a new field has emerged called mathematical modeling. A mathematical model is an equation, or set of equations, whose solution describes the physical behavior of a related physical system. In this context, we say, for example, that Maxwell's equations are a model for electrodynamical phenomena. Like most mathematical models, Maxwell's equations are based on experiment and physical observations. In general, a mathematical model is a simplified description of physical reality expressed in mathematical terms. Mathematical modeling involves physical observation, selection of the relevant physical variables, formulation of the equations, analysis of the equations and simulation, and, finally, validation of the model. In this last step information from the simulations and solutions is fed back into the model to ascertain whether indeed the model does describe the phenomenon; at this step, modifications can occur and refinements can be made. J. D. Logan, Applied Partial Differential Equations © Springer-Verlag New York, Inc. 1998

1

2

1.

The Physical Origins of Partial Differential Equations

In this book we study partial differential equation (PDE) models. That is, we examine physical phenomena that can be modeled by equations called partial differential equations (PDEs). The focus is on the origin of such models and the tools used for their analysis. Often such models are called continuous models, or infinite-dimensional models. The reader should be familiar with finite models, or physical systems governed by ordinary differential equations (ODEs). For example, a typical ODE model is the Malthus population equation du

dt = ru,

t > 0,

which is a simple model of population growth where the time rate of change of population u = u(t) is proportional to the population. Here t is time, and u = u(t) is the population of a given system of individuals. We refer to u as the state variable and say that the evolution of the state variable is governed by the model equation. The real number r is a given physical parameter that represents the relative growth rate; presumably, r could be measured for the given population under investigation. The solution to the model equation is easily found to be u(t) = uoe rt ,

t > 0,

where Uo represents the initial population. Thus, the Malthus model predicts exponential growth, which accurately describes some populations during beginning stages of growth. The Malthus model is typical of ODE models. The state variable is a function of a single independent variable (time t), and the model contains a parameter, in this case r, that characterizes the physical system that is being modeled. In general, an ODE model has the form du

dt = F(t, u; r),

t > 0,

where F is a given functional relation between t, u, and r. Often the model includes an initial condition of the form u(O) = uo, where Uo is a given state value at t = O. More generally, an ODE model may consist of a system ofn ODEs for n state variables Ul(t), ... , un(t), and there can be m parameters rl, ... , rm. A PDE model differs from an ODE model in that the state variable depends on more than one independent variable, and the resulting model equation is a PDE. Whereas an ODE may model the evolution of a system in time, and observations are made in time, a PDE may model the evolution of a system in both time and space; the system can be observed both in a time interval and in a spatial region (which may be one-, two-, or three-dimensional). PDE models may also be independent of time, but depend on several spatial variables. To fix the notion, let us consider the

1.1.

Mathematical Models

3

problem of determining the temperature in a laterally insulated metal bar oflength 1 and unit cross-sectional area, whose two ends are maintained at a constant zero degrees and whose temperature initially (at time zero) varies along the bar and is given by a fixed function O. Using this fact, construct a solution to the two-dimensional Laplace's equation in y > 0 that also satisfies the boundary conditions u(x, 0) = 1 for x > 0 and u(x, 0) = -1 for x < O. Be sure to explain which branch of the arctan function you are using. 9. Show that e-~Y sin(~x), x E R, Y > 0, is a solution to value of the parameter~. Deduce that u(x, y)

U xx

+ Uyy

= 0 for any

= 1':>0 c(~)e-~y sin(~x)d~

is a solution to the same equation for any function c(~) that is bounded and continuous on [0, (0). (The hypotheses on c allow you to differentiate under the integral sign.) This exercise shows that taking integrals of solutions sometimes gives another solution; integration is a way of superimposing, or adding, a continuum of solutions.

1.2.

1.2

Conservation Laws

9

Conservation Laws

Many PDE models come from a basic balance law, or conservation law. A conservation law is just a mathematical formulation of the basic fact that the rate at which a quantity changes in a given domain must equal the rate at which the quantity flows across the boundary plus the rate at which the quantity is created, or destroyed, within the domain. For example, consider a population of a certain animal species in a fixed geographical region. The rate of change of the animal population must equal the rate at which animals migrate into the region, minus the rate at which they migrate out, plus the birth rate, minus the death rate. Such a statement is a verbal expression of a balance, or conservation, law. One can make similar kinds of statements for many quantities-heat energy, the concentration of a chemical, the density of automobiles on a freeway, and so on. To quantify such statements we require some notation. Let the state variable u = u(x, t) denote the density of a given quantity (mass, energy, animals, automobiles, etc.); density is usually measured in amount per unit volume, or sometimes amount per unit length. For example, energy density is measured in energy units per volume. We assume that any variation in the quantity be restricted to one spatial dimension, that is, we assume a one-dimensional domain (say, a tube, as in Figure l.3) where each cross-section is labeled by the spatial variable x; we require that there be no variation of u(x, t) within the cross-section at x. Implicit is the assumption that the quantity in the tube is abundant and continuous enough in x so that it makes sense to define a density at each section of the tube. The amount of the quantity in a small section of width dx is u(x, t)Adx, where A is the cross-sectional area of the tube. Further, we let ¢ = ¢(x, t) denote the flux of the quantity at x, at time t. The flux measures the amount of the quantity crossing the section at x at time t, and its units are given in amount per unit area, per unit time. Thus, A¢(x, t) is the actual amount of the quantity that is crossing the section at x at time t. By convention, flux is positive if the flow is to the right, and

, I

:A \ \

I X

I

a

I

b

x

Figure 1.3. Thbe with cross-sectional area A shown with arbitrary section . The lateral sides are insulated, and quantities vary only in the x-direction.

1. The Physical Origins of Partial Differential Equations

10

negative if the flow is to the left. Finally, let [ = fCx, t) denote the given rate at which the quantity is created, or destroyed, within the section at x at time t. The function [ is called a source term if it is positive, and a sink if it is negative; it is measured in amount per unit volume per unit time. Thus, fCx, t)Adx represents the amount of the quantity that is created in a small width dx per unit time. A conservation law is a quantitative relation between u, ¢, and f. We can formulate the law by considering a fixed, but arbitrary, section a ~ x ~ b of the tube (Figure 1.3) and requiring that the rate of change of the total amount of the quantity in the section must equal the rate at which it flows in at x = a, minus the rate at which it flows out at x = b, plus the rate at which it is created within a ~ x ~ b. In mathematical symbols, :t

lab u(x, t)Adx = A¢(a, t) -

A¢(b, t)

+ lab fCx, t)Adx.

(1.7)

This equation is the fundamental conservation law; it is an integral expression of the basic fact that there must be a balance between how much goes in, how much goes out, and how much is changed. Note that because A is constant, it may be canceled from the formula. Equation (1.7) is an integral model. However, if the functions u and ¢ are sufficiently smooth, then it may be reformulated as a PDE model. For example, if u has continuous first partial derivatives, then the time derivative on the left side of (1.7) may be brought under the integral sign to obtain d dt

lb

a u(x, t)dx =

lb

a Ut(x, t)dx

.

(note that this is a special case of Leibniz's rule-see Section 1.1, Exercise 4). If ¢ has continuous first partials, then the fundamental theorem of calculus can be applied to write the change in flux as the integral of a derivative, or ¢(a, t) - ¢(b, t) =

-l

b

¢x(x, t)dx.

Therefore, (1.7) may be written

lb

(Ut(x, t)

+ ¢x(x, t) -

fCx, t))dx = O.

Because a ~ x ~ b can be any interval whatsoever, it follows that the integrand must vanish identically, or Ut(x, t)

+ ¢x(x, t) = [(x, t).

(1.8)

Equation (1.8) is a local version of (1.7), obtained under the assumption that U and ¢ are continuously differentiable; it is a PDE model describing

1.2.

Conservation Laws

11

the relation between the density, its flux, and the rate at which it is created. We shall call the PDE (1.8) the fundamental conservation law. The [-term is called the source term, and the ¢-term is called the flux term. In (1.7) we usually drop the understood notational dependence on x and t and just write Ut + ¢x = f. Before studying some examples, we make some general comments. The flux ¢ and source [ are functions of x and t, but their dependence on x and t may be through dependence upon the density U itself. For example, the source term [ may be given as a function of density via [ = [(U) , where, of course, U = u(x, t). Similarly, ¢ may depend on u. These dependencies may lead to a nonlinear model. Next, we observe that (1.8) is a single equation, yet there are two unknowns, U and ¢ (usually the source [ is assumed to be given). This means that another equation is required, which relates U and ¢. Such equations are called constitutive relations (or equations of state), and they arise from physical assumptions about the medium .. EXAMPLE

(Convection) A model where the flux is proportional to the density itself, that is, ¢

=

CU,

where c is a constant, is called a convection model. Notice that c must have velocity units (length per time). In this case the conservation law (1.8) becomes, in the absence of sources Cf = 0), Ut

+ CUx = o.

(1.9)

Equation (1.9) is called the convection equation. The reader should verify that the function u(x, t) = F(x - ct)

(1.10)

is a solution to (1.9) for any differentiable functionF. Such solutions (1.10) are called right-traveling waves, since the graph of F(x - ct) is just the graph of F(x) shifted to the right ct spatial units. So, as time t increases, the wave profile F(x) moves to the right, undistorted, with its shape unchanged, at speed c. Figure 1.4 shows two ways of viewing a traveling wave. Intuitively, (1.9) describes what we usually call convection. For example, a swarm of insects moving as an undistorted density wave would represent convection. Other common descriptive terms for this kind of movement are transport and advection. 0 If the flux is a nonlinear function of the density, that is, ¢ = ¢(u), then the conservation law (1.8) (again take [ = 0) takes the form Ut

+ ¢(u)x = Ut + ¢'(u)ux = o.

(1.11 )

1. The Physical Origins of Partial Differential Equations

12

(a)

(b)

u

u

F(x)

o

F(x - ct)

Yf\

~ I

o

x

x

Figure 1.4. Two views of a traveling wave: (a) time snapshots, and (b) in space-time as a ridge. If 0;

E

u(x, 0)

= e_,,2 ,x E R.

Pick c = 2 and sketch the solution surface and several time snapshots. Do you see a traveling wave? 3. Find the general solution of the convection-decay equation (1.12) by changing the independent variables to

=x -

~

ct,

r

= t.

4. Show that the decay term in the convection-decay equation (1.12) can be removed by making a change of dependent variables w = ueAt. 5. Solve the initial boundary value problem Ut

+ CU" = -AU, U(x,O)

= 0,

x, t > 0,

x > 0;

U(O, t) = get),

t > O.

In this problem you will have to treat the domains x > ct and x < ct differently; the boundary condition affects the region x < ct, and the initial condition affects the region x > ct. 6. Solve the pure initial value problem Ut

+ u" -

3u = t, x

U(x, 0)

= x2,

E

R, t > 0,

X E R.

7. 1b study the absorption of nutrients by a grasshopper we model its digestive tract by a tube of length 1 and cross-sectional area A. Nutrients of concentration n = n(x, t) flow through the tract at speed c, and they are adsorbed

1.3.

Diffusion

15

locally at a rate proportional to ..;n. What is the PDE model? If the tract is empty at t = 0 and then nutrients are introduced at the constant concentration no at the mouth (x = 0) for t > 0, formulate an initial boundary value problem for n = n(x, t). Solve this PDE model and sketch a graph of the nutrient concentration exiting the tract (at x = 1) for t > O. Physically, why is the solution u(x, t) = 0 for x > ct? 8. Explain why the function u(x, t) = G(x + ct) is called a left-traveling wave. How would you attempt to solve the convection equation Ut -CUx = F(x, t, u)? 9. The density of cars on a busy one-lane freeway with no exits and entrances is u = u(x, t) cars per mile. If ¢ = ¢(x, t) is the flux of cars, measured in cars per hour, derive a conservation law relating the density and flux. Why would ¢ = au(f3 - u) (a, f3 > 0) be a reasonable assumption? Write down the resulting nonlinear PDE for u. 10. Find a formula that implicitly defines the solution u = u(x, t) of the initial value problem for the reaction-convection equation Ut

+ vUx

= - f3

u(x, 0) = [(x),

au

+ u' X

X E

R, t > 0,

E R.

Here, v, a, and f3 are positive constants. Show from the implicit formula that you can always solve for u in terms of x and t. Hint: Use graphical methods.

1.3

Diffusion

Let us again write down the basic conservation law (1.8) with no sources: Ut

+ ¢x = o.

(1.16)

To reiterate, U = u(x, t) represents the density of a physical quantity, and ¢ = ¢(x, t) represents its flux. Equation (1.16) describes locally how changes in density are related to changes in flux. In the last section we modeled convection by assuming that the flux was proportional to the density (flux equals velocity times density). Now we want to model a simple diffusion process. To fix the notion let u denote the concentration of some chemical species, say a gas, in a tube. We expect that the random motion and collisions of the molecules will cause concentrations of the gas to spread out; the gas will move from higher concentrations to lower concentrations. The same could be said for insects in a tube, people congregated in a hallway, or heat energy in a metal bar. To model this type of random motion we make two observations: (i) the movement is from higher concentrations to lower concentrations, and (ii) the steeper the concentration gradient, the greater the flux. Therefore, the flux should depend on the x-derivative ofthe density (which measures the

16

1.

The Physical Origins of Partial Differential Equations

steepness of the density curve). Assuming a simple linear relationship, we take (1.17) where D is a constant of proportionality. The minus sign guarantees that if U x < 0, then ¢ will be positive and the flow will be, by our convention, to the right; if U x > 0, then ¢ will be negative and the flow will be to the left. We say that the flow is down the gradient. Equation (1.17) is called Pick's law, and the constant D is called the diffusion constant; D is measured in length-squared per unit time. When the constitutive equation (1.17) is substituted into the conservation law (1.16), we obtain a simple model equation Ut -

Duxx

= 0,

(1.18)

which is called the diffusion equation. This PDE model is one of the fundamental equations in applied mathematics. EXAMPLE

(Heat Flow) Let us consider heat flow in a one-dimensional bar having a constant density p and constant specific heat C. Both of these constants are physical parameters that are tabulated in engineering and physics handbooks. The specific heat is the amount of energy required to raise a unit mass of material one degree, and it is given in units of energy per mass per degree. We may apply the basic conservation law (1.16) to the bar, with u being the energy density given by u(x, t) = pCO(x, t), where o = O(x, t) is the temperature at (x, t) (the stipulation that energy is proportional to temperature is, in itself, an assumption about the medium). Therefore, pCOt

+ ¢x = 0

(1.19)

is an expression of energy balance in the bar when no sources are present. The energy flux ¢ is assumed to be given by a Fick's law-type expression (1.20) where K is thermal conductivity, another physical constant. In the context of heat flow, (1.20) is called the Fourier heat law. This law, a constitutive relation based on empirical evidence, is a statement of the fact that heat flows from hotter regions to colder regions; stated differently, heat flows down the temperature gradient. Now, we may substitute (1.20) into (1.19) to obtain a single equation for the temperature O(x, t), namely, K k= - . pC

(1.21 )

1.3.

Diffusion

17

Equation (1.21) is the heat equation; it is the diffusion equation in the context of heat flow. The constant k is called the diffusivity and is a property of the medium; values ofk for different media (metals, plastics, etc.) can be found in physical handbooks. Note that k has the same dimensions (length-squared per time) as the diffusion constant. In the sequel we shall often use u in place of e for the temperature function. 0 In some cases the thermal conductivity K in (1.20) may not be constant, but rather may depend on x if the bar is nonhomogeneous; over large temperature ranges the conductivity could also depend on the temperature e. If, for example, K = K(e) , then we obtain a nonlinear heat model

It is possible, of course, that the density and specific heat could depend on the location x or the temperature e. EXAMPLE

(Convection-Diffusion) Ifboth diffusion and convection are present, then the flux is given by

and the conservation law (1.16) becomes Ut

+ CU"

-

Du""

=

0,

which is the convection-diffusion equation. This equation would govern the density of a chemical, say, that is being convected by the bulk motion of a fluid moving at velocity c in which it is dissolved, while at the same time it is diffusing according to Fick's law. If the chemical also decays at rate A, then we include a source term, and the model is Ut

+ CU" -

Du""

=

-AU,

which is a convection-diffusion-decay equation.

o

Diffusion-like equations are often accompanied by an initial condition that specifies the density at time t = 0 at all points of the spatial domain, as well as boundary conditions that specify conditions on the density at the boundaries of the domain for all time. To fix the idea, let us consider a tube of finite length 0 ::::: x ::::: l. Then an initial condition has the form U(x, 0) = uo(x),

0::::: x ::::: 1,

where Uo is the given initial density distribution. There are three types of boundary conditions that often occur in physical problems. If the density

1. The Physical Origins of Partial Differential Equations

18

is specified at the boundary x = 0, then we have u(O, t)

= get), t >

0,

where g is given; this condition is called a Dirichlet condition. We could also specify the flux at a boundary, i.e., -Dux(O, t) = h(t),

t > 0,

which is called a Neumann condition. If h(t) == 0, then we say that the boundary is insulated; if g or h is zero, we also say that the corresponding boundary condition is homogeneous. A third type of boundary condition has the form (say, at x = 0) -Dux(O, t) = -f3(u(O, t) - 'Ijf(t)) , t > 0.

In heat flow in a bar, for example, this law expresses Newton's law of cooling, which states that the heat flux is proportional to the temperature difference between the end of the bar and the given temperature 'Ijf(t) of the environment; f3 is the constant of proportionality representing a heatloss factor. These types of conditions are also called radiation conditions, or Robin conditions. Part of the requirement of a well-posed PDE model is that the PDE be accompanied by appropriate initial and boundary conditions so that the resulting mathematical problem has a unique solution as well as a correct physical meaning.

Exercises 1. Heat flows longitudinally through a metal bar of length 10 centimeters, and the temperature u = u(x, t) satisfies the diffusion equation Ut = ku xx , where k = 0.02 square centimeters per second. Suppose the temperatures at some fixed time T atx == 4,6,8 cm are 58, 64, and 72 degrees, respectively. Estimate u xx (6, T) using a difference approximation. Will the temperature at x == 6 increase or decrease in the next instant of time? Estimate the temperature at x == 6 at T + 0.5 seconds. (Recall from calculus that 1/.

f (a) ""

f(a-h)-2f(a)+f(a+h) h2

'

where h is a small increment.) 2. Let u == u(x, t) satisfy the PDE model

o< u(O, t)

==

u(l, t)

u(x, 0)

=

uo(xJ,

==

0,

x < I, t > 0, t > 0,

0 :::: x :::: 1.

Exercises

19

Show that

Hint: Let E(t) u(x, t) if uo(x)

l

u(x, t)2dx ::::

l

uo(xidx,

t::: 0.

= f~ u(x, t)2dx and show that E'Ct)

:::: 0. What can be said about

= o?

3. Show that the problem

°

U(O, t)

= gCt),

u(x, 0)

= uo(x),

< x < I,

u(l, t)

t > 0,

= hCt),

°: : x::::

t > 0,

I,

with nonhomogeneous boundary conditions can be transformed into a problem with homogeneous boundary conditions. Hint: Introduce a new dependent variable w by subtracting from u a linear function of x that satisfies the boundary conditions at any fixed t. In the transformed problem for w, notice that the PDE picks up a source term, so you are really trading boundary conditions for source terms. 4. Show that the convection-diffusion-decay equation Ut

= Duxx -

eu" -

AU

can be transformed into the diffusion equation by a transformation of the form

= w(x, t)e -{3t. = A + c2 /(4D).

u(x, t)

Solution: Thke a

= e/(2D),

fJ

ctX

5. Heat flow in a metal rod with an internal heat source is governed by the problem Ut

u(O, t)

= ku xx + 1,

° < x < 1,

=

=

0,

u(1, t)

1,

t > 0,

t > 0.

What will be the steady-state temperature in the bar that will be reached after a long time? Does it matter that no initial condition is given? Hint: In the steady state u depends only on x. 6. A bar loses heat across its lateral boundary at a rate proportional to the temperature u. The equation is Ut

u(O, t)

°

= leuxx - au, < = 1, u(l, t) = 1,

x < 1, t > 0, t > 0.

Graph the steady-state temperature distribution in the bar and analyze how heat is flowing in the bar and through its boundaries. 7. Bacteria in a one-dimensional medium (a tube of unit cross-sectional area, length I, and capped on both ends) have a growth rate given by the logistic law ru(1 - u/K), where r is the growth constant, K is the carrying capacity, and u = u(x, t) is the density of bacteria measured in bacteria per unit

1. The Physical Origins of Partial Differential Equations

20

length. Initially, the density is given by u = ax(l- x). For t > 0 the bacteria also diffuse, with diffusion constant D. Formulate an initial boundary value problem for the density. Ifwe wait a long time, what will the density be? Use your intuition to sketch a sequence of density profiles that show how the density evolves. You may want to consider the cases aZ2 < 4K and aZ 2 > 4K separately. Hint: Think about the problem if no diffusion were present. 8. Show that the equation Ut = k(t)uxx can be transformed into the diffusion equation by changing the independent variable time tto r = f~ k( rOdr}. Show that the equation Ut = kuxx - b(t)u" can be transformed intO the diffusion equation by changing the spatial variable to ~ = x - f~ b(r})dr}.

1.4

Contaminant Transport in Aquifers

A problem of great importance in environmental science is to understand how chemical contaminants are transported through subsurface aquifer systems. In this section we derive a simple model, based on mass balance, of contaminant flow through a porous medium that incorporates the phenomena of diffusion, convection, and adsorption. We imagine that water is flowing underground through a fixed soil matrix. In any fixed volume, the fraction of space, called the pore space, available to the water is assumed to be w (0 < w < 1), which is called the porosity of the medium. In general, the porosity can vary with position, or even pressure, but for simplicity we assume that w is. constant throughout the medium. Further, we assume that the flow is saturated, which means that all of the available pore space is always filled with fluid. We will assume that the flow is one-dimensional, in the x-direction, and takes place in a tube of cross-sectional area A. Hydrogeologists often measure the discharge q of a flow, which is the actual volume of water per unit area per unit time that crosses a given section. The average velocity of the flow is defined to be v = qlw. Now, let c(x, t) be the concentration of a chemical contaminant that is dissolved in the water, and let ¢(x, t) denote its flux, or the rate per unit area that the contaminant crosses the section at x. We further assume that the chemical disappears with rate F(x, t), measured in mass per unit volume per unit time; for example, F might measure the decay rate for a radioactive contaminant, or the rate that the contaminant is adsorbed by the soil. Then mass balance for the contaminant, applied to a section a S x S b of the tube, implies the integral conservation law d dt

lb a

c(x, t)wAdx = ¢(a, t)A - ¢(b, t)A -

lb a

F(x, t)Adx.

Bringing the time derivative under the integral sign and appealing to the fundamental theorem of calculus and the arbitrariness of the section (as

1.4. Contaminant Transport in Aquifers

21

in Section 1.2), we obtain the mass balance law in differential form WCt

+ ¢x =

(l.22)

-F.

Note that this equation is, except for the porosity factor w, the same as the basic conservation law derived in Section l.2. At this point, constitutive relations must be postulated regarding the form of the flux ¢. We make the assumption that the flux contains both a Fick's law-type diffusive term and a convective term due to the bulk movement of the water. Specifically, we assume ¢

= -D'cx + qc = -D'cx + vwc.

(l.23)

The first term is basically Fick's law (Section 1.3), where D' is called the dispersion constant. The dispersion constant contains contributions from two processes, molecular diffusion and hydrodynamic dispersion. The latter is the spreading caused by the kinematic motion of the fluid through the distorted pathways in the porous matrix. The second term in the flux expresses the fact that the contaminant is transported by the groundwater, which is assumed to be flowing with average velocity v. Thus, wvc measures the actual mass of contaminant crossing the section at x, due to convection. The average velocity is the velocity that would be measured by a flow meter in the porous domain. Substituting the expression (l.23) for flux into the mass balance equation (1.22) gives WCt

= D'C

XX -

wVCx -

F,

which is a convection-diffusion (or convection-dispersion) equation with a source term. For example, for radioactive decay, F would be the decay rate given by F = Ac. A process that plays an important role in hydrogeology is adsorption. Adsorption is the process by which solutes, or chemicals dissolved in the water, attach themselves to the solid soil matrix. Exactly how a particle binds itself to the soil can be a complicated physical process, and the mechanics of this process are discussed in books on hydrogeological phenomena. 'Ib measure adsorption, we introduce a quantity s = sex, t), which is the mass of the solute, per unit mass of the solid, that is adsorbed onto the soil. Then we write F, the sorption rate, as F = (1 - w)pSt,

where P is the density of the soil. Note that the factor 1-w is the fraction of the medium that is soil. Then, upon dividing by w, the governing equation becomes 1-w (l.24) Ct = Dcxx - vCx - - - PSt, W

which is the basic transport equation studied in groundwater theory. Here, D = D'/w. If there is no flow (v = 0) and no adsorption (St = 0),

1. The Physical Origins of Partial Differential Equations

22

then in the context of solute transport, the resulting diffusion equation = Dcxx is frequently called the dispersion equation. Equation (1.24) still has two unknowns, the solute concentration c and the sorbed concentration s. Clearly, there ought to be some relation between these two quantities. Intuitively, as the solute concentration increases, the sorbed concentration should increase; but, as more and more is adsorbed, the attachment process should slow down as a maximum saturation is reached on the soil particles. There are several experimental constitutive relations that model this type of behavior, based on an instantaneous relation between c and s of the form s = f(c). Such relations are called equilibrium isothenns because they are measured at constant temperature and at equilibrium. Three of these are Ct

s

= klc -

k2C 2 ,

S

= kcl/n,

n:::

quadratic isotherm, I,

klc s- - - - 1 + k2C'

Freundlich isotherm, Langmuir isotherm.

One type of solution that is interesting in hydrogeological applications is called a wave-front solution. This is a traveling wave solution of the form c = C(x - at) that has constant limits at ±oo. It represents a concentration wave moving into a constant state at speed a. Let us illustrate how to find wave fronts by considering the quadratic isotherm with k2 < k l . In this case, equation (1.24) can be written (1.25) =' Cxx - cx , where, for simplicity, we have taken D = v = land b == kl Ps(l- w)l wand

((1

+ b)c + mc2 )t

m == k 2 Ps(l - w)/w. We now assume solutions of the form c = C(z) , z = x - at with C( -00) = land C(oo) = O. The chain rule requires that t derivatives transform via

Moreover,

Cx

=

a d at = -a dz' C'(z) and Cxx = C"(z). Therefore, (1.25) becomes -a((l + b)C + mC 2 )' = C" - C'.

This is an ordinary differential equation for the shape of the waveform C(z). Ifwe integrate with respect to z from an arbitrary value of z to z = +00, we obtain, usingthefactthatC(+oo) = oand assuming C'(±oo) = 0, -a((l

+ b)C -

mC 2 )

= C' -

C,

which can be arranged to yield C' + (a + ab - l)C

= amC2 •

Exercises

23

0.8

Figure 1.5. The wave front C(z) given by (1.26) in the fixed frame z. This profile moves to the right at speed a in the x-coordinate frame.

0.2

-4

-2

This is a Bernoulli equation that can be solved exactly for C(z). But first let us observe that if we take the limit as z --+ 00, we obtain

a + ab -1 = am,

1

a =

or

1 +b-m

.

Therefore, we have determined the wave speed a, which is positive by our assumption that kl > k z. Now, the Bernoulli equation can be transformed to a linear equation by the substitution W = C- 1 • We get

w' - amW

= -am,

which easily is solved to obtain W(z)

= 1 + Ke amz ,

where K is a constant of integration. Thus, choosing K = 1 (notice that K only shifts the wave profile but does not change its shape), we get C(z)

=

1 1

+ eamz

.

(1.26)

Therefore we have determined a wave-front solution of the form c(x,t)

1

= 1 + emn(x-at).

The function C(z) is shown in Figure 1.5.

Exercises 1. Sketch graphs of the quadratic, Freundlich, and Langmuir isotherms and discuss their qualitative similarities and differences. 2. Show that the pure dispersion equation Ct = kc xx does not admit nonconstant wave-front solutions.

1.

24

The Physical Origins of Partial Differential Equations

3. Find a formula for the wave-front solution u equation Ut

+ UUx

=

= U(z), z = x -

ct, to Burgers'

Uxx

that satisfies the conditions U( -00) = 1, U(oo) = O. Choose U(O) = 0.5. 4. Consider the pure dispersion equation Ct = kcxx for x, t > o. (a) Try to find solutions of the form c(x, t) = fez), where z = xl Nt. Find an ordinary differential equation for fez) and solve it. (b) What conditions would fez) have to satisfy if c satisfies c(O, t) = 1 for t > 0 and limr.....o+ c(x, t) = 0 for x > O? Use these two conditions to determine the arbitrary constants in your solution to part (a). Interpret this problem in the context of solute transport. 5. Consider the transport equation (1.24) with a linear isotherm s = Kc (the rate constant K is called the distnbution coefficient). Show that one obtains the governing equation RCt

= Dcxx -

vCx,

where R == 1 + pK(l - w)lw is the retardation constant. Does this equation admit wave-front solutions? Compare this equation to the transport equation with no sorption and indicate why R is called the retardation constant.

Vibrations of a String

1.5

Wave motion is one of the most commonly occurring phenomena in nature. We mention electromagnetic waves, water waves, sound waves, and stress waves in solids as typical examples. In this section we select a simple model to introduce the concept of wave motion and one of the fundamental PDEs, the wave equation, that describe such phenomena. The example is the small, transverse vibrations of a flexible, elastic string, e.g., a guitar string. Let us imagine a taut string oflength 1 fastened at its ends. Its motion is described by a function u = u(x, t) that gives the vertical displacement T(b, t) 8(b, t)

1 1 1

1 1 u(a, 1 1

u(b, t)

t)

I

a

b

Figure 1.6. Displaced string with tension forces shown.

1.5.

Vibrations of a String

25

of each point x of the string at time t. Our basic postulate is that the displacement u is small. Implicitly, we assume that the motion is in a plane and no element of the string moves horizontally-only vertically. At each instant of time we assume that the string has a density p(x, t), with dimensions mass per unit length, and the tension in the string is given by a function T(x, t), with force dimensions. By convention, T(x, t) is the force on the segment to the left of x caused by the portion of the string to the right of x, and we assume that the tension always is directed along the tangent to the profile at x. This latter assumption implies that the string does not resist bending. A priori, we do not know p or T. Figure 1.6 shows an arbitrary segment of the string between x = a and x = b. We denote the angle that the tangent makes with the horizontal by O(x, t), and we observe that tan O(x, t) = ux(x, t). 1b obtain an equation of motion of the string, we appeal to mass balance and Newton's second law. First, mass balance implies that the mass of the segment at any time t must equal its mass in equilibrium position (which we take to be at t = 0). In symbols, this means

where Po is the density of the string in equilibrium. (Note that

is an element of arc length.) But this equation must hold for every segment [a, b], and therefore we may equate the integrands to obtain p(x,

t)J

1 + ux(x, t)2 = Po(x).

(1.27)

Next let us apply Newton's second law, which states that the time rate of change of the total momentum of the segment must equal the net external force. We shall assume that the only force acting on the segment is tension caused by the rest of the string; we ignore gravity and any damping forces-see Exercises 1 and 2. Because there is no horizontal motion of the segment, the horizontal forces must balance, or T(b, t) cos O(b; t) - T( a, t) cos O(a, t) = O.

Because this equation must hold for all a and b, we must have T(x, t) cos O(x, t) = ret)

26

1.

The Physical Origins of Partial Differential Equations

for some function r. Now, the time rate of change of the total momentum of the segment must equal the net vertical force. That is,

!

lb p(x, t)Ut( x, t)J 1 + ux(x, t)2dx

= T(b, t) sin B(b, t) -

T( a, t) sin B(a, t).

Using the result of equation (1.27) and bringing the derivative inside the integral sign gives lb Po(x)Utt(x, t)dx

= T(b, t) sin B(b, t) -

T(a, t) sin B(a, t).

We observe that the right side of this equation can be written as T(b, t) cos B(b, t)[tan B(b, t) - tan B(a, t)],

or r(t)( ux(b, t) - u x( a, t)).

Therefore, lb Po(X)Utt(X, t)dx

= r(t)(ux(b, t) -

ux(a, t)).

By the fundamental theorem of calculus we can rewrite the right side as an integral, and we obtain b lb Po(x)Utt(x, t)dx = r(tyl uxx(x, t)dx .

.Again using the arbitrariness ofthe interval of integration (what we have said holds for any segment [a, b]), we can deduce Po(x)Utt(x, t) = r(t)uxx(x, t),

0 < x < 1, t > O.

As an additional constitutive assumption, we assume that the tension r(t) = ro = constant, which is a good assumption for small displacements of the string. Thus, in the framework of the assumptions that we have made, we obtain the final model equation for the transverse displacements, namely, (1.28) If we introduce the function co(x) = .)rol Po(x),

then we can write (1.28) as Utt =

Co (xi Uxx ·

(1.29)

One should observe that Co has dimensions of speed. It is called the wave speed, and we shall observe thatit measures the actual speed that waves

Exercises

27

travel along the string. In many applications the density of the string is a constant, and hence Co is constant; in that case equation (1.29) becomes (1.30) which is called the wave equation; it is one of the fundamental equations in mathematics and applications. If Co is constant, it is easy to check (Exercise 3) that u = F(x - cot) and u = G(x + cot) are solutions of the wave equation for any twice continuously differentiable functions F and G; thus the wave equation admits both right- and left-traveling wave solutions. In fact, we shall show later that the general solution ofthe wave equation is the superposition (sum) ofleft- and right-traveling waves. Our original assumption that the string is fastened at its ends x = 0 and x = 1 translates into boundary conditions u(O, t)

=

0,

=

u(l, t)

0,

t 2: O.

(1.31 )

To determine a unique displacement, we must impose initial conditions as well; that is, we prescribe both the position of the string and its velocity at time t = O. Symbolically, u(x, 0)

= f(x) ,

Ut(x, 0)

= g(x) ,

0:::: x :::: 1.

(1.32)

In summary, to determine the displacement u(x, t), we must solve the initial boundary value problem consisting of the PDE (1.30) subject to the boundary conditions (1.31) and initial conditions (1.32). We point out that this model is a small-displacement theory. A more complicated nonlinear model is required if we want to describe large displacements.

Exercises 1. In the derivation of the wave equation (1.30) we assumed no gravitational force on the string. What would the model equation be if gravity were included? Note that gravity acts at each point of the string and derive the equation Utt

= c~uxx -

g,

where g is the constant acceleration due to gravity. 2. Repeat the derivation in this section when, in addition, the vertical motion is retarded by a damping force proportional to the velocity of the string. Obtain the damped wave equation Utt

= co(xiu.yx

-

ku l ,

where k is the constant damping coefficient.

28

1.

The Physical Origins of Partial Differential Equations

3. Let Co be constant. Verify that u = F(x - cot) and u = G(x + cot) are solutions of (1.30) for any twice-differentiable functions F and G. If F(x) = 11(1 + x 2 ), sketch time profiles of the wave u(x, t) =

1

"2 (F(x

- t)

+ F(x + t)).

4. Consider the displacements of a string governed by the wave equation

where

C

is constant, and subject to the boundary conditions

=

u(O, t)

=

u(l, t)

0,

0,

t::: o.

Show that for any positive integer n there is a solution of the form mrct mrx un(x, t) = cos -1- sin -1- .

What are the initial conditions? Describe the motion ofthe string in the cases n = 1 and n = 2 (sketch several fixed time snapshots of the string profile). In general, for any n, what is the temporal frequency of the oscillations represented by these solutions? These frequencies are called the fundamental frequencies, and the corresponding displacements are called the fundamental modes. How do the frequencies change as the length 1 of the string is changed? How do the frequencies depend on the tension? Discuss your answers in the context of a vibrating guitar string. These special solutions are called standing waves. 5. The total energy of the string governed by equation (1.28) with boundary conditions (1.31) is defined by E(t) =

1

1 2 ( - POUt

1

0

1

2

+ - roux)dx.

22

Show that the total energy is constant for all t ::: O. Hint: Carry out this calculation by multiplying (1.28) by Ut and noting that (u;)t = 2ututt and (UtUx)x = UtU xx + UtxUx. Then derive the formula

dt

2

1

dt Jo POUt dx = rOUtUx 10 -

1.6

dt Jo rouxdx.

dt

2

Quantum Mechanics

One of the fundamental equations of mathematical physics is the Schrbdinger equation, which is the basic equation of quantum mechanics. In the next few paragraphs we give a brief description of Schrbdinger's equation in one dimension. Let us consider a particle of mass m moving on the x-axis under the influence of a continuous, conservative force given by F(x) that depends only on the position x. According to the canon of classical particle mechanics, the motion x = xU) of the particle is governed by the dynamical

1.6. Quantum Mechanics

29

equation (1.33) which is Newton's second law of motion (mass times acceleration equals force). If the initial position x(O) and initial velocity dxl dt at t = 0 are known, then one can in theory solve the ODE model (1.33) subject to the initial conditions and actually determine the state of the particle, that is, its position and velocity, for all times t > O. In this sense, classical mechanics is a deterministic model-knowledge of the initial state determines the future history of the particle. In the early 1900s it was found that this classical, deterministic model of motion fails on the atomic scale. Quantum mechanics, which is a probabilistic model, grew out of the atomic physics revolution. Quantum theory dictates that the particle has no definite position or velocity; rather, one postulates a statistical, or probabilistic, interpretation of the state of the particle in terms of a wave function 'II(x, t). The square of the wave function,l'II1 2, is a probability density; that is,

lb

1'II(x, t)1 2 dx

is interpreted as the probability of the particle being in the interval a :s x :s b at time t. Thus J~oo 1'II(x, t)1 2 dx = 1, since the particle is located somewhere on the x-axis. From elementary probability theory it is known that the probability density 1'II(x, t)12 coritains all of the statistical information for a given problem (for example, the mean and variance of the position). Note that'll is complex-valued, and thus 1'1112 = 'II'll, where the overbar denotes complex conjugation. So the question is how to find the wave function. The equation that governs the evolution of a quantum-mechanical system (the analogue of (1.33) for a classical system) is the Schrodinger equation, a second-order partial differential equation having the form Fz2 in'llt = - - 'IIxx 2m

+ V(x) 'II,

x

E

R, t > 0,

(1.34)

where V = Vex) is the potential energy, m is the mass, andFz = hl(2rr), where h = 6.625 . 10-34 kg m 2Is is Planck's constant. (Recall that associated with the force F(x) is a potential function Vex) defined by the equation F(x) = - V'(x); that is, the force is the negative gradient of the potential.) One can motivate the Schrodinger equation from momentum and energy conSiderations, but here our goal is only to note the form of the equation. The reader is referred to one of the many excellent treatises on quantum mechanics for a complete discussion. A popular equation studied in detail in the mathematical literature is the free Schrodinger equation (1.35)

1. The Physical Origins of Partial Differential Equations

30

where the potential V is taken to be zero (hence a free particle with no forces acting); the constants are taken to be unity. Note that this equation resembles the diffusion equation, but it has a complex coefficient; this makes solutions of the two equations quite different. One method used to find solutions of (1.34) (we shall observe later that this method is basic in PDEs) is to assume that the variables separate, i.e., the solution has the form of a product W(x, t) = y(x)¢(t). Substituting this into (1.34) gives, after rearrangement, iJi¢'(t) ¢(t)

=

- frny"(x)

+ V(x)y(x)

y(x)

Here, prime denotes the derivative. The left side of this equation depends only on t, and the right side depends only on x. The only way equality can occur for all t and x is if both sides are equal to the same constant, which we shall call E. Therefore, we obtain two equations, one for ¢,

d:

(1.36)

= (-iEIn)¢,

and one for y, namely,

n2

- - y" 2m

+ (V(x) -

E)y

= O.

(1.37)

The solution of (1.36) is easily found to be ¢

= C exp( -iEtln) = C(cos( -Et/Ji) -

i sin(Et/Ji)),

where C is a constant. Thus, the temporal part is oscillatory (periodic). Equation (1.37), whose solution y(x) gives the spatial part of the wave function, is called the time-independent Schrodinger equation, and it is one of the fundamental ODE models of mathematical physics. The values of E for which (1.37) has a nontrivial solution y(x) with f~oo y(xidx = 1 are interpreted as the allowable energy levels. This is in contrast to the situation in classical mechanics, where any energy level is possible.

Exercises 1. Show that the ordinary differential equation y" + p(x)y' + q(x)y = 0 can be transformed into the Schrodinger-like equation u" + r(x)u = 0, r = q p'/2 - p2 / 4, without a first derivative, by the Liouville transformation u = y exp(! p(~)~).

J:

2. Show that J~oo 1\II(x, t)1 2 dx = J~oo y(xldx. 3. If Vex) = ~ x 2 , then for a special value of E, equation (1.37) admits a solution of the form y(x) = Ce- ax2 for some a > 0 and any constant C. Determine

3l

1.7. Heat Flow in Three Dimensions

a and E, and sketch a graph of the probability density function 1\lII2 in the case k = m = 1i (use a calculator to compute C). Numerically calculate the probability that the particle is in the interval [0, 1]. 4. Suppose \lI(x, t), x E R, t > 0, is a twice continuously differentiable function that satisfies (1.35), and both 1\lII, 1\lI,,1 ~ 0 as Ixl ~ 00. Show that

i:

1\lI(x, t)1 2 dx

=

constant, t > O.

Hint: Thke the time derivative of the left side, pulling it under the integral; use the Schrodinger equation and its complex conjugate, and finally integrate the terms like ii1",,\lI by parts. 5. A free particle of mass m is confined to the interval 0 < x < 1!, and \lI(O, t) = \lI(1!, t) = 0 for all t. Show that the associated time-independent problem is y"

2mE

+7

y

= 0,

0 < x
0, B = 0 (and so D < 0) graphs as an ellipse in the xt-plane. Similarly, Ax2 - Ct 2 = 1 graphs as a hyperbola, and so forth. As it turns out, all parabolic equations are diffusion-like, all hyperbolic equations are wave-like, and all elliptic equations are static. Later we shall show why this is true. If A, B, and C are not constants, but rather functions ofthe independent variables x and t, then the discriminant D depends on x and t. We make the same definitions in this case, but now the sign of D can change, depending upon the domain. Problems that change type from one region to another can be difficult, and we do not consider them in this text. Classification of systems of PDEs is also possible, as well as classification of PDEs with

46

1. The Physical Origins of Partial Differential Equations

nonlinear principal parts where the coefficients A, B, and C depend on the solution U itself. The reader is referred to a more advanced treatment for a discussion of these matters. Now we show that the principal part Lu in equation (1.49) can be simplified for each of the three types by introducing a new set of independent variables. This strategy of searching for a change of variables that simplifes a problem is common in differential equations. We examine the case where A, B, and C are constant, and we seek a linear transformation ~

= ax + bt,

i

= CX + dt

that simplifies Lu. Here ~ and i are new independent variables, and a, b, c, and d are to be determined. We assume that ad - bc =J 0, so that the transformation is invertible, that is, we may solve for x and t in terms of ~ and i. The dependent function u in the new variables will be denoted by U = U(~, i); that is, u(x, t) = U(ax + bt, cx + dt). Then, by the chain rule for derivatives,

= UI;~x + U,ix = aUI; + cU" Ut = UI;~t + U,it = bUI; + dU,.

Ux

1b get the second partials requires another application of the chain rule. We have Uxx

=

a

ax (aUI;

+ cU,)

+ cUI;,) + c(aUI;, + cU,,) a 2 UI;I; + ZacUI;, + c2 U".

= a(aUI;I; =

The other two second partials are computed similarly, and we get Utt Uxt

= b2 UI;I; + ZbdUI;, + d 2 U", = abUl;1; + (ad + cb)UI;, + cdU".

Substituting these quantities into the principal part and collecting terms, we get Auxx

+ BUxt + CUtt

+ Bab + Cb 2 )UI;I; + (ZacA + B(ad + bc) + ZCbd)UI;, + (Ac2 + Bcd + Cd 2 )U".

= (Aa 2

(1.50)

Now we can select a, b, c, d so that some of the second partials in the new variables disappear. This process must be handled differently depending on the sign of the discriminant D.

Hyperbolic case, D > 0 We have some flexibility, so let lis choose a = c = 1. Then the coefficients of UI;I; and U,,' which have the same form, become quadratic expressions

1.10.

Classification of PDEs

47

in band d, respectively, and we can force those coefficients to vanish by choosing

-B+JD

b=

-B- JD .

d=

2C'

(1.51 )

2C

Here we have just used the quadratic formula. The remaining coefficient, that of U~r, is nonzero (we leave this as an exercise). In summary, the transformation ~

= x + ( -B+JD) t, 2C

•

= x + (-B-JD) 2C

t

(1.52)

transforms the PDE (1.49) into a simpler equation of the form U~r

+ G(~, ., U,

U~,

Ur) = 0,

where only the mixed second partial derivative appears. Thus there is a significant simplification over (1.49), where all the second partial derivatives occur. This latter equation is called the canonical form of a hyperbolic equation, and the coordinates ~ and • defined by (1.52) are called the characteristic coordinates. Transformation to these coordinates is almost always a preferred strategy for hyperbolic PDEs. Finally, if C = 0, then (1.51) is not valid; in this case select b = d = 1 and

a=

-B+JD

c=

2A

-B-JD 2A

Parabolic case, D = 0 Now equations (1.51) give b = d, and the resulting transformation ~ = x + bt, r = x + bt is not invertible. So we must proceed differently. Observe that if we choose a = c = I, d = -BI(2C), and b = 0, then now the coefficients of UTI and U~r in (1.50) vanish (another exercise!). Therefore, the transformation ~

= x,

B • = x - -t 2C

transforms equation (1.49) into U~~

+ H(~, ., U,

U~,

Ur) = 0,

where only a double partial appears, and this is the canonical form in the parabolic case. Observe its Similarity to the heat equation.

Elliptic case, D < 0 Now band d in (1.51) are conjugate complex numbers, i.e., d Selecting a = c = I, we obtain a complex transformation ~ =

x

+ bt,

• = x

+ bt.

b.

1.

48

The Physical Origins of Partial Differential Equations

But then, a real transformation can be found by taking real variables 1

= 2: (~ + r), fJ =

cx

1 2i (~ - r).

It is an exercise to show that this transforms (1.49) into the equation

Uaa

+ Uf3f3 + K(cx, fJ,

U, Ua , U(3)

=

0,

where both double partials are present and the mixed partial is absent; the last equation is the canonical form for elliptic equations. We recognize the combination of second partials as the Laplacian operator. Generally, characteristic coordinates are useful only for hyperbolic equations, and they do not playa particularly important role in elliptic and parabolic equations.

Exercises 1. Classify the PDE

+ 2kux t + kZUtt

Uxx

= 0,

k

f.

O.

Find a transformation ~ = x + bt, r = x + dt of the independent variables that transforms the equation into a simpler equation of the form U~~ = O. Find the solution to the given equation in terms of two arbitrary functions. 2. Find the solution of the PDE 2u."" - 4u x t

+ Ux

= 0

in terms of two arbitrary functions. Hint: Make a transformation. 3. Classify the PDE XU xx -

4u x t = 0

in the region x > O. Observe that the PDE has variable coefficients. Solve the equation by making the nonlinear change of variables ~ = t, r = t + 4ln x. The solution is U

= e- t/ 4[Ct + 4ln x) + get),

where [ and g are arbitrary functions. 4. Show that the equation Utt -

cZuxx

+ aUt + bux + du

= [(x, t)

can be transformed into an equation of the form W~r

+ kw

=

get r),

by first making the transformation ~ U = we ClHfJr for some choice of a, {J.

=

w = w(~, r)

x - ct, r

=

x

+ ct and then letting

Exercises

49

5. ClassifY the PDE U xx

- 6uxy

+ 12uyy

= O.

Find a transformation of independent variables that changes it to Laplace's equation.

Partial Differential

Equations on

CHAPTER

2.1

Unbounded Domains

Cauchy Problem for the Heat Equation

In this chapter we investigate problems on unbounded spatial domains. In some sense problems on infinite domains are easier than problems that involve boundaries. We begin with the heat equation, or diffusion equation, on the real line. That is, we consider the initial value problem x

E

R, t > 0,

(2.1)

U(x, 0) = (x) , x

E

R.

(2.2)

Ut

= ku xx ,

Physically, this problem is a model of heat flow in an infinitely long bar where the initial temperature (x) is prescribed. Notice that there are no boundaries in the problem, so we do not prescribe boundary conditions explicitly. However, for problems on infinite domains, conditions at infinityare sometimes either stated explicitly or understood. Such a condition might require boundedness of the solution or some type of decay of the solution to zero as x ~ ±oo. In mathematics, a pure initial value problem like (2.1 )-(2.2) is often called a Cauchy problem. Deriving the solution of (2.1)-(2.2) is accomplished in two steps. First we will solve the problem for a special step function (x), and then we will construct the solution to (2.1)-(2.2) using that special solution. So, first let us consider the problem Wt

= kwxx , x E R, t > 0,

50 J. D. Logan, Applied Partial Differential Equations © Springer-Verlag New York, Inc. 1998

(2.3)

2.1.

Cauchy Problem for the Heat Equation

w(x, 0)

=0

for

X

< 0;

51

w(x, 0)

= Uo

for x > 0,

(2.4)

where we have taken the initial condition to be a step function with jump Uo·

We motivate our approach to the solution of (2.3)-(2.4) with a simple idea from the subject of dimensional analysis. Dimensional analysis deals with the study of units (seconds, meters, kilograms, and so forth) and dimensions (time, length, mass, and so forth) of the quantities in a problem and how they relate to each other. Equations must be dimensionally consistent (one cannot add apples to oranges), and important conclusions can be drawn from this fact. The cornerstone result in dimensional analysis is called the pi theorem. The pi theorem guarantees that whenever there is a physical law relating dimensioned quantities Q1, ... , qm, then there is an equivalent physical law relating the independent dimensionless quantities that can be formed from qr, ... ,qm. By a dimensionless quantity we mean one in which all the dimensions (length, time, mass, etc.) cancel out. As a simple example take the law h

=

1 -zgt

2

+vt

that gives the height h of an object at time t when it is thrown upward with initial velocity v; the constant g is the acceleration due to gravity. Here the dimensioned quantities are h, t, v, and g, having dimensions length, time, length per time, and length per time-squared. This law can be rearranged and written equivalently as

!::=-~(gt)+l vt 2 v in terms of the two dimensionless quantities ][1

h vt

== -

and

][2

== gt V

For example, h is a length and vt, a velocity times a time, is also a length; so ][1, or h divided by vt, has no dimensions. Similarly, ][2 = gtlv is dimensionless. A law in dimensioned variables can always be reformulated in dimensionless quantities. So the physical law can be written as ][1

= - ~][2

+ l.

We use similar reasoning to guess the form of the solution of the initial value problem (2.3)-(2.4). First we list all the variables and constants in the problem: x, t, w, uo, k. These have dimensions length, time, degrees, degrees, and length-squared per time, respectively. We notice that wluo is a dimensionless quantity (degrees divided by degrees); the only other independent dimensionless quantity in the problem is xl.../4kt (the "4" is included for convenience). We expect, therefore, that the solution can be

2. Partial Differential Equations on Unbounded Domains

52

written as some combination of these dimensionless variables, or W x Uo = f( ../4kt) for some function f to be determined. In fact, this is the case. So let us substitute w=f(z),

X

Z=

~

v4kt

into the PDE (2.3). We have taken Uo = 1 for simplicity. The chain rule allows us to compute the partial derivatives as

,Ix

=f

Wx

I = f, (z)zx =f ~'(z),

Wxx

(z)Zt

=- -

,

Wt

r::17>f (z),

2 v4kt3

v4kt

a ax

1" = -wx = -k f (z). 4 t

Substituting these quantities into (2.3) then gives, after some cancellation, an ordinary differential equation, f"(Z)

+ 2zf'(z) = 0,

for fez). This equation is easily solved by multiplying through by the integrating factor e 2 and integrating to get . f'(z) = cle-z2 ,

where

Cl

is a constant of integration. Integrating from 0 to z gives fez)

l

= Cl

z

e- r2 dr

+ C2,

where C2 is another constant of integration. Therefore we have determined solutions of (2.3) of the form w(x, t) =

t/,J4'fi

Cl

10

e- r2 dr

+ C2·

Next we apply the initial condition (2.4) (with Uo = 1) to determine the constants Cl and C2. For a fixed x < 0 we take the limit as t ~ 0 to get

o = w(x, 0) = Cl

11

For a fixed x > 0 we take the limit as t 1 = w(x, 0) =

Cl

00

e- r2 dr

~

00

+ C2.

0 to get

e- r2 dr

+ C2.

2.1.

Cauchy Problem for the Heat Equation

53

Recalling that

100

00

../ii 2 '

-r2d

e

r=-

we can solve the last two equations to get Cl Therefore, the solution to (2.3)-(2.4) with Uo = 1 is

w(x, t)

= -1 + 2

1

r.;;

yJr

l/../ii,

Cz

1/2.

1o"/.J4ki e_r2 dr.

(2.5)

(~))

(2.6)

0

This solution can be written nicely as

w(x, t)

=

~

(1 +

erf

in terms of a special function called the "erf" function, which is defined by erf (z)

=

2

t

2

../ii 10 e- r dr.

Figure 2.1 shows a graph of several time snapshots ofthe solution (2.6). Now we will use (2.5) and a physical argument to deduce a solution to the Cauchy problem (2.1)-(2.2). Later, in Section 2.7, we present an analytical argument based on Fourier transforms. We make some observations. First, if a function w satisfies the heat equation, then so does w", the partial derivative of that function with respect to x. This is easy to see because

o=

(Wt - kw""),, = (w,,)t - k(w")",,.

Therefore, since w(x, t) solves the heat equation, the function G(x, t)

-8

-6

-4

-2

==

w,,(x, t)

Figure 2.1. Thmperature profiles for different times when the initial temperature is a step function.

54

2. Partial Differential Equations on Unbounded Dommns

solves the heat equation. By direct differentiation we find that G(x, t) = _1_ e-x2/ (4kt)

.J4nkt

(2.7)

The function G is called the heat kernel or fundamental solution to the heat equation; the reader will note that for each t > 0 it graphs as a bell-shaped curve (see Exercise 1, Section 1.1), and the area under the curve for each t > 0 is one; that is,

i:

G(x, t)dx = 1,

t > O.

G(x, t) is the temperature surface that results from an initial unit heat source, i.e., injecting a unit amount of heat at x = 0 at time t = O. We further observe that shifting the temperature profile again leads to a solution to the heat equation. Thus, G(x - y, t), which is the temperature surface caused by an initial unit heat source at y, solves the heat equation for any fixed, but arbitrary, y. If ¢(y), rather than unity, is the magnitude of the source at y, then ¢(y)G(x - y, t) gives the resulting temperature surface; the area under a temperature profile is now ¢(y), where y is the location ofthe source. Now, let us regard the initial temperature function ¢ in (2.2) as a continuous distribution of sources ¢(y) for each y E R. Then, superimposing all the effects ¢(y)G(x - y, t) for all y gives the total effect of all these isolated sources; that is, u(x, t) = =

i:

1

00

-00

¢(y)G(x - y, t)dy 1 ¢(y) _ _ e-(x-Yi l(4kt)dy

.J4nkt

(2.8)

is a solution to the Cauchy problem (2.1 )-(2.2) for reasonable assumptions on the initial condition ¢. More preCisely, if ¢ is a bounded continuous function, then u(x, t) given by (2.8) is a solution to the heat equation and u(x, t) -+ ¢(x) as t -+ 0+. If ¢ is only piecewise continuous (i.e., it has finitely many jump discontinuities in any bounded interval [a, b] in R), then u(x, t) still solves (2.1), but as t -+ 0+ the solution approaches the average value of the left and right limits at a point of discontinuity of ¢; in symbols, u(x, t) -+ (¢(x-) + ¢(x+)) as t -+ 0+. (A precise definition of piecewise continuity is given in Section 3.3). This discussion of the Cauchy problem for the heat equation has been intuitive, and it is a basis for understanding why the solution has the form it does. There is another way to write the solution (2.8) to the Cauchy problem (2.1 )-(2.2). If we change variables in the integral using the substitution

!

Exercises

55

r = (x - y)!.J4ki, then dr = -dy/.J4ki, and (2.8) becomes u(x, t) =

1r:;;

y Tr

1

00

(2.9)

e- r 2 ¢(x - r.J4kt)dr.

-00

This formula is called the Poisson integral representation. We make several observations. First, the solution (2.8) or (2.9) of the Cauchy problem is an integral representation; although the formula is not complicated, for most initial conditions ¢(x) the integration cannot be performed analytically. Therefore, numerical or computer evaluation of the integral will ultimately be required if temperature profiles are desired. Also, notice that the temperature u(x, t) is nonzero for every real x, even if ¢ is zero outside a small interval about the origin. Thus, a signal propagated by the heat, or diffusion, equation travels infinitely fast; according to this model, if odors diffuse, a bear would instantly smell a newly opened can of tuna ten miles away. Next, although we will not prove it, the solution given by (2.8) is very smooth; that is, u in infinitely differentiable in both x and t in the domain t > 0; this is true even if ¢ is piecewise continuous. Initial signals propagated by the heat equation are immediately smoothed out. Finally, we note that the heat kernel G(x, t) defined in (2.7) is also called the Green's function for the Cauchy problem. In general, the Green's function for a problem is the response of a system, or the effect, caused by a point source. In heat flow on the real line, G(x, t) is the response, i.e., the temperature surface caused by a unit, point heat source given to the system atx = 0, t = O. Some ofthe references discuss the construction of Green's functions for a variety of problems. Because of the basic role this function plays in diffusion problems, G(x, t) is also called the fundamental solution to the heat equation.

Exercises 1. Solve the Cauchy problem (2.1 )-(2.2) for the following initial conditions. (a) O.

Hint: Change variables as in the derivation of Poisson's integral representation.

2.2

Cauchy Problem for the Wave Equation

The one-dimensional wave equation is Utt -

c2 uxx = O.

(2.10)

We observed in Section 1.5 that it models the amplitude of waves on a string. It also arises in acoustics, in electromagnetic wave propagation, and in the mechanical vibrations of elastic media, as well as in other problems. It is a hyperbolic equation and is one of the three fundamental equations in PDEs (along with the diffusion equation and Laplace's equation). Under the transformation of variables (to characteristic coordinates) ~

= x - ct,

T

=

X

+ ct,

the wave equation is transformed into the canonical form U,~

= 0,

U

= U(~, T),

which can be integrated twice to obtain the general solution U(~,

T}=

F(~)

+ G(T),

2.2. Cauchy Problem for the Wave Equation

57

where F and G are arbitrary functions. Thus, the general solution to (2.10) is u(x, t)

= F(x -

ct)

+ G(x + ct).

(2.11)

Hence, solutions of the wave equation are the superposition (sum) of right- and left-traveling waves moving at speed c. The Cauchy problem for the wave equation is = 0,

cZuxx

Utt -

= fcx),

u(x,O)

R, t > 0,

X E

Ut(x, O)

(2.12)

= g(x),

X E

R.

(2.13)

Here, f defines the initial displacement of an infinite string, and g defines its initial velocity. The equation is second-order in t, so both the position and velocity must be specified initially. There is a simple analytical formula for the solution to the Cauchy problem (2.12)-(2.13). It is called d'Alembert's fonnula, and it is given by u(x, t) = -1 CfCx

2

-

ct)

+ fcx + ct)) + -1

2c

l

x+ct

x-ct

g(s)ds.

(2.14)

If f" and g' are continuous, then it is a straightforward exercise in differential calculus, using Leibniz's formula, to verify that this formula solves (2.12)-(2.13). The formula can be derived (see Exercise 1) by determining the two functions F and G in (2.11) using the initial data (2.13). Much insight into the behavior of solutions comes from examining the special case where the initial velocity is zero and the initial displacement is a bell-shaped curve. Specifically, we consider the problem (with c = 2) Utt -

4uxx = 0,

u(x, 0)

X E

= e_x , 2

R, t > 0,

Ut(x,O)

= 0,

X E

R.

The exact solution is, by d'Alembert's formula, u(x, t) =

~ (e-Cx-Zti + e-CX+2t)\ 2

Either the solution surface or wave profiles can be graphed easily using a computer algebra package. Figure 2.2 shows the solution surface; observe how the initial signal splits into two smaller signals, and those travel off in opposite directions at speed c = 2. In the exercises the reader is asked to examine the case where f = 0 and g ::j:. 0; this is the case where the initial displacement is zero and the string is given an initial velocity, or impulse, by, say, striking the string with an object. Close examination of d'Alembert's formula reveals a fundamental property of the wave equation. If the initial disturbance is supported in some interval a ::: x ::: b (this means that it is zero outside that interval, so the signal is located only in a ::: x ::: b), then the signal is always zero outside the region bounded by the two straight lines x + ct = a and x - ct = b.

58

2.

Partial Differential Equations on Unbounded Domains

Figure 2.2. The ~0lution surface.

See Figure 2.3. This region is called the region of influence of the interval [a , b]. An initial signal in [a, b] can never affect the solution outside this region. The lines x + ct = constant are paths in space-time along which signals are propagated at velocity -c; the lines x - ct = constant are paths in space-time along which signals are propagated with velocity c. These two families of lines are called the negative and positive characteristics, respectively. If the interval [a, b] is shrunk to a point, then the region of influence shrinks to a cone, which is called the light cone. Looking at the situation in reverse, we can ask what initial data can affect the solution at a point (xo , to) . From the d'Alembert formula, only the initial values in the interval [xo - cto, Xo + cto] will affect the solution at (xo, to). This interval is called the domain of dependence, and it is found by tracing the characteristics emanating from the point (xo, to) backward in time to the x-axis. In summary, there are important points to note regarding the characteristic curves. First, they are curves that carry the signals forward in space-time with velocity c and -c. Second, they define a special coordinate system ~ = x - ct, r = x + ct under which the wave equation Utt - c2 u xx = 0 is reduced to the simple canonical form u;r = O. In hy(Xo. to)

x domain of dependence

Figure 2.3. Region of influence and domain of dependence.

Exercises

59

perbolic problems there is always a set of characteristic curves that play these roles. Even first order PDEs, which are actually wave-like, have one family of such curves that carry signals and provide a distinguished coordinate sytem where the problem simplifies (recall the examples in Section 1.2). Finally, we point out again the important differences between parabolic and hyperbolic problems. Hyperbolic, or wave-like, equations propagate signals at a finite speed along characteristics; there is coherency in the wave form as it propagates, and therefore information is retained. Parabolic, or diffusion, equations propagate signals at infinite speed; because the signals diffuse or smear out, there is a gradual loss of information. A good way to understand how different equations propagate information is to determine how a signal is propagated in a special case. For example, suppose the initial signal is a Gaussian function or bellshaped curve exp( _x 2 ). Think of this signal as being a bit of information. The convection equation Ut + CUx = 0, which is a wave-like equation, propagates this signal via U(x, t)

= e-(X-ct)'.

That is, it moves it at speed c without distortion. The wave equation Utt = c2 u xx moves it via u(x, t) = O.S( e-(x-ct)'

+ e-(x+ct)\

So the signal breaks into two pieces, and they propagate in opposite directions at speed c. The diffusion equation Ut = kuxx propagates the signal via U(x, t) =

1

.J1

+ 4kt

e- X2/ (1+4kt) .

So the signal stays at the same place, but it spreads out and decreases in amplitude. Any information in the signal is eventually lost.

Exercises 1. Derive d'Alembert's formula (2.14) by determining the two arbitrary func-

tions F and G in the general solution (2.11) using the initial conditions (2.13). 2. Calculate the exact solution to the Cauchy problem when c = 2, the initial displacement isf(x) = 0, and the initial velocity isg(x) = 1/(1 + 0.25x2 ). Plot the solution surface and discuss the effect of giving a string at rest an initial impulse. Contrast the solution with the case when f #- 0 and g = O. 3. Solve the outgoing signal problem Utt -

c2 u"x

= 0,

X

> 0, -00 < t
(x) ,

x

E

R, t > 0,

x

E

R,

= kvxx , = 1/I(x) ,

x

E

R, t > 0,

Ut

and Vt V(x, 0)

x

E

R,

where 4> and 1/1 are continuous, bounded functions and close in the sense that 14>(x) - 1/I(x) 1::: 8 for all x, where 8 is a small number. We would like to show that the corresponding solutions u(x, t) and v(x, t) are close. Let us define w(x, t) = u(x, t) - vex, t) and note that w satisfies the Cauchy problem Wt

= kwxx ,

x

E

R, t > 0,

w(x, 0) = 4>(x) - 1/I(x) , x E R.

i:

The solution formula (2.8) gives w(x, t) =

(p(y) - 1/I(y))G(x - y, t)dy.

2.4.

Semi-Infinite Domains

63

f: : : f:

Therefore, for each t > 0,

I ¢(y) - 1/I(y) II G(x - y, t) I dy

lu(x, t) - vex, t)1 ::::

oG(x - y, t)dy = 0,

since J G(x - y, t)dy = l. Therefore, in the sense interpreted above, closeness of the initial data implies closeness of the solution.

Exercises 1. Show that the Cauchy problem for the backward diffusion equation,

+ U xx = 0, x U(x, 0) = f(x),

Ut

E R,

x

t > 0,

E R,

is unstable by considering the solutions u(x, t)

1

= 1 + -n e"

2

t

sin nx

for large n. 2. Let u

= u(x, y). Is the problem U xy = 0,

0 < x, y < I,

on the unit square, where the value ofu is prescribed on the boundary of the square, a well-posed problem? Discuss. 3. Consider two Cauchy problems for the wave equation with different initial data: i

U tt

ui(x, 0)

= C2 Uxxi , =

fi(x),

X E R,

u;(x,O)

0 < t < T,

= gi(X),

X E R,

for i = I, 2, where fl, f2, gl, and g2 are given functions (the superscripts are indices and not exponents). Iffor all x E R we have

1flex) - f2(X) I:::: Dl,

1gl(X) - g2(x) I:::: 8z,

show that 1 u l (x, t) - u 2(x, t) 1 :::: Dl + D2T for all x this mean with regard to stability?

2.4

E R,

0 < t < T. What does

Semi-Infinite Domains

In Sections 2.1 and 2.2 we solved the heat equation and the wave equation,

respectively, on the domain-oo < x < 00. Now we study these problems when the domain is semi-infinite, i.e., on the interval 0 < x < 00.

2. Partial Differential Equations on Unbounded Domains

64

This means that there is a boundary in the problem, at x = 0, and one fully expects that it is necessary to impose a boundary condition there. For example, to determine how the temperature distribution evolves in a semi-infinite bar, one should know the temperature in the bar initially, as well as the temperature at x = o. Therefore, we consider the initial boundary value problem for the heat equation Ut

= kuxx ,

U(O, t) = 0,

x > 0, t > 0,

(2.17)

t > 0,

(2.18)

u(x,O) = ifJ(x) , x > 0,

(2.19)

where we have specified the temperature to be zero at x = 0 for all time. Th solve this problem we use the method of reflection through the boundary. The idea is to extend the problem (2.17)-(2.19) to the entire real axis by extending the initial data ifJ to an odd function 1/1 defined by 1/I(x)

= ifJ(x) if

x > 0;

1/I(x)

=

-ifJ( -x) if x < 0;

1/1(0)

= o.

We then solve the extended problem by formula (2.8) and then restrict that solution to the positive real axis, which will then be the solution to (2.17)-(2.19). Figure 2.5 shows the initial data for the extended problem and a resulting odd solution profile vex, t). Physically, we are attaching a bar occupying the space -00 < x < 0 and giving it an initial temperature equal to the negative of that in the original bar. Therefore, let us consider the Cauchy problem Vt

= kvxx ,

x

E

R, t > 0,

(2.20) (2.21)

v(x,O) = 1/I(x) , x E R,

where 1/1 is the odd extension of the function ifJ as defined above. By formula (2.8) the solution to (2.20)-(2.21) is given by vex, t) =

i:

G(x - y, t)1/I(y)dy.

v

x V(x, t) (= u(x, t) for x> 0)

Figure 2.5. The initial data 1/I(x) and solution profile for the odd, extended problem.

2.4.

Semi-Infinite Domains

L: L:

65

Breaking up this integral into two parts, y < 0 and y > 0, we obtain

=

vex, t)

G(x - y, t)ljI(y)dy

==

-1

00

+

1 +1 00

00

G(x - y, t)¢( -y)dy G(x

= 1°O[G(X -

+ y, t)¢(y)dy + y, t) - G(x

G(x - y, t)ljI(y)dy

1

00

G(x - y, t)¢(y)dy

G(x - y, t)¢(y)dy

+ y, t)]¢(y)dy.

We restrict this solution to x > 0, and therefore the solution to (2.17)(2.19) is u(x, t)

= 1°O[G(X -

y, t) - G(x

+ y, t)]¢(y)dy.

The wave equation on a semi-infinite domain can be solved in the same manner. Consider the problem of the transverse vibrations of a string occupying x > 0 when the end at x = 0 is held fixed. The initial boundary value problem is (2.22) t > 0,

u(O, t) = 0, u(x,O)

= [(x),

(2.23)

= g(x),

Ut(x, 0)

(2.24)

x > O.

For x > ct (Le., ahead of the leading signal from the boundary) the interval of dependence lies in (0, 00), where the initial data are given; therefore, in this domain, the solution is given by d'Alembert's formula: u(x, t)

= -1 ([(x 2

ct)

+ [(x + ct)) + -1

2c

l

x ct

+ g(s)ds,

x-ct

x > ct.

(2.25)

The data given along the x = 0 boundary cannot affect the solution in the region x > ct, since signals travel outward from the boundary at speed c. See Figure 2.6. To solve the problem in the region 0 < x < ct we proceed as we did with the heat equation and extend the initial data [ and g to odd functions on the entire real axis. Therefore, let us consider the problem (2.26) v(O, t) v(x,O)

= =

0,

t > 0,

F(x),

(2.27)

Vt(x,O)

=

G(x),

x > 0,

=

-[(-x),

(2.28)

where F(x)

= [(x),

x > 0;

F(O)

=

0;

F(x)

x < 0,

66

2. Partial Differential Equations on Unbounded Domains

O 0 where the problem (2.22)-(2.24) is defined. The region x > ct is affected only by the initial data f and g and can be found by d'Alembert's formula.

and G(x)

= g(x),

x > 0;

G(O)

=

0;

G(x)

=

-g( -x), x < O.

By d'Alembert's formula the solution to (2.26)-(2.28) is vex, t) = -1 (F(x - ct) 2

+ F(x + ct)) + -1

2c

1

x+ct

x-ct

G(s)ds.

In the region x < ct this becomes (since x - ct < 0) vex, t)

= -21 (-f( -x + ct) + f(x + ct)) + -1 . 2c + -1

2c

l

0

1 0

x~ct

-g( -s)ds

x ct

+ g(s)ds.

If the variable of integration s in the first integral is replaced by -s, then the two integrals can be combined, and we may write 1 u(x, t) = vex, t) = - (f(x 2

o
ct or x < ct. Why does this reflection method work for the heat equation and wave equation? 1b reiterate, the solutions to the Cauchy problems for these two equations are odd functions if the initial data is odd. And the restriction of an odd solution to the positive real axis is the solution to the given initial boundary value problem. If this intuitive reasoning leaves the reader perplexed, then one can always verify analytically that the solutions we have obtained by this reflection method are, in fact, solutions to the given problems.

Exercises

67

If the boundary condition (2.18) along x = 0 in the heat flow problem (2.17)-(2.19) is replaced by a Neumann condition t > 0,

ux(O, t) = 0,

then the problem can be solved by extending the initial data to an even function. The same is true for the wave equation. We leave these calculations as exercises.

Exercises 1. Solve the problem Ut

= kuxx ,

U",(O, t) = 0, u(x,O)

x > 0, t > 0,

t > 0,

= cjJ(x),

x > 0,

with an insulated boundary condition by extending cjJ to all of the real axis as an even function. The solution is U(x, y)

=

i:

[G(x - y, t)

+ G(x + y, t)]cjJ(y)dy.

2. Find a formula for the solution to the problem Ut u(O, t)

= kuxx ,

= 0,

x > 0, t > 0,

t > 0,

U(x, 0) = 1, x >

o.

Sketch a graph of several solution profiles with k

= 0.5.

3. Find the solution to the problem

= 0, t u(x,O) = xe-", u(O, t)

> 0,

Ut(x,O)

= 0,

x >

o.

Pick c = 0.5 and sketch several time snapshots of the solution surface to observe the reflection ofthe wave from the boundary. 4. Solve the problem Ut = kuxx , x > 0, t > 0, U(O, t) = 1, u(x, 0)

= 0,

t > 0, x >

o.

Hint: 'Transform the problem to one of the form (2.17)-(2.19) and use Exercise 2.

2. Partial Differential Equations on Unbounded Domains

68

5. Use Lord Kelvin's argument, given in the mid 1860s, to estimate the age of the earth using a measurement of the geothermal gradient at the surface. Treat the earth as flat with x > 0 measuring the depth from the surface x = O. Assume that the diffusivity is k = 0.007 square centimeters per second, the initial temperature was 7000 degrees Fahrenheit (molten rock), and the temperature of the surrounding atmosphere has always been 0 degrees. Suppose the current geothermal gradient is 3.7 x 10-4 degrees per centimeter. About what percent of its original heat has the earth lost?

2.5

Sources and Duhamel's Principle

How do we proceed if the PDE contains a source term? For example, consider the heat-flow problem Ut

= kuxx

u(X, 0) = 0,

+ fCx, t),

x

E

R, t > 0,

X E R,

(2.30) (2.31)

where f is a given heat source. The key to the analysis of this problem can be discovered by examining an ordinary differential equation with a source term. For example, consider the initial value problem y' (t)

+ ay =

F(t),

t > 0;

yeO) =

o.

(2.32)

Multiplying by the integrating factor eat makes the left side a total derivative, and we obtain d _ (eaty)

dt Integrating from 0 to t then gives

= eatF(t).

eaty(t) = 1t ea'F(r)r,

which can be rewritten as yet) = 1t e-a(H)F(r)r.

(2.33)

Now let us consider another problem, where we put the source term in as the initial condition. Let w = wet; r) be the solution to the problem w'(t; r)

+ aw(t; r) = 0,

t > 0;

w(O; r)

= F(r) ,

where a new parameter r has been introduced. It is straightforward to see that the solution to this problem is wet; r) = F(r)e- at .

2.5.

Sources and Duhamel's Principle

69

So the solution to (2.32), the problem with a source, is the integral of the solution w(t, r) (with t replaced by t - r) ofthe associated homogeneous problem where the source is included as an initial condition. The fact that a particular solution of a linear equation can be deduced from the solution of the homogeneous equation is called Duhamel's principle. For ODEs we state the principle as follows:

Duhamel's principle for ODEs The solution of the problem

+ ay = F(t),

y'(t)

is given by

=

yet)

it

t > 0;

yeO)

=

0

wet - r, r)dr,

where w = w(t, r) solves the homogeneous problem w'(t; r)

+ aw(t; r) =

t > 0;

0,

w(O; r)

= F(r).

The same type of result is true for second-order ODEs as well; the reader may recall that the variation of parameters method uses the homogeneous solutions to construct a particular solution. See the Appendix on ODEs for the formula. Now let us extrapolate this idea and apply it to the heat flow problem (2.30)-(2.31). If Duhamel's principle is valid in this case, then the solution of (2.30)-(2.31) should be

it

u(x, t) =

w(x, t - r, r)dr,

where w(x, t; r) solves the homogeneous problem Wt

= kwxx ,

w(x, 0; r) = f(x, r),

x

E

x

R, t > 0, E

R.

(2.34) (2.35)

In fact, we can write down the explicit formula; by (2.8) the solution to (2.34)-(2.35) is w(x, t; r) =

i:

G(x - y, t)f(y, r)dy,

where G is the heat kernel. Therefore, the solution to (2.30)-(2.31) should be given by u(x, t) =

t 1 G(x -

Jo

00

y, t - r)f(y, r)dydr.

-00

Indeed, one can verify that this is the case. Observe again that r is being regarded as a parameter in the discussion above.

2.

70

Partial Differential Equations on Unbounded Domains

It is not too surprising that the solution turned out to be an integral. The PDE (2.30) has the form Hu = f, where H = k ;t22 is a differentialoperator. Ifwe formally write u = H-1f (as we might do in matrix theory if H were a matrix and u and f vectors), then we might expect H- 1 , the inverse of H, to be an integral operator, since integration and differentiation are inverse processes. We may now write down the formula for the solution to the problem

ft -

Ut = ku xx

+ f(x, 0,

u(x, 0) = ¢(x) , x

E

X E

(2.36)

R, t > 0,

R,

(2.37)

where the initial condition is no longer zero. By linearity, the solution to (2.36)-(2.37) is the sum of the solutions to the two problems Ut = ku xx , x

E

R, t > 0,

(2.38)

u(x,O) = ¢(x) , x

E

R,

(2.39)

and Ut

= ku,yx + f(x, t), x

u(x, 0) = 0,

i:

X E

R, t > 0,

(2.40) (2.41 )

R.

Thus, the solution to (2.36)-(2.37) is u(x, t) =

E

G(x - y, t)¢(y)dy + fut

i:

G(x - y, t - r)f(y, r)dydr. (2.42)

This is the variation of parameters formula for the problem (2.36)-(2.37), which is the analogue of (2.33) for the ODE (2.32). Duhamel's principle can also be applied to the wave equation. The solution to the problem (2.43)

u(x,O) = Ut(x, 0)

=

0,

X E

(2.44)

R,

is u(x, t)

where w

= fut w(x, t -

r, r)dr,

= w(x, t; r) is the solution to w(x, 0; r) = 0, Wt(x, 0; r) = f(x, r),

x

E

R.

We put the source term f into the initial condition for Wt rather than for W because f, like W t , is an acceleration; note that the displacement u is

Exercises

a time integral of w, so w must be a velocity, making W t an acceleration. From d'Alembert's formula, w(x, t; r)

= -1

2c

l

x ct

+ f( s, r)ds,

x-ct

and therefore by Duhamel's principle the solution to (2.43)-(2.44) is given by the formula u(x, t)

1t l

1

= -

2c

x +c(t-r)

f(s, r)dsdr.

(2.45)

x-c(t-r)

0

Source terms also arise in PDEs when problems are transformed in order to homogenize the boundary conditions. For example, consider the diffusion problem Ut = ku xx ,

u(x, 0)

=

x > 0, t > 0,

(2.47)

¢(x) , x > 0,

u(O, t) = g(t),

t >

(2.46)

o.

(2.48)

We solved this problem in the last section wheng(t) = O. So let us attempt to transform the problem into one where the boundary condition is zero. To this end, let vex, t) = u(x, t) - g(t), or u(x, t) = vex, t) + get); then substituting into (2.46)-(2.48) gives Vt

x > 0, t > 0,

= kv xx - g'(t),

vex, 0) = 1/r(x) v(O, t) = 0,

==

¢(x) - g(O),

x > 0,

(2.49) (2.50) (2.51 )

t > O.

Therefore, transformation of the dependent variable has changed the problem into one with a homogeneous boundary condition, but a price was paid-an inhomogeneity, or source term -g'(t), was introduced into the PDE. In general, we can always homogenize the boundary conditions in a linear problem, but the result is an inhomogeneous PDE; so inhomogeneous boundary conditions can be traded for inhomogeneous PDEs. We can solve (2.49)-(2.51) for vex, t) by formulating a Duhamel's principle. However, in Section 2.5 we will observe that Laplace transform methods can also be applied to find the solution.

Exercises 1. Write a formula for the solution to the problem Utt -

c2 u xx = sin x,

U(X, 0) = Ut(x, 0) = 0,

x

E R,

x E R.

t > 0,

2.

72

Partial Differential Equations on Unbounded Domains

Graph the solution surface when e

= 1.

2. Write a formula for the solution to the problem

= sin x, x E u(x, 0) = 0, X E R.

R, t > 0,

Ut - ku xx

3. Using Duhamel's principle, find a formula for the solution to the initial value problem for the convection equation u,

+ eux = f(x, t),

x E R, t > 0;

u(x,O)

= 0,

X E

R.

Hint: Look at the problem Wt(x, t; r)

+ ew(x, t;

= 0,

r)

x E R, t > 0;

w(x, 0; r)

= f(x, r),

x E R.

Solve the problem Ut

+ 2ux = xe- t ,

x E R, t > 0;

u(x, 0)

= 0,

X E

R.

4. Formulate Duhamel's principle and solve the initial boundary value problem u, u(x, 0) u(O, t)

The solution is u(x, t)

2.6

=[

[>0

= ku xx + f(x, t), = 0, x > 0, = 0, t > 0.

x > 0, t > 0,

(G(x - y, t - r) - G(x

+ y, t -

r))f(y, r)dydr.

Laplace Transforms

Laplace transforms are first encountered in elementary differential equations courses as a technique for solving linear, constant-coefficient ordinary differential equations; Laplace transforms convert an ODE into an algebra problem. The ideas easily extend to PDEs, where the the operation of Laplace transformation converts PDEs to ODEs. Let U = u(t) be a piecewise continuous function on t 2: 0 that does not grow too fast; for example, assume that u is of exponential order, which means that lu(t)1 ::::: C exp(at) for t sufficiently large, where a, C > O. Then the Laplace transform of u is defined by (LU)(S)

==

U(s) =

1

00

u(t)e-stdt.

(2.52)

The Laplace transform is an example of an integral transform; it takes a given function u(t) in the time domain and converts it to a new function U(s) , in the so-called transform domain. U and s are called the transform variables. The transform is linear in that L(CIU + C2V) = CILU + CzLV,

2.6.

Laplace'Ii'ansforms

where Cl and Cz are constants. If the transform U(s) is known, then u(t) is called the inverse transform of U(s) and we write C- 1 U = u. Pairs of Laplace transforms and their inverses are tabulated in books of tables, and many software packages contain commands that yield transform pairs. A short table is given at the end of this text. In Maple (version V4), the command laplace (sin{t) ,t,s); returns 1

1 + sz'

which is the Laplace transform of sin t. 1b obtain the inverse transform use invlaplace (1/{1+sA2),s,t); which returns sin t. The importance of the Laplace transform, like other transforms, is that it changes derivative operations to multiplication operations in the transform domain. In fact, we have (Cu')(s) = sUes) - u(O),

(2.53)

(Cu")(s) = szU(s) - su(O) - u'(O).

(2.54)

Formulae (2.53) and (2.54) are readily derived using integration by parts, and they are the basic operational formulae for solving differential equations. EXAMPLE

Solve the initial value problem u"

+u =

0,

t > 0;

u(O) = 0, u'(O) = l.

Thking Laplace transforms of both sides of the differential equation and using (2.54) gives szU(s) - 1

+ U(s) =

0,

or U(s)

1

= 1 + Sz •

So, the ordinary differential equation has been transformed into an algebraic equation, and we have solved the problem in the transform domain. 1b recover u(t) from its transform, we can look up the inverse transform

2. Partial Differential Equations on Unbounded Dpmains

74

in the table of Laplace transforms or use a computer algebra program to find u(t) = sin t,

o

which is the solution.

In summary, taking the Laplace transform of an ordinary differential equation for u(t) results in an algebraic equation for U( s) in the transform domain. We solve the algebraic equation for U(s) and then recoveru(t) by inversion. The same strategy applies to partial differential equations where the unknown is a function of two variables, for example u = u(x, t). Now we transform on t, as before, with the variable x being a parameter unaffected by the transform. In particular, we define the Laplace transform of u(x, t) by (.cu)(x, s)

==

= LX! u(x, t)e-stdt.

U(x, s)

(2.55)

Then time derivatives transform as in (2.53) and (2.54); for example, (.cUt)(x, s) = sU(x, s) - u(x, 0). On the other hand, spatial derivatives are left unaffected; for example,

=

1 -a

= -a

1

u(x, t)e-stdt = Ux(X, s). ax axo Therefore, upon taking Laplace transforms, a PDE in x and t is reduced to an ordinary differential equation in x; all the t-derivatives are turned into multiplication in the transform domain. (.cUx)(X, s)

00

o

u(x, t)e-stdt

00

EXAMPLE

Use Laplace transforms to solve the following initial boundary value problem associated with the diffusion equation. Let u = u(x, t) denote the concentration of a chemical contaminant dissolved in a fluid in the semiinfinite domain x > 0. Initially, assume that the domain is free from contamination. For times t > we impose a constant unit concentration of a contaminant on the boundary x = 0, and we ask how this contaminant diffuses into the region. Assuming a unit diffusion constant, the mathematical model is

°

= 0, u(x,O) = 0, u(O, t) = 1,

Ut - Uxx

x > 0, t > 0, x > 0, t > 0;

u(x, t) bounded.

Taking Laplace transforms of both sides of the PDE gives sU(x, s) - Uxx(x, s) = 0.

2.6.

Laplace Transforms

75

This is an ordinary differential equation with x as the independent variable, and the solution is

U(x, s) = a(s)e- Jsx

+ b(s)e Jsx .

Because we want bounded solutions, we must set b(s) =

o. Then

U(x, s) = a(s)e- Jsx . Now we take a Laplace transform of the boundary condition to get U(O, s) = 1/s, where we have used .eel) = l/s. Therefore, a(s) = l/s, and the solution in the transform domain is

U(x, s)

1

= - e- Jsx .

s Consulting the table or a computer algebra program, we find that the solution is

u(x, t) = erfc

(.Ju) ,

where erfc is the function defined by the formula erfc(y) = 1 -

2

fo

[Y

2

Jo e- r dr.

o

Observe that erfc(y) = 1 - erf(y).

In both preceding examples we were able recover a function from its Laplace transform by looking up the transform pair in a table or finding it in a software package. One might ask, in general, how to determine u(t) from knowledge of its transform U(s). It would take us far afield to give a thorough discussion. However, we can indicate the general formula used to compute the inverse transform, that is, to compute u(t) from its transform U( s). The inversion formula is

u(t)

1 = (.c- 1 U)(t) = -. 27ft

l

a ioo

+

a-icc

U(s)estds.

The integral in this formula is a complex contour integral taken over the infinite straight line (called a Bromwich path) in the complex plane from a - ioo to a + ioo. The number a is any real number for which the resulting Bromwich path lies to the right of any singularities (poles, essential singular points, or branch points and cuts) of the function U(s). A thorough discussion of the inversion formula can be found in the references. Another result that is useful in calculations is the convolution theorem.

Theorem Let u and v be piecewise continuous on t :::: 0 and of exponential orda Then .c(u

* v)(s) =

U(s)V(s),

2.

76

Partial Differential Equations on Unbounded Domains

where (u

* v)(t) ==

1t

is the convolution ofu and v, and U

u(t - t)v(r)dr

= Cu, V = Cv.

We note that the Laplace transform is additive, but it is not multiplicative. That is, the Laplace transform of a product is not the product of the Laplace transforms. The convolution theorem tells what to take the transform ofin order to get a product of Laplace transforms, namely, the convolution. EXAMPLE

In the contaminant transport model in the last example let us change the boundary condition to a function of time and consider Uxx

= 0,

x > 0, t > 0,

u(x, 0)

= 0,

x > 0,

Ut -

u(O, t) = fl.t),

t > 0;

u(x, t) bounded.

Thking Laplace transforms of both sides of the PDE gives, as in the example, sU(x, s) - Uxx(x, s)

= 0,

which has solution U(x, s) = a(s)e-~sx.

Now we take a Laplace transform of the boundary condition to get U(O, s) = F(s). Therefore, a(s) = F(s) and the solution in the transform domain is U(x, s)

= F(s)e-~sx.

Consulting a table or Maple, we find that C- 1 (e-~sx) = _X_e-X2/4t .

.J41l't3

Therefore, we can use the convolution theorem to write the solution as

u(x, t) =

1J t

o

X

41l'(t - r)3

2

e- x 1(4(t-r))f(r)dr.

o

Contained in the Maple software program there is a package (the integral transform package) that can be adapted to solve PDEs using Laplace transforms. We discuss this procedure in Section 2.8.

Exercises

77

A very complete reference for the Laplace transform that includes an extensive table, theory, and applications is Churchill (1958). A short table appears in the back of this text.

Exercises 1. Solve the following problem using Laplace transforms.

u(O, t) = 0,

t > 0,

U(x,O) = Ut(x, 0) = 0,

x > O.

The solution shows what happens to a falling cable lying on a table that is suddenly removed. Draw some time snapshots of the solution. 2. In the quarter plane x, y > 0, where the temperature is initially zero, heat flows only in the y-direction; along the edge y = 0 heat is convected along the x-axis, and the temperature is constantly 1 at the point x = y = O. The boundary value problem for the temperature u(x, y, t) is Ut = Uyy ,

Ut(x, 0, t)

x,

t,y > 0,

U(x,y,O) = 0,

x,y > 0,

u(O, 0, t) = I,

t > 0,

+ ux(x, 0, t)

= 0,

x, t > O.

Find a bounded solution to this problem using Laplace transforms. 3. Show that £ (J~ f(r)dr) = F(s)/s. 4. Show that £ (H(t - a)f(t - a)) = e-asF(s), where H is the unit step function defined by H(x) = 0 for x < 0, and H(x) = 1 for x ::: O. 5. A very deep container ofliquid is insulated on its sides. Initially, its temperature is a constant Uo degrees, and for t > 0 heat radiates from its exposed top surface according to Newton's law of cooling (see Section 1.3). The air temperature is zero degrees. Formulate an initial boundary value problem for the temperature of the liquid, and find a formula for the temperature at various depths at various times. Use p = c = K = fJ = 1 and use a computer algebra package to graph some temperature profiles. 6. Derive the solution u(x, t)

=

H(t - xl c)g(t - xl c) to the problem

U(x,O) = Ut(x, 0) = 0, u(O, t)

= get),

t > O.

x > 0,

78

2.

Partial Differential Equations on Unbounded Domains

2.7 Fourier Transforms The Fourier transform is another integral operator with properties similar to the Laplace transform in that derivatives are turned into multiplication operations in the transform domain. Thus the Fourier transform, like the Laplace transform, is useful as a computational tool in solving differential equations. In PDEs the Laplace transform is usually applied to the time variable; the Fourier transform is usually applied to the spatial variable when it varies over (-00, 00). First let us begin with functions of one variable. The Fourier transform of a function u = u(x), x E R, is defined by the equation

(Fu)(~) == um =

i:

u(x)ei~xdx.

(2.56)

If u is absolutely integrable, Le., J~oo luldx < 00, then u can be shown to exist. However, in the theory of Fourier transforms, it is common to work with a a smaller set of functions; the integrals involved are improper integrals, and so the functions must decay rapidly enough in order for the integrals to exist. We define S as the set of rapidly decreasing functions on R that have continuous derivatives of all orders; the rapidly decreasing functions are functions that, along with all their derivatives, decay to zero as x -+ ±oo faster than any power function (functions like lIx 2 and llx 6 ). The function exp( _x2 ) defining the bell-shaped curve is a such a rapidly decreasing function. More technically, if the set of functions that have continuous derivatives of all orders on R is denoted by Coo, then S = {u

E Coo:

k

dku dx k

I-

I< -

M-

1

IxlN

as Ixl

-+ 00

'

= 0,1,2, ... ; for all integersN}.

The set S is called the Schwartz class of functions, and one can show that if U E S, then E S, and conversely. So S is a closed set under Fourier transformation, which makes it a good set to work with. There is one important remark about notation. There is no standard convention on how to define the Fourier transform (2.56); some put a factor ofl/(2:rr) or 11.j2ii in front of the integral, and some have a negative power in the exponential, or even a factor of 2:rr. One should be aware of these variations when consulting other sources. A basic property of the Fourier transform is that the kth derivative UCk) (k = 1,2, ... ) transforms to an algebraic expression. That is,

u

(2.57) confirming our comment that derivatives are transformed to multiplication (by a factor of (-i~l). This formula is easily proved by integration by parts (as for the Laplace transform); all the boundary terms generated

2.7.

Fourier 'fiansforms

79

in the integration by parts are zero, since u and all its spatial derivatives vanish at ±oo. For functions of two variables, say u = u(x, t), the variable t acts as a parameter, and we define

(Fu)(~, t) == u(t t) =

f : u(x, t)eii;xdx.

Then, under Fourier transformation, x-derivatives turn into multiplication, and t derivatives remain unaffected; for example, (Fux)(~, t)

= (-i~)u(~, t),

(Fuxx)(t t) = (-i~iu(~, t), (Fut)(~, t)

=

Ut(t t).

Solving a differential equation for u involves first transforming the problem into the transform domain and then solving for U. Then one is faced with the inversion problem, or the problem of determining the u for which Fu = U. Another nice property ofthe Fourier transform is the simple form of the inversion formula, or inverse transform. It is (F-1u)(x)

==

u(x)

= -1

2n

1

00

. u(~)e-ll;x~.

(2.58)

-00

This result is called the Fourier integral theorem; it dictates how to get back from the transform domain. Some Fourier transforms can be calculated directly; many require complex contour integration. In the following example we calculate the transform of the Gaussian function u(x) = e- ax2 , a > 0, using a differential equation technique. EXAMPLE

We want to calculate U, where

[!(~) = Differentiating with respect to u'm

=

~

f : e- ax' eil;xdx.

and then integrating by parts gives

i f : xe- ax2 eil;xdx -i 2a

1

00

-00

~ e-ax2 eil;xdx dx

-~ U(~). 2a

Therefore, we have a differential equation [t' = ~! variables and integrating gives the general solution um

= Ce-1;2!(4a).

u for U.

Separating

2.

80

Partial Differential Equations on Unbounded Domains

The constant C can be determined by noticing that u(O) =

1

00

V~. -;.

e- ax2 ax =

-00

Consequently, we have (2.59)

So, the Fourier transform of a Gaussian function is itself a Gaussian; likewise, the inverse transform of a Gaussian is a Gaussian. 0 A convolution relation also holds for Fourier transforms. If u, v then we define their convolution, which is in S, by

(u

* v)(x) =

i:

E

S,

u(x - y)v(y)dy.

Then we have the following convolution theorem:

Theorem Ifu,v E S, then :F(u

* v)(~) = umvm.

By the Fourier integral theorem it follows immediately that (u

* v)(x) =

:F-l(U(~)V(~)).

This formula states that the inverse transform of a product is a convolution; this is a useful relationship in solving differential equations. EXAMPLE

Let f

E

S and determine u for which u" - u

= f(x),

X E

R.

Thking the transform of both sides yields (_i~)2U -

u = 1,

or il(~)

=-

1

A

- - 2 f(~)'

l+~

In the transform domain the solution is a product of transforms, and so we apply the convolution theorem. From Exercises 2 and 3 we have :Fe ~ e-iXi) 2

=

1

1

+ ~2

.

2.7.

Fourier'Ii'ansfonlls

81

Therefore, u(x)

= - -1 e- 1xl * rex) = - -1 2

2

1

00

e-1x-y1rey)dy.

o

-00

The strategy in applying transform methods to solve differential equations is to proceed formally, making any assumptions that are required to obtain an answer; for example, assume that all the data is in S. When a solution is obtained one can then attempt to verify that it is indeed a solution to the problem. Often one can prove that the solution obtained holds under less severe conditions than are required in application of the transform method. Now we apply the Fourier transform method to the Cauchy problem for the heat equation. We will derive the same solution formula (2.8) that we obtained in Section 2.1 by a different method. EXAMPLE

Use Fourier transforms to solve the pure initial value problem for the heat equation: Ut - kuxx

= 0,

X E

Again we assume that [

E

R, t > 0;

u(x, 0)

= [(x),

X E

R.

(2.60)

S. Thking Fourier transforms of the PDE gives

Ut =

-~2ku,

which is an ordinary differential equation in t for parameter. Its solution is

u(~,

t), with

~

as a

= Ce-1;2kt. But the initial condition gives u(~, 0) = 1m, and so C = lc~). Therefore, U(~, t)

u(~, t) = e-1;2ktlc~).

Replacing a by 1/(4kt) in formula (2.59) gives

:F

(_1_ e-X2/(4kt)) .J4:rrkt

= e-1; 2kt .

Thus, by the convolution theorem we have u(x, t) =

1

00

-00

1 2 - - e-(x-y) 1(4kt)rey)dy.

.J4:rrkt

(2.61)

This solution was derived under the assumption that [ is in the Schwartz class. But now that we have it, we can attempt to show that it is a solution under milder restrictions on f. For example, one can prove that (2.61) is 0 a solution to (2.60) if [ is a continuous, bounded function on R.

2. Partial Differential Equations on Unbounded Domains

82

EXAMPLE

Now we work a problem involving Laplace's equation in the upper half plane. Consider U xx

+ Uyy = 0,

X E

R, Y > 0;

u(x,O)

= [(x),

X E

R.

We also append the condition that the solution U stay bounded as y ~ 00. This example is similar to the last example. Thking the transform (on x with y as a parameter) of the PDE, we obtain A U YY

c -5

U= 0 ,

2A

which has general solution

+ b(~)e~Y.

u(~, y) = a(~)e-~Y

The boundedness condition on U forces ~ < O. So we take u(~, y)

b(~)

= 0 if ~ > 0 and a(~) = 0 if

= c(~)e-I~IY.

Upon applying Fourier transforms to the boundary condition, we get c(~) = lc~). Therefore, the solution in the transform domain is

u(~, y) = e-'~'Ylm. Therefore, using the convolution theorem and Exercise 2, we obtain the solution U

x y =}£

(,)

7T:

* [=

1

x 2 + y2

1

}£

00

7T:

-00

JC:r)d"C

(x - "C)2

+ y2 .

o

An excellent introduction to Fourier transforms can be found in Strichartz (1994).

Exercises 1. Find the convolution of the functions [(x) = x and g(x) = e-"z. The answer

is xI./2.

2. Show that the inverse Fourier transform of e-al~l, a > 0, is

a 1 -; x2 + a2

.

3. Verify the following properties of the Fourier transform: (a) (Fu)(~) = 21l'(F- 1 u)( -~). (b) F(eiaxu)(~) = u(~ + a). (c) F(u(x + a)) = e-ia~um.

Exercises

83

Formula (a) states that if a transform is known, so is its inverse, and conversely. 4. Find the Fourier transform of the function u defined by u(x) and u(x) = 0 if x .::: o. 5. Compute F(xe- ax

2

= e-ax if x

> 0,

).

6. Solve the following initial value problem for the inhomogeneous heat

equation: Ut

= Uxx + F(x, t),

x E R, t > 0 u(x,O)

= 0,

X E R.

7. Find a formula for the solution to the following initial value problem for the

free Schrodinger equation: Ut = iuxx ,

u(x, 0) = e- x

x E R, t > 0;

2

,

x E R.

8. Find a bounded solution to the Neumann problem

+ U yy = 0, x uy(x, O) = g(x),

U xx

Hint: Let v is

=

R, Y > 0,

E

x E R.

uy and reduce the problem to a Dirichlet problem. The solution

u(x, Y)

= -1

1

00

2:rr

g(x -~) In(y2

-00

+ e)d~ + c.

9. Solve the boundary value problem

u xx

+ Uyy =

0,

U(x, 0) = I,

X E

Ixl .:::

R, Y > 0,

Ixl

u(x, 0) = 0,

I;

> I,

where b > 0 is a constant. 10. Use integration by parts to verify (assume (Fux)(t t)

=

U E

S)

(-i~)u(~, t),

= (-i~iu(~, t). u(x) = 1 if Ixl .::: a and u(x)

(Fuxx)(~, t) 11. Let u( x) be a square wave, i.e,

Show that

(Fu)(~)

=

o if Ixl

> a.

2 sin a~ -~-.

12. Solve the Cauchy problem for the following convection-diffusion equation

using Fourier transforms: Ut

= Duxx

-

CU x ,

X E

R, t > 0;

U(x, 0)

13. This exercise explores the role of a term PDE by examining the equation Ut

+ U xxx

=

U xxx

=

0): assume (k > 0) : invlaplace(U(x,s),s,t)i gives u(x, t) = erfc

(~) .

o

Of course, using a computer algebra software package like Maple is not the only way to use computers to solve PDEs. The reader should recall from an elementary course in ODEs that there are methods to numerically solve ordinary differential equations (Euler's. method, the Runge-Kutta method, and so on). Similar methods, based on differencequotient approximations of the partial derivatives, can be developed for partial differential equations. We shall examine some of these methods in later sections of this text.

Exercises 1. Use pdesol ve in Maple, or the appropriate command in other packages, to find the general solution of the following PDEs: (a) Utt - 5u = o. (b) Utt - Uxt - 6uxx = o. (c) Ut + x + U = o. Cd) Uxt + Ux = xte- t •

-d:xzu

2. Use PDEplot to graph the solution surface of the PDE Ut

+ cos(2t)ux =

- sin x

that satisfies the initial condition u(x, 0) = 1 + x 2 .

2.

90

Partial Differential Equations on Unbounded Domains

3. Consider an age-structured population where U = u(x, t) is the density of people at age x at time t. In a model where no person survives past age x = d, the PDE for u is given by c Ut + Ux = - d-- u, < x < d, t > 0.

°

-x

Justify this model and give an interpretation of the positive constant c. Use pdesolve to find the general solution to the PDE. Next, use analytical techniques to find the solution that satisfies the initial and boundary conditions U(x,O)

= p(x),

0::::: x ::::: d;

u(O, t)

= bet),

t > 0,

where p(x) is the density initially and bet) is the birth rate. Hint: Determine the arbitrary function differently for the separate regions x > t and x < t. Note that Maple was not careful in choosing the correct branch of the logarithm for the domain of interest. 4. The Fourier transform of a function u(x, t) on x is defined in Maple by the command four i er (u (x ,t) X ~) where ~ is the transform variable. Use Fourier transforms, in the same way as we used Laplace transforms in the example in the text to solve the initial value problem I

Ut - ku xx = 0,

I

x E R, t > 0;

,

u(x,O) = e

5. Use pdesol ve to find the general solution to Ut the correct answer?

_x 2

, x

+ Ux =

E R.

u 2. Did Maple give

Orthogonal Expansions CHAPTER

3.1

The Fourier Method

The basic technique for solving PDEs on a bounded spatial domain is the Fourier method, named after Joseph Fourier (1768-1830). In this section we take a nineteenth-century perspective and make some comments about the origin of the method. Our discussion will motivate one of the fundamental topics in analysis and in PDEs, namely orthogonal expansions. To fix the notion let us consider heat flow in a finite bar of length n and unit diffusivity (k = 1), where the ends are held at zero degrees. We choose n as the length to make the constants come out simpler at the end. From preceding sections we know that the temperature u = u(x, t) in the bar must satisfy the model equations Ut

u(O, t)

= U xx , 0 < x = u(n, t) = 0,

< n;

> 0,

t > 0,

(3.1) (3.2)

consisting of the heat equation and the boundary conditions. It is easy to check that for any positive integer n a solution to (3.1) and (3.2) is given by 2

un(x, t) = e- n t sin nx.

These solutions clearly satisfy the initial conditions un(x,O)

= sin nx, 91

J. D. Logan, Applied Partial Differential Equations © Springer-Verlag New York, Inc. 1998

3.

92

Orthogonal Expansions

which give temperature profiles at t = O. But what if the initial temperature profile is not sin nx, but rather an arbitrarily given function rex)? That is, what is the solution to (3.1) and (3.2) augmented by the initial condition u(x,O) = rex),

0 < x < n?

(3.3)

Fourier argued that the solution u(x, t) to (3.1)-(3.3) could be taken as an infinite linear combination of the un(x, t), that is,

L a e00

u(x, t) =

n

n2t

sin nx

(3.4)

n=l

for appropriately chosen constants an. 1b satisfY the initial condition, (3.3) would require that

L an sin nx, 00

u(x, 0) = [(x) =

(3.5)

n=l

that is, the initial temperature function [(x) would have to have a representation in terms of the periodic functions sin x, sin 2x, sin 3x, .... 1b determine such an one could proceed formally and multiply both sides of (3.5) by sin mx, for some fixed but arbitrary index m, and then integrate from x = 0 to x = n to obtain

i rr o

[(x) sin mx

ax = !orr L an sin nx sin mx ax. 00

0

n=l

Now if we interchange the order of integration (we alert the reader that this is not always a valid operation), we obtain

!o rr rex) sin mx ax = o

L an i rr sin nx sin mx ax. 00

n=l

0

Now we observe that the infinite series collapses to a single term. This is because of the integration formula

!orr sin nx sin mx ax =

0,

n

=1=

m.

(3.6)

So the only term that survives in the infinite series is the one where the summation index n hits m, the fixed index. Then we get

r [ex) sin mx ax = am 10r

10

sin 2 mx

ax = :: am. 2

But m is an arbitrary index, and so we have shown that the coefficients in (3.4) are given by an

2irr

=-

n

0

rex) sin nx

ax,

n = 1,2, ....

(3.7)

Exercises

93

So, in a formal sense (a formal calculation in mathematics is one done without complete rigor) the solution to the heat flow problem (3.1)-(3.3) is given by (3.4), where the coefficients an are given by (3.7). Indeed, one can verify that this is a solution for reasonable initial temperature distributions [(x). What we have just described is Fourier's method; namely, we solve a general initial boundary value problem by superimposing the solutions to the pure boundary value problem and then choosing the constants such that the initial condition is satisfied as well. One key to the calculation was that the functions sin nx satisfied a relation like (3.6), which is called an orthogonality condition. It enabled us to easily calculate the coefficients in the series (3.5). The infinite series (3.5) for [(x) in terms of the orthogonal functions sin nx is called an orthogonal expansion, or Fourier series. In the sequel we shall examine these ideas in detail. On another historical note, this calculation was troublesome in Fourier'S time, in the early 1800s. Concepts like convergence of series were not well understood. L. Euler, D. Bernoulli, and d'Alembert, all of whom had addressed similar problems in the mid-1700s regarding the wave equation and vibrating strings, had wondered about the possibility of expanding a nonperiodic function in terms of periodic functions like sines and cosines. It was L. Dirichlet who in 1829 established conditions under which such series representations are valid.

Exercises 1. Consider the initial boundary value problem for the wave equation (3.8) U(O, t)

= u(Jr, t) = 0,

U(x,O) = f(x),

(3.9)

t > 0,

Ut(x,O) = 0,

°< x
~ = IIfII2.

n=l

This equation is called Parseval's equality. Pointwise convergence results, or even stronger uniform convergence results, are more difficult to obtain. We can extend the modal interpretation of the generalized Fourier series still further. We can think of f as a signal and the orthonormal set fn as fundamental modes. The Fourier coefficient Cn determines the contribution of the nth mode, and the generalized Fourier series is the decomposition of the signal into fundamental modes. The sequence of squared coefficients, cr, c~, c~, ... , is called the energy spectrum, and c~ is called the energy of the nth mode; by Parseval's equality, the total energy in the signal is IIf 112 .

Exercises

101

Exercises 1. Using a table of integrals or a computer algebra package, verify that the set of functions I, cos x, cos 2x, ... form an orthogonal set on the interval [0, Jr). Next verify that the set of functions cos(nJrx/l), n = 0, 1,2 ... , form an orthogonal set on the interval [0, I). If

L Cn cos(nJrx/l) 00

fCx) =

n=O

in the mean-square on [0, I), what are the formulae for the cn ? 2. (Gram-Schmidt orthogonalization) We know from elementary linear algebra that any set of linearly independent vectors may be turned into a set of orthogonal vectors by the Gram-Schmidt orthogonalization process. The same process works for functions in LZ[a, b). Let h, fz, h ... be an independent set of functions in L z. Define the sequence gn by gl = h, gz = fz - ~;;W gl, g3 = F. - "jjg2jj2 (iJ.gzl gz - flg;"jjT (h,gJ) gl, . . .. Sh ow th at )3 gn·IS an 0 rth ogonaI sequence. 3. The functions I, x, xZ, x 3 are independent functions on the interval [-1, I). Use Exercise 2 to generate a sequence of four orthogonal polynomials on [-I, I). Denote the polynomials by Po (x) , ... ,P3(X) (they are called Legendre polynomials). Find an approximation of eX on [-1, 1) of the form C ~ coPo(x)

+ CIPI(X) + czPz(x) + C3P3(X)

that is best in the mean-square sense, and graph eX and the approximation on a set of coordinate axes. What is the pointwise error? What is the maximum pointwise error over [-1, I). What is the mean-square error? C

4. For which powers r is the function [(x)

= xr in LZ[O, 1J? In LZ[O, oo)?

5. Let [be defined and integrable on [0, I). The orthogonal expansion nJrx L bn sin -1-' 00

bn =

211

1

n=1

nJrx [(x) sin -1- dx,

0

is called the Fourier sine series for [ on [0, I). Find the Fourier sine series for [(x) = cos x on [0, Jr12). What is the Fourier sine series of fCx) = sin x on [0, Jr)? 6. (Hermite functions) Consider the differential equation -y"

+ xZy = Ey,

x E R,

E

= constant,

(3.16)

and the functions n xZ

Hn(x)=(-I)e

d n _Xz dxne ,n=O,I,2, ....

(a) Find Ho(x), ... ,H4(X). Why is Hn(x) always a polynomial? (These are called Hermite polynomials.) (b) Verify that vn(x) = Hn(x)e- xz / z is a solution of(3.16) when E = 2n + 1. Hint: You wiU need to show H:, = 2nHn-l.

3.

102

Orthogonal Expansions

(c) Show that f~oo v"vmdx = 0, m =1= n, and thus the v" are orthogonal on the interval (-00,00). Hint: Use the fact that both v" and Vm satisfy (3.16) for appropriate E. (d) If a function f(x) can be represented by f(x) = .L~=o cnv,,(X), how would you expect to find the c,,? Assume uniform convergence. Take f(x) = 11 ~ and use your calculator to find Co, ... , C4. (e) The differential equation (3.16) is a time-independent Schr6dinger equation (see Section 1.5). Thus E = Zn + 1 can be interpreted as the energy associated with the probability density v,,(xi, appropriately normalized so that f v~dx = 1. Sketch the normalized probability densities for n = 0, 1, Z. 7. (The Haar wavelets) In this exercise we explore the properties of a complete orthonormal system on (-00,00) called the Haarwavelets. Let¢ be a function such that ¢(x) = 1 for x E [0,1), and ¢(x) = otherwise. Let o/(x) = ¢(Zx) ¢(Zx - 1). Then the Haar wavelets are the functions

°

o/m,,(x)

=

Zm!2o/(Zm x - n),

for m, n = 0, ±l, ±Z, ., .. Sketch a graph of o/(x) , and then sketch a graph of o/nJ"(x) for m, n = 0, ±1, ±Z. What is the graph of o/mn(x)? If f is square integrable and [(x)

=

L L 00

00

Cmno/mn(x),

m=-oo n=-oo

find a simple formula for the coefficients Cm,,' 8. For f, g

E

L2[a, bj, prove the Cauchy-Schwarz inequality

I (f,g) I::: 11[llllgll· Hint: Proceed as follows. Define q(t) = Cf + tg, [ + tg) for any real number t. Using the rules of scalar product, expand this expression and obtain a quadratic polynomial in t; because q(t) 2: (why?), the quadratic polynomial can have at most one real root. Look at the discriminant of the polynomial.

°

3.3

Classical Fourier Series

In the last section we introduced the idea of representing a given function [(x) in terms of an infinite series of orthogonal functions [" (x). Now we want to focus those concepts on a special set of orthogonal functions where the [" are given by sines and cosines. The resulting series is called a (classical) Fourier series. Such series were introduced by Joseph Fourier in the early 1800s. Classical Fourier series are a special case of the orthogonal series discussed in the last section. We will usually drop the adjective "classical." Equation (3.11) is an example of a Fourier series on the interval [-If, lfJ. We will work on an arbitrary symmetric interval [-I, IJ about the origin.

3.3.

Classical Fourier Series

103

If [ is an integrable function on [-1, I], then its Fourier series is ao

nnx

00

nnx

2" + L(an cos -1- + bn sin -1-)'

(3.17)

n=1

where the Fourier coefficients an and bn are given by the formulae an bn

/1 [(x) cos -1nnx dx, /1 nnx = T _/(X) sin -1- dx, 1 = T

-I

1

n

= 0, 1,2, ... ,

(3.18)

n

= 1,2, ....

(3.19)

n=I,2, ... ,

(3.20)

Here, the set of functions 1

2'

nnx. nnx cos-I- , sln-I- ,

are orthogonal on the interval [-I, I]; they are playing the role of the [,,(x) in formula (3.12), and the an and bn are playing the role of the coefficients cn . It is shown in more advanced texts that the set of functions (3.20) is complete, and therefore the Fourier series (3.17) converges in the mean-square sense to [ when [ E L2[ -I, I]. EXAMPLE

Let [(x) = x for -I computed to be

:s

x

T1 /1-I X cos

an

=

bn

= -1 I

:s

I. Then the Fourier coefficients are easily

nnx -1- dx = 0,

n

= 0, 1,2, ... ,

/1 x sin -nnxd x = -(-1),,+1, 21 -I

I

nn

n = 1,2, ....

So the Fourier series for [ is 21 nx 1 2nx . 3nx - (sin - - - sin + -1sm - ... ). n 1 2 1 3 1 We make two observations. First, at x = ±1 the series clearly converges to zero and thus does not converge to [(x) at these two points. So the series does not converge pointwise to [(x) on the interval [-I, 1]. It does, however, converge pointwise to [(x) = x on (-1, l). Second, the derived series, obtained by differentiating term by term, does not converge at all, much less to the derivative f' (x) = 1. So the series cannot be differentiated term by term. We do know that the series converges to [ in the mean-square &~.

For the Fourier series (3.17) we define the nth harmonic to be nnx nnx an cos -1- + bn sin -1-'

D

3. Orthogonal Expansions

104

Notice that each harmonic in the series has a higher frequency (and thus more oscillations) than the preceding harmonic. The frequency spectrum of the series is the sequence of numbers Yn defined by

_ laol

Yo -

../2'

_ /

2

Yn - V an

+

b2 n

(n :::: 1).

The frequency spectrum is a measure of the contribution of the various harmonics in the decomposition of f. The numbers n:::: 0, form the energy spectrum. The reader is invited to show that Parseval's formula takes the form

Y;,

The frequency spectrum for [(x) = x in the last example is Yo = 0, Yn = 2l1(mr) , n:::: 1. It is often graphed, as in Figure 3.1, to show visually the contribution of each harmonic. A Fourier series simplifies considerably if [ is an even or an odd function. First observe that sin n7x is an odd function and cos n7x is an even function; moreover, an even function times an odd function is odd. Therefore, if [ is an even function, then the product [(x) sin n7x, which is the integrand in the expression for bn in (3.19), is an odd function; so all of the coefficients bn are zero because an odd function integrated over a symmetric interval about the origin is zero. Hence, if [ is an even function, then its Fourier series reduces to a series of the form aD

~

mrx

2" + .i...J an COS -z- . n=l

I.

frequency

3

n

Figure 3.1. Frequency spectrum for the function [(x) = x on [-Jr, Jr].

3.3. Classical Fourier Series

105

Similarly, iff is an odd function, then the coefficients an in (3.18) vanish, and the Fourier series reduces to ~

mrx

~bnsin -Z-. n=l

This fact is illustrated in the last example, where we found the Fourier series for f(x) = x, an odd function. We commented earlier that Fourier series are a valuable tool in signal processing and data storage. Let us explore this idea in more detail. For this discussion it is helpful to think of the variable x as representing time. We observe that the Fourier series (3.17), although computed only on the interval [-Z, Z], is a 2Z-periodic function* because the sines and cosines are 2Z-periodic. Thus, it repeats its values in every interval of length 21. Therefore, if f is defined on all of R and is 2Z-periodic, we cap. represent it by its Fourier series on all of R. Hence, iff is some given signal or data set 21-periodic in time (say a signal from an electrical circuit, an electrocardiogram, or an extraterrestial), then we could digitize the signal and save it by storing the Fourier coefficients or some finite subset of them; knowledge of those coefficients would allow us reproduce the signal via (3.17)-(3.19). EXAMPLE

Consider the 2-periodic signal shown in Figure 3.2, which is called a triangular wave. Analytically it is given by f(x) = x + 1 if -1 < x :s: 0; f(x) = 1 - x if 0 < x :s: 1; and otherwise 2-periodic. We compute its Fourier series. Here, f is an even function, and so bn = 0 for all n. The coefficients an are given by an =

11

f(x) cos mrxdx

-1

=

1 0

-1

(x

+ 1) cos mrxdx +

t 10

(1 - x) cos mrxdx.

When n = 0, we easily get ao = 1. For n ~ 1 we can calculate an by hand (using integration by parts), use a calculator, or use a computer algebra package. The next five coefficients are a1

= 0.405, az = 0,

a3

= 0.045,

a4

= 0,

as

= 0.016.

Therefore, a four-term Fourier approximation to the triangular wave is f(x)

~

0.5

+ 0.405 cos:n:x + 0.045 cos 3:n:x + 0.016 cos S:n:x.

Figure 3.2 shows a plot off and its approximation. From the last section we know that this approximation is the best, in the mean-square sense. D 'Recall that a function f is P-periodic, or periodic of period P, if it repeats itself in every interval oflength P; that is, fex + P) = fex}for all x.

3.

106

Orthogonal Expansions

Figure 3.2.

A

triangular wave.

-3

Finally, we make some brief comments about convergence of Fourier series. We have already noted that if f is square integrable, then meansquare convergence is automatic. Tb obtain pointwise convergence results some additional assumptions must be made regarding the smoothness of f. Our assumption is that the graph off is made up of finitely many, not necessarily continuous, smooth pieces; this will cover most of the interesting functions encountered in science and engineering. A function f is piecewise continuous on [a, b] if it is continuous except possibly at a finite number of of points in [a, b] where it has simple jump discontinuities; f has a simple jump discontinuity at x = c ifboth one-sided limits of fex) exist and are finite at c. The function mayor may not be defined at a jump discontinuity. We say f is piecewise smooth on [a, b] if both f and f' are piecewise continuous on [a, b], and we say f is piecewise smooth on (-00,00) if it is piecewise smooth on each bounded subinterval [a, b] of the real line. Then the basic pointwise convergence theorem is the following: Theorem Iff is piecewise smooth on [ - 1,1] and otherwise 21-periodic, then its Fourier series (3.17) converges pointwise for all x E R to the value f(x) iff is continuous at x, and to the average value of its left and right limits at x, namely ~ ([(x -) + fex + )), iff is discontinuous at x. Tb get stronger convergence results, like uniform convergence, additional smoothness conditions on f are required. Note that continuity off

0.& 0.60.4-

o.~

-8

-4

-2

00

Figure 3.3.

Graph of the 271:periodic square wave.

Exercises

107

is not enough to guarantee pointwise convergence; incredible as it may seem, there are continuous functions whose Fourier series diverge at every point!

Exercises 1. Find the Fourier series for the 2Jr-periodic square wave shown in Figure 3.3.

Sketch a two-term, a four-term, and a six-term approximation. 2. Find the Fourier series for f(x) = x 2 on [-Jr, Jrj. Sketch a graph ofthe function defined on all ofR to which the Fourier series converges. Is the convergence pointwise? Use the Fourier series to show Jr2 1 1 -=1--+ 12 4 9

1

+ ....

16

Graph the frequency spectrum. 3. Consider the function (signal) defined by f(x) = x + 1 for -2Jr < x :::: 0 and f(x) = x for 0 < x :::: 2Jr, and otherwise 4Jr periodic. Sketch a graph of the signal and find its Fourier series. Find and graph the frequency spectrum. To what value does the Fourier series converge at x = O? At x = 2Jr? At x = Jr? Graph the sum of the first four harmonics and observe how well it approximates f. 4. Show that the Fourier series for [(x) integer, is cos ax

=

2a sin aJr

Jr

1 ( 2a2

=

cos ax on [-Jr, Jrj, where a is not an

~

+ ~(-1)

cosnx )

11

a2

_

n2

.

Show csc aJr =

+ -2a Jr

2:: (-1) 00

n=1

5. Let f(x) = - ~ on -Jr < x < 0 and f(x) Fourier series for [ is

=

1 -.

n2

-

a2

~ on 0 :::: x :::: Jr. Show that the

2:: (2n -2 l)Jr sin(2n 00

11

1)x.

11=1

If SN(X) denotes the sum of the first N terms, sketch a graph of SI (x), S3(X), S7(X), and SIO(X) and compare with fex). Observe that the approximations overshoot fex) in a neighborhood of x = 0, and the overshoot is not improved regardless of how many terms are taken in the approximation. This overshoot behavior of Fourier series near a discontinuity is called the Gibbs

phenomenon.

108

3.4

3.

Orthogonal Expansions

Sturm-Liouville Problems

Let us revisit the basic idea of Fourier that we discussed in Section 3.1. We were able to represent the solution of the initial boundary value problem Ut

U(O, t)

=

(3.21)

0 < x < n, t > 0,

= U xx ,

=

u(n, t)

(3.22)

0 < x < n,

0,

(3.23)

t > 0,

u(x,O) = f(x),

as

=L 00

u(x, t)

(3.24)

cn e- n2t sin nx

n=l

for suitably chosen constants Cn . We note that this solution is an expansion in terms of orthogonal functions sin nx. We recall from Section 3.1 that orthogonality was essential in computing the Cn from the initial condition. How do we know that this particular set of orthogonal functions is the appropriate choice for this problem? In answering this question we now describe one of the fundamental methods, called the Fourier method or separation of variables, to solve PDEs on bounded domains. We will observe that every such PDE leads in a natural way to a boundary value problem for an ODE whose solutions are the appropriate orthogonal functions for the problem; this will mean that the solution to the PDE can be represented in terms of those orthogonal functions. We restrict our attention to (3.21)-(3.22). The separation of variables method consists in looking for solutions ofthe form of products, i.e., (3.25)

u(x, t) = y(x)g(t).

Substituting this into the PDE (3.21) and boundary conditions (3.22), we obtain y(x)g'(t)

= yl/(x)g(t) ,

y(O)g(t)

= 0,

y(n)g(t)

= 0,

which we can write as g'(t) get)

yl/(x)

-y(x) = -A,

yeO)

= yen) = O.

Here, -A is some unknown constant (we may set g' /g and yl/ /y equal to a common constant because the only way a function of t can equal a function of x for all t and x is ifboth are equal to a constant); the minus sign on Ais for convenience. Therefore, we obtain an ordinary differential equation for g, namely, g'(t)

=

-Ag(t) ,

3.4.

Sturm-Liouville Problems

109

and we obtain a boundary value problem for y, namely, -y"(x)

= AY(X) ,

0 < x < n,

(3.26)

yeO) = 0, yen) = o.

(3.27)

So, the PDE problem separated into two ODE problems, for get) and y(x). We can easily solve the g-problem; it has a solution get)

= e- At •

As it turns out (our goal is to discuss this in detail), the boundary value problem (3.26)-(3.27) for y will have a solution only when the constant A takes on the special values A = An = n Z , n = 1, 2, .... The corresponding solutions are y = Yn(x) = sin nx, n = 1,2, .... Then we will have g = gn(t) = e- n2t , n = 1,2, .... Thus we will have constructed infinitely many solutions U = un(x, t)

= gn(t)Yn(X) = e- n2t sin nx

that satisfy the PDE (3.21) and the boundary conditions (3.22). Now we can clearly see where the terms in the solution (3.24) arise. The orthogonal functions sin nx come from the solution of the ODE boundary value problem (3.26)-(3.27) for y. What we have just briefly described is the separation of variables method, and the procedure can be imitated on a large number of PDE models. Each will give rise to a boundary value problem for an ODE that has infinitely many solutions that form an orthogonal system. In this section we focus attention on the ODE boundary value problem, saving a full discussion of the separation of variables method to Chapter 4. When we separate variables, i.e., assume (3.25), for an equation of the form Ut

= (P(x)ux)x - q(x)u,

a < x < b,

(here, a and b are finite), we obtain an ODE for Y -(P(x)y')'

+ q(x)y =

AY,

t

> 0

= y(x) of the form

a < x < b,

(3.28)

where A is some constant (see Exercise 1). This ODE is called a SturmLiouville differential equation. The boundary conditions on the PDE will usually lead to boundary conditions on y(x) ofthe form aly(a)

+ azY'(a) = 0,

(3.29)

fhy(b)

+ flzy'(b) =

(3.30)

O.

Here the constants al and az are not both zero, and the constants ,Bland flz are not both zero (i.e., the boundary condition at an endpoint does not collapse). TWo special cases of the boundary conditions (3.29)-(3.30) are yea)

= 0,

y(b)

= 0 (Dirichlet conditions)

(3.31)

3. Orthogonal Expansions

110

and y'(a) = 0, y'(b) = 0

(Neumann conditions).

(3.32)

Ifp, p', and q are continuous functions on the interval [a, b] and p is never zero in [a, b], then the ODE boundary value problem (3.28)-(3.30) is called a regular Sturm-Liouville problem (SLP). Otherwise it is called singular. Singular problems usually arise when the interval [a, b] is infinite or when p(xo) = 0 for some Xo E [a, b]. SLPs are named after J.C.F. Sturm and J. Liouville, who studied such problems in the mid 1830s. A regular SLP will not have a solution for every value of the constant A. (Clearly, the zero function y(x) == 0 is always a trivial solution, but we have no interest in it). A value of A for which there is a nontrivial solution of (3.28)-(3.30) is called an eigenvalue, and the corresponding solution is called an eigenfunction; observe that any constant multiple of an eigenfunction gives another (not independent) eigenfunction. The interesting fact about regular SLPs is that they have an infinite number of eigenvalues, and the corresponding eigenfunctions form a complete, orthogonal set, which makes orthogonal expansions possible. This is a key idea in applied mathematics. Before introducing some general facts for SLPs, we study some examples. EXAMPLE

The SLP

-y"(x) = Ay(X) ,

0 < x < n,

yeO) = 0, yen) = 0,

(3.33) (3.34)

has eigenvalues A = An = n 2 and corresponding eigenfunctions y = yn(x) = sin nx, n = 1,2, .... Thus (3.33)-(3.34) has solution y = sin x when A = I, Y = sin 2x when A = 4, and so on. One way to find the eigenvalues and eigenfunctions is to go through a case argument by separately considering A = 0, A > 0, and A < O. We will prove later that A cannot be a complex number. We examine different cases because the solution of (3.33) has a different form depending on the sign of A. If A = 0, then the ODE (3.33) has the form y" = 0, whose general solution is the linear function y(x) = Ax + B. But yeO) = 0 implies B = 0, and yen) = 0 implies A = O. Thus we get only the trivial solution in this case and so A = 0 is not an eigenvalue. If A < 0, say for definiteness A = _k2, then the ODE (3.33) has the form y" - k2y = 0; this has general solution consisting of exponentials,

3.4.

111

Sturm-Liouville Problems

The boundary conditions force yeO) yen)

= A + B = 0, = Aekrr + Be- krr

= O.

It is easy to see that these two equations for A and B force A = B = 0, and therefore we obtain only the trivial solution in this case. Thus, there are no negative eigenvalues. Finally, let us consider the case A > 0, or A = 1 0, and it satisfies the boundary conditions, for only mild restrictions on the initial temperature distribution[(x) (e.g., [is a piecewise smooth function on 0 :::: x :::: z). More severe conditions are required on [ if one wishes to prove that the solution surface we obtained is continuous on 0 :::: x :::: z, t:::: O. Justification and verification can be found in more advanced texts on PDEs. Second, one can show that the solution we obtained is the only solution to the problem. Next, for the separation of variables method to be successful, the PDE and boundary conditions must both be linear and homogeneous. The method would not be expected to work for nonlinear equations. If we are faced with a problem having nonhomogeneous boundary conditions, then that problem must be transformed into one with homogeneous conditions; as we observed earlier, that transformation leads to a source term in the PDE. Later we shall learn how to deal with sources. Finally, once an infinite series representation for the solution is found, the work is not over if we want temperature profiles. Generally, we cannot sum an infinite series, and therefore we often take the first few terms as an approximation. Such approximations will be illustrated in an example below and in the exercises. The separation of variables procedure described above can be imitated for a large number of problems. Now we illustrate the method for the

4.

120

Partial Differential Equations on Bounded Domains

wave equation, but with less detail. Consider the problem (4.10) u(O, t) u(x, 0)

= 0, u(l, t) = 0, t > 0, = F(x) , Ut(x, O) = G(x) ,

(4.11) 0 < x < Z.

(4.12)

Assuming u(x, t) = y(x)g(t) and substituting into the PDE (4.10) yields y(x)g"(t)

= e2 y"(x)g(t),

or, upon dividing by e2 yg, g"(t)

y"(x)

e 2g(t)

y(x)

Setting each term equal to a constant -A, we obtain the two ODEs g"(t)

+ e2 Ag(t) = 0,

Next we substitute u(x, t) to obtain

-y"(x)

= Ay(X).

= y(x)g(t) into the boundary conditions (4.11)

y(O)g(t)

= 0,

y(l)g(t)

= 0,

which gives yeO) = 0 and y(Z) = O. Therefore, we are led to the boundary value problem -y"(x)

= AY(X) ,

0 < x < Z;

yeO)

= 0,

y(l) = O.

This is the Sturm-Liouville problem for y, and it is exactly the same as in the heat flow example above. The eigenvalues and eigenfunctions are (see (4.6)-(4.7)) yn(x)

= sin

mrx

-Z- , n

= 1,2 ....

(4.13)

Solving the ODE for get), nnet get) = C sin -Z-

+ D cos

nnet -Z-·

Therefore, we have constructed infinitely many product solutions of the PDE (4.10) having the form un(x, t)

=

nnet ( en sin -Z-

+ dn cos

nnet ) nnx -Z- sin -Z-,

n

= 1,2, ... ,

where we have replaced the arbitrary constants C and D by arbitrary constants en and d n depending on n. These product solutions un(x, t) represent modes of vibrations. The temporal part is periodic in time with period 2nll(nen); the spatial part is a "standing wave" with frequency nnll. These product solutions also satisfy the boundary conditions (4.11)

4.1. Separation of Variables

121

but do not satisfy the given initial conditions (4.12). So we form the linear combination U(X,

~ (mfct Cn sin -z-

t) = ~

en

and select the constants Thus we require

mfct ) + lin cos -z-

mfX -z-

(4.14)

and dn such that the initial conditions hold.

0) = F(x) =

U(X,

sin

nJl'X L dn sin -z. 00

n=l

The right side is the Fourier sine series of the function Fex) on the interval (0, D. Therefore, the coefficients dn are the Fourier coefficients given by dn =

211

1

0

nJl'X

F(x) sin -Z-dx,

n = 1,2, ....

(4.15)

Th apply the other initial condition Ut(x, 0) = Gex) we need to calculate the time derivative ofu(x, t). We obtain

L 00

Ut(x, t) =

n=l

ncJl' ( n J l ' c t

nJl'ct )

Cn cos -Z- - lin sin -Z-

-Z-

nJl'X

sin -Z-·

Thus Ut(x,O)

~ ncJl'

= G(x) = ~ n=l

n1l'X

-Z-Cn sin -Z-·

Again, the right side is the Fourier sine series of G on (0, D, so the coefficients of sin n7x , which are n~7I' cn, are the Fourier coefficients; that is, ncJl'

-Z-en =

or Cn

2

1

2 = -ncJl'

11 0

11 0

nJl'X

G(X) sin -Z-dx,

nJl'X

G(x) sin -Z- dx,

n

= 1,2, ... ,

n = 1,2, ....

(4.16)

Therefore, the solution of the initial boundary value problem is given by the infinite series (4.14) where the coefficients are given by (4.15) and (4.16). Computer algebra packages can be of tremendous help in developing the formulas found by the separation of variables method to obtain approximate solution surfaces.

4.

122

Partial Differential Equations on Bounded Domains

EXAMPLE

Consider the following initial boundary value problem for the heat equation: Ut

U(O, t)

=

U xx ,

0 < x < I, t > 0,

= u(I , t) = 0,

u(x,O) = 10x3 (1 - x),

> 0, 0 < x < 1.

The solution is given by the formulae (4.8) and (4 .9). In Maple the following sequence of commands defines the functionf(x), computes the Fourier coefficients en, and then finds the first six terms of an approximate solution. Here k = 1 and 1 = 1. f:=x -+ 10*x A 3*(1-x) i coef:=n -+ 2*int(f(x)*sin (n*Pi*x) ,x=O .. 1)

i

c:=seq(coef(n) ,n=1 .. 6) i

The approximate solution surface can be plotted for 0 command plot3d(u(x,t) ,x=O .. 1,t=0 .. 0.1)

:s t :s

0.1 via the

i

The temperature surface is shown in Figure 4.1. The initial temperature profile f(x) and the temperature profile at t = 0.1 are shown in Figure

no

.b

.. nl

n.

" 0'

I

Figure 4.1. The temperature surface.

Exercises

123

0.8

0.6

0.4

0.2

o0

0.2

0.4

x

Figure 4.2. Thmperature profiles at t = 0 and t = 0.1.

0.6

4.2, which was generated by the commands pl:=plot(f(x),x=O .. l) : p2:=plot(u(x,O.l) ,x=O .. l): display ( [pl,p2]); These graphs show how the bar cools down.

D

Exercises 1. In the heat flow problem (4.1)-(4.3) take k = 1, I = 77:, and [(x) = 0 if o < x < 77:12, [(x) 1 if77:12 < x < 77:. Find an infinite series representation

=

of the solution. Use the first four terms in the series to obtain an approximate solution, and on the same set of coordinate axes sketch several time snapshots of the approximate temperature distribution in the bar in order to show how the bar cools down. Estimate the error in these approximate distributions.

2. In the vibration problem (4.10)-(4.12) take c = 1, I = 77:, and F(x) = x if o < x < 77:/2, and F(x) = 77: - x if 77:12 < x < 77:, and take G(x) == O. Find an infinite series representation of the solution. Use the first four terms in the series to obtain an approximate solution, and on the same set of coordinate axes sketch several time snapshots of the wave. Can you estimate the error in these approximate wave forms? 3. Consider the pure boundary value problem for Laplace's equation given by

+ Uyy

= 0,

0
0,

Ut(x,O) = 0,

0 < x < I,

governs the displacement of a string immersed in a fluid. The string has unit length and is fixed at its ends; its initial displacement is f, and it has no initial velocity. The constant k is the damping constant. Use separation of variables to find the solution in the case k < 2nc.

4.2

Flux and Radiation Conditions

In the last section we solved the initial boundary value problem for the heat equation: Ut

= kuxx ,

0 < x < I, t > 0,

u(O, t) = u(l, t) = 0, u(x,O) = [(x),

t > 0,

0 < x < 1.

The boundary conditions in this problem, where we have fixed the temperature to be zero at the ends of the bar, are examples of Dirichlet, or fixed endpoint, conditions. Now we wish to consider two other types of boundary conditions; the first is a Neumann, or flux, condition, and the second is a Robin, or radiation, condition. An insulation-type flux condition at, say, x = 0 is a condition on the derivative of the form u;\{O, t) = 0, t > 0; because the flux is proportional to the temperature gradient (Fourier's heat flow law states flux = - Ku x, where K is the conductivity, see Section 1.4), the insulation condition requires that no heat flow across the boundary x = O. A radiation condition, on the other hand, is a specification of how heat radiates from the end of the bar, say at x = 0, into its environment, or how the end absorbs heat from its environment. Linear, homogeneous radiation conditions take the form -Kux(O, t) + bu(O, t) = 0, t > 0, where b is a constant. If b > 0, then the heat flux is negative, which means that heat is flowing from the bar into its surroundings (radiation); ifb < 0, then the flux is positive, and heat is flowing into the bar (absorption). A typical problem in heat conduction may have a combination of Dirichlet, insulation, and radiation boundary conditions. In other physical contexts, e.g., contaminant transport, biological diffusion, these three types of boundary conditions also play an important role.

4.2.

Flux aud Radiation Conditions

125

The method of eigenfunction expansions applies to problems with these different types of homogeneous boundary conditions. The technique is exactly the same as in the last section; the only change will be the form of the boundary conditions in the Sturm-Liouville problem that we obtain upon separating variables. EXAMPLE

Consider the diffusion problem Ut

Ux(O, t)

0 < x < I, t > 0,

= kuxx ,

= 0,

u(l, t)

= 0,

(4.17)

t > 0,

(4.18) (4.19)

0 < x < l,

u(x,O) = fCx),

where the left end is insulated and the right end is fixed. We assume that u(x, t) = y(x)g(t) and substitute this into the PDE to obtain, exactly as in Section 4.1, g'(t) yl/(x) = - - =-A kg(t) y(x) ,

where -A is the separation constant. This gives the two ODEs g'(t)

= -Akg(t),

-yl/(x)

= Ay(X).

Now substituting the product u = yg into the boundary conditions gives y' (O)g(t) = 0,

y(l)g(t) = O.

Since we do not want get) to be zero, we are forced to take y'(O) = 0 and y(l) = O. Therefore, we have obtained the Sturm-Liouville problem -yl/ = AY,

0 < x < l;

y'(O) = 0,

y(l) = O.

We solved this problem in Section 3.4 to obtain eigenvalues and eigenfunctions given by (the reader should review that calculation) An

= (2n: 1

y7:,

yn(X) = cos (2n :/)1rX ,

n = 0, 1,2 ....

By the results in Section 3.4 we know that the eigenfunctions yn(x) are orthogonal on 0 < x < l, which means that

1/

(2n

+ 1)1rX 1

(2m

+ 1)1rX ~.. _

I o 2 2 The solution to the time equation for get) is cos

cos

get)

(..U -

0 for m:j:. n.

= ce-J..kt,

and since there are infinitely many values of A, we have infinitely many such solutions,

4.

126

Partial Differential Equations on Bounded Domains

Summarizing, we have obtained infinitely many product solutions ofthe form

cne

-A kt n

cos

+ l)nx

(2n

'

2l

n=O,I,2 ....

These functions will solve the PDE (4.17) and the boundary conditions (4.18). Now we superimpose these to form u(x t) = ,

2: cn e00

Ankt

n=O

(2n

+ 1)nx

cos -'------'2l'

(4.20)

and we properly select the coefficients en such that u will satisfy the initial condition (4.19). 1b this end we have ~ (2n + l)nx u(x, 0) = f(x) = ~ en cos 2l . n=O

This equation is an expansion of the function f(x) in terms of the eigenfunctions Yn(x). We know from Section 3.2 that the coefficients en are given by the formula

t

1

en = IIYnl1 2 Jo f(x) cos

(2n

+ l)nx 2l

ax,

n

= 0, I, 2 ....

A straightforward integration (either by hand or using a computer algebra package) gives

llYn II

2

=

11 o

COS

2

(2n

+ l)nx ax = -, l 2l

2

n

= 0, I, 2 ....

Consequently, the solution of(4.17)-(4.19) is given by (4.20), where the coefficients en are given by 211 (2n + 1 )nx f(x) cos ax, l 0 2l

en = -

n = 0, I, 2 ....

o

Now we solve a problem with a radiation boundary condition. EXAMPLE

Consider the diffusion problem Ut

=

kuxx ,

u(O, t) = 0, u(x,O) =f(x),

(4.21 )

0 < x < I, t > 0,

u(1, t)

+ u x (I, t)

0 < x < I,

= 0,

t > 0,

( 4.22) ( 4.23)

where the left end is fixed and the right end satisfies a radiation condition. The separation of variables method dictates that we assume u(x, t) =

4.2.

Flux and Radiation Conditions

127

y(x)g(t). Substituting this into the PDE, we obtain, exactly as above, g'(t) yl/(x) - - = - - =-A kg(t) y(x) ,

where -A is the separation constant. This gives the two ODEs

= -Akg(t),

g'(t)

-yl/(x)

= Ay(X).

Substituting the product u = yg into the boundary conditions gives. y(O)g(t)

= 0,

y(l)g(t)

+ y'(l)g(t) = O.

Since get) =1= 0, we are forced to take yeO) = 0 and y(l) we are left with the Sturm-Liouville problem -yl/

= AY,

0 < x < 1;

yeO)

= 0,

gel)

+ y'(l) = 0,

and

+ y'(l) = O.

Now we determine the eigenvalues, i.e., the values of A for which this problem has a nontrivial solution. In the case A = 0 the solution of the ODE is y(x) = mx + b; the two boundary conditions force m = b = 0, so zero is not an eigenvalue. To show that there are no negative eigenvalues, we argue as in Section 3.4, using an energy argument. We multiply the ODE by y(x) and integrate from 0 to 1 to get

11

-yyl/ ax

=A

11

y 2ax.

On the left side we integrate by parts to obtain -(y(x)y'(x)) 16

+

11

y,2ax = A

11

y 2ax.

Because the integrands are positive, the integral on the right is positive, and the integral on the left is positive. The boundary term is nonnegative since -(y(x)y'(x)) 1& = -y(l)y'(l)

+ y(O)y'(O)

= y(li :::

o.

Therefore, A cannot be negative. For A > 0 the solution to the ODE is y( x) = A cos ,J).x + B sin ,J).x. The boundary condition yeO) = 0 forces A = 0, and the boundary condition y(l) + y'(l) = 0 forces the condition sin .Ji.. + .Ji.. cos.Ji.. = 0, or .Ji.. = - tan .Ji... This equation, where the variable A is tied up in a trigonometric function, cannot be solved analytically. But graphs of the functions ,J). and - tan,J). versus A (see Figure 4.3) show that the equation has infinitely many positive solutions AI, A2, ... , represented by the intersection points

4. Partial Differential Equations on Bounded Domains

128

80

40

100

Figure 4.3. Graphical representation of the eigenvalues of the SLP.

of the two graphs. These values An are the eigenvalues for the Sturm-Liouville problem, and the corresponding eigenfunctions are Yn = sin ~x,

n = I, 2, ....

Using a solving routine available in computer algebra programs (for example, f sol ve in Maple) we obtain the first three eigenvalues: Al = 4.115858, A2 = 24.13934, A3 = 63.65911. By the results in Section 3.4 we know that the eigenfunctions Yn(x) are orthogonal on 0 < x < I, which means that

t

Jo

sin ~x sin ~xdx = 0 for m =f. n.

The solution to the time equation is

gU)

= ce- Akt ,

and since there are infinitely many values of A, we have infinitely many such solutions,

gnU) = cne-Ankt. Summarizing, we have obtained infinitely many product solutions of the form

These functions will solve the PDE (4.21) and the boundary conditions (4.22). Now we superimpose these to form

= L cne- Ankt sin ~x, 00

u(x, t)

ii=1

(4.24)

4.2.

Flux and Radiation Conditions

and we select the coefficients condition (4.23). We have

129

Cn

such that u will satisfy the initial

= [(x) = L 00

u(x, 0)

Cn

sin ~x.

n=l

This equation is an expansion of the function ftx) in terms of the eigenfunctions Yn(x). We know from Section 3.2 that the coefficients Cn are given by the formula Cn

1 = -112

Ynll

or Cn

=

1

io sin

11

21

Fnxdx

0

= 1,2 ... ,

[(x) sin ~xdx, n

t [(x) sin ~xdx,

Jo

n = 1,2, ....

(4.25)

Consequently, the solution of (4.21 )-(4.23) is given by (4.24), where the coefficients Cn are given by (4.25). 0 Computer algebra packages can be extremely helpful in obtaining approximate solutions and their graphs. Let us use Maple to calculate a three-term approximation to the solution in the last example in the case k = 0.1 and the initial temperature is ftx) = 10x(1 - x). As noted above, the first eigenvalue laml can be determined by the commands sol1:=fsolve({sqrt (z)+tan(sqrt(z) )=O},{z},1..5 ): assign(soll)

i

Zi

laml:=z

In the f sol ve routine the expression 1..5 gives a range of values where a known solution lies. Additional eigenvalues can be obtained in a similar manner. The first three eigenvalues are laml=4.115858,lam2=24.13934, lam3=63.65911. The coefficient C1 is calculated via cl:=int(lO*x*(l-x)*sin(sqrt(laml)*x),x=O .. l)/ int(sin(sqrt(laml)*x)A2,x=O .. 1)i The coefficients c2 and c3 are calculated similarly using lam2 and lam3. We obtain c1

= 2.13285,

c2

= 1.040488,

c3

= -0.219788.

Then a three-term approximate solution is u(x, t) ~ 2.13285e-4.115858t sin J4.115858x

+ 1.040488e-24.13934t sin .J24.13934x - 0.219788e-63.65911t sin J63.65911x.

4.

130

Partial Differential Equations on Bounded Domains

"

o ,

"., 1

.,

-

0 -1

~_-• _ - - - 01

Figure 4.4. Approximate solution surface.

01 1

"

The solution surface is shown in Figure 4.4 . Notice that the approximate solution does not match the initial condition accurately at x = 1.

Exercises 1. Find an infinite series representation for the solution to the wave problem

U(O, t) = ux(l, t) = 0, U(X,O)

= [(x),

t > 0,

Ut(X,O)

=

°°

< x < 1.

Interpret this problem in the context of waves on a string. 2. Find an infinite-series representation for the solution to the equilibrium problem

° < x < a, ° < y < b, Ux(O, y) = UX(a, y) = 0, ° < y < b,

Uxx

+ Uyy = 0,

U(x, 0) = [(x),

u(x, b) = 0,

°

< x < a.

Interpret this problem in the context of steady heat flow. 3. Find an infinite-series representation for the solution to the heat absorptionradiation problem

Ux(O, t) -

° °

= Uxx , < x < 1, t > 0, aou(O, t) = 0, ux(l, t) + aj u(l, t) = 0, Ut

u(x,O)

= f(x),

< x < 1,

t > 0,

4.3.

131

Laplace's Equation

Figure 4.5.

Circular ring.

where ao < 0, a] > 0, and ao + a] > -aOa]. State why there is radiation at x = 1, absorption at x = 0, and the radiation greatly exceeds the absorption . Choose ao = -0.25 and a] = 4 and find the first four eigenvalues; if [(x) == x(1 - x), find an approximate solution to the problem and graph the approximate solution surface. 4. Consider a large, circular, tubular ring of circumference 21 that contains a chemical of concentration c(x, t) dissolved in water. Let x be the arc-length parameter with < x < 21. See Figure 4.5. If the concentration of the chemical is initially given by Co (x) , then c(x, t) satisfies the initial boundary value problem

°

Ct

== Dc,,,,,

c(O, t) = c(21, t), c(x,O) == [(x),

°

< x < 21, t > 0,

Cx(O, t) == cx (21, t),

°< x < 2l.

t > 0,

These boundary conditions are called periodic boundary conditions, and D is the diffusion constant. Apply the separation of variables method and show that the associated Sturm-Liouville problem has eigenvalues An = (mr/li for n == 0,1,2, . . . and eigenfunctions Yo(x) = 1, Yn(x) == An cos(n:n:xl l) + En sine n:n:xll) for n == 1, 2, .... Show that the concentration is given by c(x, t) ==

Ao

"2 +

00

L(An cos(n:n:xll)

•

2 2

+ En sm(n:n:xl l))e- n "

Dtll

2

n=O

and find formulae for the All and En.

4.3

Laplace's Equation

Now we investigate some problems associated with Laplace's equation on a bounded domain. In the first two sections of this chapter we have presented a general method, the Fourier method, to solve Laplace's equation on rectangles. Now we examine circular domains. Moreover, we shall study some general properties of harmonic functions (i.e., solutions to liu = 0) in both two and three dimensions. Suppose we know the temperature on the boundary of a circular, laminar plate of radius R. Can we find the equilibrium temperature inside

4. Partial Differential Equations on Bounded Domains

132

the plate? We will now find a particularly nice formula for the solution to this problem. When the region of interest is a circle, we might guess that polar coordinates are more appropriate than rectangular coordinates, so let us formulate the steady-state heat flow problem in polar coordinates r, (), where x = r cos () and y = r sin (). Then a circular plate of radius R can be represented simply as r :::: R with 0 :::: () :::: 2rr. The unknown temperature inside the plate is u = u(r, 8), and the given temperature on the boundary of the plate is u(R, 8) = fCe), where [ is a known function. We know from prior discussions (Section 1.8) that u must satisfy Laplace's equation D-u = 0 inside the plate. Therefore, upon representing the Laplacian Din polar coordinates (see Section 1.8), we have the following boundary value problem for u(r, 8): U rr

1

1

+ -r U y + "2 UIJ(} = 0, r u(R, 8)

0:::: r < R, 0:::: () :::: 2rr,

= [(8),

0:::: () :::: 2rr.

(4.26) (4.27)

Implicit are the periodic boundary conditions u(r, 0)

= u(r, 2rr),

Ue(r, 0)

= Ue(r, 2rr).

(4.28)

The separation of variables method works the same way as it did in rectangular coordinates. We assume the product solution u(r, 8) = y(r)g(8) and substitute into the PDE and boundary conditions. When the expression is substituted into the PDE (4.26), we obtain y"(r)g(())

1

+ -r y'(r)g(8) +

1

"2 y(r)g"(8) r

= O.

This can be written

-

r 2Y"(r) + ry'(r) y(r)

gil (e) =- = -A, g(8)

where -A is the separation constant. Therefore, we obtain two ordinary differential equations for y and g, namely, r 2y"(r)

+ ry'(r)

= -Ay(r), g"(8) = -Ag(8).

The periodic boundary conditions force g(O) = g(2rr) and g'(O) Therefore, we have the following problem for g(8): g"(8)

= -Ag(8);

g(O)

= g(2rr),

g'(O)

= g'(2rr).

= g'(2rr). (4.29)

The boundary conditions, which are called periodic, are not quite the type that we defined for a Sturm-Liouville problem. But many of the same results are true. First, it is clear that A = 0 is an eigenvalue with corresponding eigenfunction goer) = 1. Moreover, there are no negative eigenvalues; if A is

4.3.

Laplace's Equation

133

negative, then the ODE has exponential solutions, and exponential solutions cannot satisfy periodicity conditions (or an energy argument could be used). Therefore, let us assume A = p2 > O. The differential equation for g has general solution g(Y) = a cos pO

+ b sin pO.

The boundary conditions force (cos2np - l)a

(sin 2np)a

+ (1

+ (sin2np)b

= 0,

- cos 2np)b = O.

This is a system of two linear homogeneous equations for a and b. From a basic result in matrix algebra, we know that it has a nontrivial solution ifthe determinant ofthe coefficients is zero, i.e., (cos2np - 1)(1 - cos2np) - sin 2 2np = 0,

or, simplifying, cos 2np = O. This means that p =

../I must be a positive integer, that is, A = An = n 2 ,

n = 1,2, ....

Along with AO = 0, these are the eigenvalues of the problem (4.29). The eigenfunctions are go(O)

=

I, gn(O)

= an cos n(j + bn sin nO,

n

= 1,2, ....

Now let us solve the y-equation. Of course, we want bounded solutions. For A = 0 the only bounded solution is yo(r) = 1 (the other solution in this case is In r, which is unbounded). For A = n 2 the equation is r 2y"(r)

+ ry'(r) + n2y(r)

= 0,

which is a Cauchy-Euler equation (see the Appendix) with general solution (4.30)

We can set en = 0 because again we seek bounded solutions for 0 Thus, setting d n = 1 for all n, we have yn(r) = r n,

:s r :s R.

n = 1,2, ....

In summary, we have constructed solutions of the given boundary value problem of the form uo(r, 0) = constant = aol2, u = un(r,O) = rn(a n cos nO + bn sin nO), n = 1,2, .... Now, to satisfy the boundary condition at r = R we form the linear combination ao 00 (4.31) u(r, 0) = - + rn(a n cos nO + bn sin nY). 2

L

n=1

4. Partial Differential Equations on Bounded Domains

134

The boundary condition u(R, t!) fCt!)

=

a

= fCt!) then yields

+ L Rn(an cos ne + bn sin nt!), 00

20

n=l

which is the full-range Fourier series for f(t!) (see Section 3.3). Therefore, the coefficients are given by 1 an = rrRn

1

bn = rrRn

121r f(t!) cos ne de,

n

= 0, 1,2, ... ,

( 4.32)

121r fee) sin ne de,

n

= 1,2, ....

(4.33)

0

0

So the solution to the BVP (4.26)-(4.28) is given by (4.31) with the coefficients given by (4.32)-(4.33). As it turns out, this series solution can be manipulated to obtain a simple formula, called Poisson's integral formula. Let us substitute the coefficients given in formulae (4.32) and (4.33) (after changing the dummy integration variable from e to ¢) into the solution formula (4.31) to obtain

r21r f(¢)d¢ rn 121r f( ¢)(cos n¢ cos ne + sin n¢ sin ne)d¢ +L

u(r, t!) = Jo

00

n=l

-n

rrR

r

0

21r 1 = 2rr Jo f(¢)

( 1 + ~ (r)n R cos n(e 00

¢)

) d¢.

But the infinite sum in the integrand can be determined exactly as follows. Recalling the identity cos a = ~ (eia + e- i"), we can write 1

(r)n cosn(e - ¢) = 1 + ~ (r)n R .

+~ R 00

00

etn((j- 0 and n is a bounded region. Show that the only solution is u = O. (Hint: Use an energy argument. Multiply the PDE by u, integrate, and then use Green's first identity.) Use this result to show that solutions to the boundary value problem tJ.u - cu = g(x, y, z), n . grad u

+ au =

[(x,

(x, y, z) E

n,

y, z), on an,

must be unique. 7. Solve the boundary value problem tJ.u = 0,

u(R,

e)

= 4

r < R,

0:::::

e
0, w(O, t, r) = w(77:, t, r) = 0, t > 0,

Wt -

kwxx

w(x, 0, r) = rex, r),

0 < x
0,

Ut

t > 0,

u(O, t) = fU),

x > 0,

u(x, 0) = 0,

(4.69) (4.70) (4.71)

which is a well-posed problem. We are asking whether we can determine K from a single flux measurement

-Kux(O, to) = a,

(4.72)

where a is known. Using Laplace transforms we solved the direct problem (4.69)-(4.71) in Section 2.6. The solution is u(x, t) =

r j4KnU x

Jo

r)3

e-x2!(4K(t-r))f(r)dr.

(4.73)

It appears that the strategy should be to calculate the flux at (0, to) from

the solution formula (4.73). Indeed this is the case, but calculating the x derivative of u is not a simple matter. The straightforward approach of pulling a partial derivative a/ ax under the integral sign fails because one of the resulting improper integrals cannot be evaluated at x = 0; it does not exist. The reader should verify this statement. Therefore, we must be more clever in calculating u x . To this end we note that (4.73) can be written as u(x, t) = -2K

1 axa t

o

-

G(x, t - r)f(r)dr,

where G(x, t) is the heat kernel G(x, t) = _1_ e- X2 /(4Kt)

J4Kt

Therefore, ux(x, t) = -2K

= 2

1 a2 t

o

1 ara t

o

-2

ax

G(x, t - r)f(r)dr

- G(x, t - r)f( r)dr,

Exercises

155

since -G, = Gt = KGxx (the heat kernel satisfies the heat equation). Now we integrate by parts to obtain ux(x, t) =

-21

t

G(x, t - 1:)f'(1:)d1:.

The boundary terms generated by the integration by parts are zero. Consequently, we have -Kux(O, to) =

../K

l

to

['(1:)

o Jrr(to - 1:)

d1: = a.

This equation uniquely determines K and solves the parameter identification problem. For example, if fCt) = f3t, Le., the temperature is increased linearly, then the integral can be calculated exactly (Exercise 6), and we obtain rra z K= - - .

4f3 Zto

(4.74)

The exercises contain additional examples. Often, inverse problems, and in particular distributed parameter identification problems, lead to an integral equation for the unknown function. These equations are difficult to resolve, and stability is often a problem.

Exercises 1. In this exercise we examine the sensitivity of determining the growth constant Y from a measurement of the population in the example at the beginning ofthis section. Suppose two different measurements PI and P2 ofthe population are made at the same time t*. Each leads to growth constants YI and Yz. Show that if PI and pz are close, then YI and Y2 are close. In particular, show that iYI - Y2i ::: CiPI - P2i,

where C is a constant. Thus, this recovery strategy is stable in the sense that a small error in measuring the population will not drastically change the growth constant. 2. This problem deals with determining the unknown density of a nonhomogeneous vibrating string by observing one fundamental frequency and the corresponding mode of vibration. Consider a stretched string of unit length and unit tension that is fastened at both ends. The displacement is governed by (see Section l.5) p(X)Utt = Uxx , 0 < x < I, t > 0, U(O, t)

== u(1, t)

= O. t > 0,

4. Partial Differential Equations on Bounded Domains

156

where p(x) is the unknown density. If we assume a solution of the form u(x, t) = y(x)g(t), then we obtain a Sturm-Liouville problem _y"

= Ap(X)Y,

0 < x < 1;

yeO)

= y(l) = O.

Suppose we observe (say, using a strobe light) a fundamental frequency A = Af and the associated mode (eigenfunction) y = Yf(X). Show that the unknown density must satisfy the integral equation

1 1

y'(O) = _f_.

(1 - X)Yf(x)p(x)dx

a (Hint: Integrate the ODE from x = 0 to x Determine the density if it is a constant.

Af

= s and then from s = 0 to s = 1).

3. Consider heat flow in a rod oflength 7r and unit diffusivity whose ends are held at constant zero temperature and whose initial temperature is zero degrees. Suppose further that there is an external heat source ftx) supplying heat to the bar. Address the question of recovering the heat source ftx) from a temporal me.asurement of the temperature U(t) made at the midpoint of the rod. Proceed by formulating the appropriate equations and show that ftx) must be the solution to the integral equation UCt) =

f 11
0,

t > 0,

(4.75) (4.76)

0 < x < 1.

(4.77)

The first step to discretize the region of space-time where we want to obtain a solution. In this case the region is 0 S x S I, 0 S t S T. We have put a bound on time because in practice we only solve a problem up until a finite time. Discretizing means defining a lattice of points in this space-time region by Xj

= jh,

tn

= nk,

j

= 0, I, ... , J;

n

= 0, I, ... , N,

4.8. Finite Difference Methods

159

wr---+---+---,r--~---+--~--~

~+1r---+---+---,r--~---+--~--~

~r---+---~---r---+---+--~--~

xJ

x

Figure 4.7. The discrete lattice and the computational atom for the diffusion equation.

where the fixed numbers hand k are the spatial and temporal step sizes, respectively. Here, h = III and k = liN. The integer 1 is the number of subintervals in 0 ::::: x ::::: I, and N is the number of time steps to be taken. Figure 4.7 shows the lattice of points (the lattice points are also called grid points or nodes). At each node (Xj, tn) of the lattice we seek an approximation, which we call to the exact value u(Xj, tn) of the solution. Note that the superscript n refers to time, and the subscript j refers to space. We regard as a two-dimensional array, or matrix, where n is a row index and j is a column index. 1b obtain equations for the we replace the partial derivatives in the PDE by their difference approximations. From calculus we know

Up,

Up

Up,

Ut(Xj, tn)

~

u(Xj, tn+l )

-

k

u(Xj, tn)

and uxx(Xj, tn)

U(Xj-I, tn) - 2u(xj, tn) h2

~

+ u(Xj+l' tn)

The first of these formulas is the usual forward difference approximation for a first derivative, and the second is a difference approximation for a second derivative; the latter follows from a Taylor approximation exactly as derived in Section l.8. So the PDE (4.75) at the point (Xj, tn ) is replaced by the difference equation

u n+ l

}

_

k

Un

} =

~~l

-

2~n

+ ~~l

h2

or, upon solving for ~n+l,

~n+l

=

Ur + :2 (~~l -

2~n + ~~l).

(4.78)

Observe that this equation relates the approximate values of the solution at the four points (Xj-l, tn), (Xj, tn), (Xj-l, tn), (Xj, tn+l). These four points form the computational atom for the difference scheme (see Figure 4.7).

4. Partial Differential Equations on Bounded Domains

160

The differ~nce equation (4.78) gives the approxi~at~ solution at the node (Xj, tn+1) in terms of approximations at three earlier nodes. Now we see how to fill up the lattice with approximate values. We know the values at t = 0 from the initial condition (4.77). That is, we know ~o

= fCxj),

j

= 0, I, ... ,J.

From the boundary conditions (4:76) we also know U~

= 0,

UF

= 0,

n

= 1,2, ... , N.

The difference formula (4.78) can now be applied at all the interior lattice points, beginning with the values at the t = ql~vel, to compute the values at the t = t1 level, then using those to compute the values at the t = t2 level, and so on. Therefore, using the difference equation, we cim march forward in time, continually updating the previous temperature profile. Think of filling out the array ~n, row by row. A finite difference scheme, or algorithm, like (4.78) is called an explicit scheme because it permits the explicit calculation of an approximation at the next time step in terms of values at a previous time step. Because of the error (called the truncation error) that is present in (4.78) due to replacing derivatives by differences, the scheme is not always accurate. If the time step k is too large, then highly inaccurate results will be obtained. It can be shown that we must have the stability condition k 1 0, u(x,O) = fcx), Ut(x,O) = g(x),

(4.80)

u(O,t)

0

=::

x

=::

1.

(4.81)

4.8.

Finite Difference Methods

161

01

Figure 4.8a. A graph of the temperature surface when [(x) =

x 3 (1 - x). Set the discretization: J:=10: h:=l/J: k:=0.004: N:=50: r:=k/h-2:

Define the array and set the initial condition: u:=array(l .. J+1,1 .. N+1) : f:=x- >x-3*(l-x): for j from 1 to J+1 do U[j,1]:=f«j-1))*h): od:

Set the Dirichlet boundary conditions: for n from 1 to N+1 do U[l,n] :=0.0: U[J+1,n] :=0.0: od:

Fill up the array using the difference equation: for n from 1 to N do: for j from 2 to J do: U [j ,n+1l : =U [j ,nl +r* (U [j -1, nl - 2*U [j ,n] +U [j +1, nl ) : od: od:

Figure 4.8b. Maple program listing to solve (4.75)-(4.77) when [(x) = x 3 (1 - x).

Plot the solution array: with(linalg): with(plots): matrixplot(U,labels=['x', 't', 'u'l);

Using the same lattice as for the diffusion problem, we can approximate the PDE by

un- 1 _ ]

2U n ]

+ U n+ 1 ]

which gives Uj

n+l

=

n

2~ - Uj

n-l

+ (Ck)2 h (Ujn_ 1 -

n

2~

+ Ujn+1)·

The computational atom is shown in Figure 4.9. 1b compute the value at the (n + 1)st time step we now need values from the previous two time steps. Therefore, to start the marching scheme we require values from the first two (to = 0 and tl = k) time rows. The t = 0 row is given by the initial condition u(x, 0) = [(x); so we have UjO = [(Xj),

j = 0,1, ... ,J.

4. Partial Differential Equations on Bounded Domains

162

The t = t} time row can be computed using the initial velocity condition Ut(x, 0) = g(x), which we approximate by U1

_

] k

UO

]

= g(Xj),

j

= 1,2, ... ,J -

1.

Now we have all the ingredients to start the scheme and march forward in time to fill up the lattice with values. Again there is a stability condition, namely, c::S:

h

"k.

This is called the Courant-Friedrichs-Lax (CFL) condition. Physically, it states that the speed hlk of which lattice points are being calculated must exceed the speed c of which waves are propagated in the system. Thus the calculated lattice points will contain all of the information about the wave. In the exercises we ask the reader to compute the solution to a wave-like problem. Now we solve a BVP for Laplace's equation in the unit square with given values on the boundary of the square. This Dirichlet problem is a pure boundary value problem, and we should not expect a marchingtype scheme to work; there is no guarantee that when we march from one boundary in the problem to another we will reach the specified values at that boundary. This means that we must solve for the approximate values at all of the lattice points simultaneously. Tb illustrate the procedure we consider the problem U xx ·

+ Uyy = 0,

0