An Introduction to Partial Differential Equations

TL;DR

This paper introduces partial differential equations (PDEs) from an applied mathematics perspective, focusing on physical context, methods, and applications, suitable for advanced students with basic programming and differential equations background.

Contribution

It provides an accessible overview of PDEs emphasizing methods and applications, bridging physical intuition with analytical and numerical techniques.

Findings

Introduces basic PDEs like heat, wave, and Laplace equations.

Discusses analytical and numerical solution methods.

Provides educational material for advanced students.

Abstract

The field of partial differential equations (PDEs) is vast in size and diversity. The basic reason for this is that essentially all fundamental laws of physics are formulated in terms of PDEs. In addition, approximations to these fundamental laws, that form a patchwork of mathematical models covering the range from the smallest to the largest observable space-time scales, are also formulated in terms of PDEs. The diverse applications of PDEs in science and technology testify to the flexibility and expressiveness of the language of PDEs, but it also makes it a hard topic to teach right. Exactly because of the diversity of applications, there are just so many different points of view when it comes to PDEs. These lecture notes view the subject through the lens of applied mathematics. From this point of view, the physical context for basic equations like the heat equation, the wave equation and the Laplace equation are introduced early on, and the focus of the lecture notes are on methods, rather than precise mathematical definitions and proofs. With respect to methods, both analytical and numerical approaches are discussed. These lecture notes has been succesfully used as the text for a master class in partial differential equations for several years. The students attending this class are assumed to have previously attended a standard beginners class in ordinary differential equations and a standard beginners class in numerical methods. It is also assumed that they are familiar with programming at the level of a beginners class in informatics at the university level.

Equations

$$u_x+u_y$$

$$u_x+y{u}_y$$

$$(u_x{)}^2+(u_y{)}^2$$

$$u_x\equiv {\partial }_{x}u\equiv \frac{\partial u}{\partial x}.$$

$$F(x,y,u(x,y),u_x(x,y),u_y(x,y))=0,$$

$$(\text{or short}F(x,y,u,u_x,u_y)=0)$$

$$F(x,y,u,u_x,u_y,u_{xx},u_{xy},u_{yy})=0.$$

$$u_{xx}+u_{yy}$$

$$u_{tt}-u_{xx}$$

$$u_t+u{u}_x-u_{xx}$$

$$\mathcal{L}={\partial }_x+{\partial }_y.$$

$$\mathcal{L}u=0.$$

$$\mathcal{L}(u+v)$$

$$={\partial }_{x}u+{\partial }_{x}v+{\partial }_{y}u+{\partial }_{y}v$$

$$={\partial }_{x}u+{\partial }_{y}u+{\partial }_{x}u+{\partial }_{y}v$$

$$=\mathcal{L}u+\mathcal{L}v,$$

$$\mathcal{L}(cu)$$

$$=c{\partial }_{x}u+c{\partial }_{y}u$$

$$=c({\partial }_{x}u+{\partial }_{y}u)$$

$$=c\mathcal{L}(u).$$

$$\mathcal{L}(u)=g,$$

$$\mathcal{L}(u)=({\partial }_{x}u{)}^2+({\partial }_{y}u{)}^2.$$

$$\mathcal{L}(u)=1.$$

$$\mathcal{L}(cu)$$

$$=c^2(({\partial }_{x}u{)}^2+({\partial }_{y}u{)}^2)$$

$$=c^2\mathcal{L}(u).$$

$$({\partial }_{x}u{)}^2+({\partial }_{y}u{)}^2=1,$$

$$a(x,y)u_x+b(x,y)u_y+c(x,y)u=d(x,y).$$

$$\mathcal{L}u$$

$$\mathcal{L}v$$

$$\mathcal{L}(w)$$

$$\stackrel{\text{1}}{=}\mathcal{L}(au)+\mathcal{L}(bu)$$

$$\stackrel{\text{2}}{=}a\mathcal{L}(u)+b\mathcal{L}(u)$$

$$=0.$$

$$w=\sum _{i=1}^n{a}_i{u}_i,a_i\in \mathbb{R}(\text{or}\mathbb{C}),$$

$$y=c_1{y}_1+\cdots +c_n{y}_n,$$

$$u_{xx}=0.$$

$$u_x=f(y),$$

Full text

1. 1 Introduction

2. 2 First notions

3. 3 PDEs as mathematical models of physical systems

4. 3.1 A traffic model

5. 3.2 Diffusion

6. 3.3 Heat conduction

7. 3.4 A random walk

8. 3.5 The vibrating string

9. 4 Initial and boundary conditions

10. 4.1 Well posed Problems

11. 5 Numerical methods

12. 5.1 The finite difference method for the heat equation

13. 5.2 The Crank-Nicholson scheme

14. 5.2.1 Finite difference method for the heat equation: Neumann boundary conditions

15. 5.3 Finite difference method for the wave equation

16. 5.4 Finite difference method for the Laplace equation

17. 5.4.1 Iterations in $\mathbb{R}^{n}$

18. 5.4.2 Gauss-Jacobi

19. 5.4.3 Gauss-Seidel

20. 6 First order PDEs

21. 6.1 Higher order equations and first order systems

22. 6.2 First order linear scalar PDEs

23. 6.3 First order Quasilinear scalar PDEs

24. 7 Classification of equations: Characteristics

25. 7.1 Canonical forms for equations of hyperbolic type

26. 7.2 Canonical forms for equations of parabolic type

27. 7.3 Canonical forms for equations of elliptic type

28. 7.4 Characteristic curves

29. 7.5 Characteristic curves again: A different point of view

30. 7.6 Stability theory, energy conservation and dispersion

31. 8 Adjoint differential operators

32. 9 Initial and boundary value problems in bounded domains

33. 9.1 Separation of variables

34. 9.2 Self-adjoint and positive operators

35. 10 The Sturm-Liouville problem and Fourier series

36. 10.1 Properties of eigenvalues and eigenfunctions

37. 10.1.1 Trigonometric eigenfunctions

38. 10.1.2 Bessel eigenfunctions

39. 10.1.3 Legendre polynomial eigenfunctions

40. 10.2 Series solutions to initial boundary value problems

41. 10.2.1 Vibrations of a string with fixed endpoints

42. 10.2.2 Heat conduction in a finite rod

43. 10.2.3 The Laplace equation in a rectangle

44. 10.3 Inhomogenous equations - Duhamel’s principle

45. 10.3.1 The inhomogeneous wave equation

46. 10.3.2 The inhomogeneous heat equation

47. 10.4 The finite Fourier transform

48. 10.4.1 Hyperbolic equations

49. 10.4.2 The Poisson equation in a disk

50. 10.5 Nonlinear stability theory

51. 10.5.1 The infinite case

52. 10.5.2 The finite case

53. 11 Integral transforms

54. 11.1 The 1D Fourier transform

55. 11.1.1 A boundary value problem for an ODE

56. 11.1.2 The Cauchy problem for the heat equation

57. 11.1.3 The Cauchy problem for the wave equation

58. 11.1.4 The Laplace equation in a half-plane

59. 11.2 The Fourier sine- and cosine transform

60. 11.2.1 The heat equation on a semi-infinite interval

61. 11.3 The Laplace transform

62. 11.3.1 The wave equation on a half interval

63. 11.3.2 The heat equation on a finite interval

64. 11.4 Group velocity 1

65. 11.5 Method of stationary phase

66. 11.6 Group velocity 2

67. 12 Projects

68. 12.1 Project1

69. 12.2 Project2

70. 12.3 Project3

1 Introduction

The field of partial differential equations (PDEs) is vast in size and diversity. The basic reason for this is that essentially all fundamental laws of physics are formulated in terms of PDEs. In addition, approximations to these fundamental laws, that form a patchwork of mathematical models covering the range from the smallest to the largest observable space-time scales, are also formulated in terms of PDEs. The diverse applications of PDEs in science and technology testify to the flexibility and expressiveness of the language of PDEs, but it also makes it a hard topic to teach right. Exactly because of the diversity of applications, there are just so many different points of view when it comes to PDEs. These lecture notes view the subject through the lens of applied mathematics. From this point of view, the physical context for basic equations like the heat equation, the wave equation and the Laplace equation are introduced early on, and the focus of the lecture notes are on methods, rather than precise mathematical definitions and proofs. With respect to methods, both analytical and numerical approaches are discussed.

These lecture notes has been succesfully used as the text for a master class in partial differential equations for several years. The students attending this class are assumed to have previously attended a standard beginners class in ordinary differential equations and a standard beginners class in numerical methods. It is also assumed that they are familiar with programming at the level of a beginners class in informatics at the university level.

While writing these lecture notes we have been influenced by the writings of some of the many authors that previously have written testbooks on partial differential equations [3, 6, 5, 8, 7, 2, 4]. However, for the students that these lecture notes are aimed at, books like [3],[6],[5],[8] and [7] are much too advanced either mathematically [3],[6],[8],[8] or technically [5],[7]. The books [2],[4] would be accessible to the students we have in mind, but in [2] there is too much focus on precise mathematical statements and proofs and less on methods, in particular numerical methods. The book [4] is a better match than [2] for the students these lecure notes has been written for, but it is a little superficial and does not reach far enought. With respect to reach and choise of topics, the book “Partial Differential Equations of Applied Mathematics”, written by Erich Zauderer[1] would be a perfect match for the students we have mind. However, the students taking the class covered by these notes does not have any previous exposure to PDEs, and the book [1] is for the most part too difficult for them to follow. The book is also very voluminous and contain far to much material for a class covering one semester with five hours of lecture each week. The lecture notes for the most part follow the structure of [1], but simplify the language and makes a selection of topics that can be covered in the time available. The lecture notes deviate from [1] at several points, in particular in the section covering physical modeling examples, integral transforms and in the treatment of numerical methods. After mastering these lecture notes, the students should be able to use [1] as a reference for more advanced methods that are not covered by the notes. In order for students to pass the master class these notes has been written for, the students must complete the three projects included in the last section of the lecture notes and pass an oral exam where they must defend their work and also answer questions based on the content of these lecture notes.

2 First notions

A partial differential equation(PDE), is an equation involving one or more functions of two or more variables, and their partial derivatives, up to some finite order. Here are some examples

[TABLE]

where

[TABLE]

The highest derivative that occur in the equation is the order of the equation. All equations (1,2,3) are of order one.

A PDE is scalar if it involves only one unknown function. In this course we will mostly work with scalar PDEs of two independent variables. This is hard enough as we will see!

Most important theoretical and numerical issues can be illustrated using such equations. It is also a fact that many physical systems can be modeled by equations of this type. We will also for the most part restrict our attention on equations of order two or less. This is the type of equations that most frequently occur in modeling situations. This fact is linked to the use of second order Taylor expansions and to the use of Newton’s equations in one form or another.

The most general first order scalar PDE, in two independent variables, is of the form.

[TABLE]

[TABLE]

where $F$ is some function of five variables.

Let $D\in\mathbb{R}^{2}$ be some region in the xy-plane. Then $u(x,y)$, defined on $D$, is a solution of equation (4) if $u_{x}$, $u_{y}$ exists and are continuous on $D$ and if $u$ satisfy the equation for all $(x,y)\in D$. PDEs typically have many different solutions. The different solutions can be defined on different, even disjoint domains.

We will later relax this notion of solutions somewhat. In more advanced textbooks what we call a solution here is called a classical solution.

The most general second order scalar PDE of two independent variables is of the form

[TABLE]

A solution in $D\in\mathbb{R}^{2}$ is now a function where all second order partial derivation are continuous and where the function satisfy equation (5) for all $(x,y)\in D$.

Here are some examples

[TABLE]

Equations (1),(2),(6) and (7) are linear equations whereas equations (3) and (8) are not linear. Equations that are not linear are said to be nonlinear. The importance of the linearity/nonlinearity distinction in the theory of PDEs can hardly be overstated. We have already seen the importance of this notion in the theory of ODEs.

The precise meaning of linearity is best stated using the language of linear algebra. Let us consider equation (1).

Define an operator $\mathscr{L}$ by

[TABLE]

Equation (1) can then be written as

[TABLE]

The operator $\mathscr{L}$ has the following important properties:

1)

For all differentiable functions $u$,$v$

[TABLE]

2)

For all differentiable functions $u$ and $c\in\mathbb{R}$ (or $\mathbb{C})$

[TABLE]

Properties $1$ and $2$ show us that $\mathscr{L}$ is a linear differential operator.

Definition 1**.**

A scalar PDE is linear if it can be written in the form

[TABLE]

for some linear differential operator $\mathscr{L}$ and given function $g$. The equation is homogeneous if $g=0$ and inhomogeneous if $g\neq 0$.

Let us consider equation (3) from page 1. Define a differential operator $\mathscr{L}$ by

[TABLE]

Then the equation can be written as

[TABLE]

Note that

[TABLE]

$\mathscr{L}$ fails property $2$ and is thus not a linear operator (It also fails property $1$). Thus the equation

[TABLE]

is nonlinear. It is not hard to see that, basically, a PDE is linear if it can be written as a first order polynomial in $u$ and its derivatives. All other equations are nonlinear.

The most general first order linear scalar PDE of two independent variables is of the form

[TABLE]

The most important property of homogeneous linear equations is that they satisfy the superposition principle.

Let $u,v$ be solutions of a linear homogeneous scalar PDE. Thus

[TABLE]

where $\mathscr{L}$ is the linear differential operator defining the equation. Let $w=au+bv$ where $a,b\in\mathbb{R}$ (or $\mathbb{C}$). Then

[TABLE]

Thus $w$ is also a solution for all choices of $a$ and $b$. This is the superposition principle. We have already seen this principle at work in the theory of ODEs.

Generalizing the argument just given, we obviously have that if $u_{1},\dots,u_{n}$ are solutions then

[TABLE]

is also a solution for all choices of the constants $\{a_{i}\}$.

In the theory of ODEs we found that any solution of a homogeneous linear ODE of order $n$ is of the form

[TABLE]

where $y_{1},\dots,y_{n}$ are $n$ linearly independent solutions of the ODE and $c_{1},\dots,c_{n}$ are some constants. Thus the general solution has $n$ free constants. For linear PDEs, the situation is not that simple.

Let us consider the PDE

[TABLE]

Integrating once with respect to $x$ we get

[TABLE]

where $f(y)$ is an arbitrary function of $y$. Integrating one more time we get the general solution

[TABLE]

where $g(y)$ is another arbitrary function of $y$. Thus the general solution of equation (9) depends on two arbitrary functions. The corresponding ODE

[TABLE]

has a general solution $y=ax+b$ that depends on two arbitrary constants. The equation (9) is very special but the conclusion we can derive from it is very general. (But not universal! )

“ The general solution of PDEs depends on one or more arbitrary functions.”

This fact makes it much harder to work with PDEs than ODEs.

Let us consider another example that will be important later in the course. The PDE that we want to solve is

[TABLE]

Integrating with respect to $x$ we get

[TABLE]

where $f(y)$ is arbitrary. Integrating equation (11) with respect to $y$ gives us the general solution to equation (10) in the form

[TABLE]

where $G(y)=\int f(y)\mathrm{d}y$ and $F(x)$ are arbitrary functions of $y$ and $x$.

3 PDEs as mathematical models of physical systems

Methods for solving PDEs is the focus of this course. Deriving approximate description of physical systems in terms of PDEs is less of a focus and is in fact better done in specialized courses in physics, chemistry, finance etc.

It is however important to have some insight into the physical modeling context for some of the main types of equations we discuss in the course. This insight will make the methods we introduce to solve PDEs more natural and easy to understand.

3.1 A traffic model

We will introduce a model for cars moving on a road. In order for the complexity not to get out of hand, we start by making some assumptions about the constructions of the road.

1. i)

There is only one lane, so that all cars move in the same direction (no passing is allowed). 2. ii)

There are no intersections or off/on ramps. Thus cars can not enter or leave the road.

Because of i) and ii), we can represent the road by a line that we coordinatize by a variable $x$. Let

[TABLE]

From the definition of $\Delta N$ and $f$ and assumption ii) we get the identity

[TABLE]

Equation 12 express the conservation of the number of cars on the road.

Both $\Delta N$ and $f$ are integer valued and thus change by whole numbers. After all, they are counting variables! Define the density of cars, $n(x,t)$, by

[TABLE]

We now make the fundamental continuum assumption: $n(x,t)$ and f(x,t) are smooth functions. This is a reasonable approximation if the number of cars in each interval $\left[x,x+\Delta x\right]$ is very large and change slowly when $x$ and $t$ vary. $f(x,t)$ is also assumed to be slowly varying. (Imagine observing a crowded road from a great height, say from a helicopter). Rewrite equation (12) as

[TABLE]

Taylor expanding to first order we get

[TABLE]

Letting $\Delta t$ and $\Delta x$ approach zero, we get a PDE

[TABLE]

This type of PDE is called a conservation law. It express the conservation of cars in differential form. The notion of conservation laws is fundamental in the mathematical description of nature. Many of the most fundamental mathematical models of nature express conservation laws.

Note that equation (14) is a single first order PDE for two unknown functions, $n$ and $f$. In order to have a determinate mathematical model we must have the same number of equations and unknowns.

Completing the model involves making further physical assumptions. We will assume that the flux, $f$, of cars at $(x,t)$ depends on $x$ and $t$ only through the density of cars, $n(x,t)$. Thus

[TABLE]

where $\alpha$ is some function whose precise specification requires further modeling assumptions.

As noted above, the quantity $f$ in this model is called a flux. Such quantities are very common in application and typically measures the flow of some quantity through a point, surface or volume.

In equation (15), we assume that the flux of cars through a point only depends on the density of cars. This is certainly not a universal law, but it is a reasonable assumption (is it?). The chain rule now gives

[TABLE]

Thus our model of the traffic flow is given by the PDE

[TABLE]

This equation describes transport of cars along a road and is an example of a transport equation. Transport equations occur often as mathematical models because conserved quantities are common in our description of the world.

The model is still not fully specified because the function $c(n)$ is arbitrary at this point. If we assume that $c(n)$ has a convergent power series expansion, we have

[TABLE]

If the density of cars is sufficiently low we can use the approximation

[TABLE]

and our transport equation is a linear equation

[TABLE]

For somewhat higher densities we might have to use the more precise approximation

[TABLE]

which gives us a nonlinear transport equation

[TABLE]

For the linear transport equation (19), let us look for solutions of the form

[TABLE]

The chain rule gives

[TABLE]

Thus

[TABLE]

is a solution of equation (19) for any function $\varphi$. Let us assume that the density of cars is known at one time, say at $t=0$. Thus

[TABLE]

where $f$ is a given function. Using equation (22) and (21) we have

[TABLE]

[TABLE]

[TABLE]

[TABLE]

[TABLE]

Thus the function

[TABLE]

is a solution of

[TABLE]

(24) is an example of an initial value problem for a PDE. Initial value problems for PDEs are often, for historical reasons, called Cauchy problems.

Note that at this point we have not shown that (23) is the only solution to the Cauchy problem (24). Proving that there are no other solutions is obviously important because a deterministic model (and this is!) should give us exactly one solution, not several.

The solution (23) describe an interesting behavior that is common to very many physical systems. It describe a wave. From (23), it is evident that

[TABLE]

Thus we get the values of $n$ at time $t$ by translating the values of $n$ at $t=0$ along the $x$-axis.

We are used to observing waves when looking at the surface of the sea, but might not have expected to find it in the description of traffic on a one-lane road.

Traffic waves are in fact common and play no small part in the occurrence of traffic jams.

3.2 Diffusion

Consider a narrow pipe oriented along the $x$-axis containing water that is not moving. Assume there is a material substance, a dye for example, of mass density $u(x,t)$ dissolved in the water. Here we assume that the mass density is uniform across the pipe so that the density only depends on $x$ which by choice is the coordinate along the length of the pipe.

The total mass of the dye between $x_{1}$ and $x_{2}$ is

[TABLE]

If the pipe is not leaking, the mass inside the section of the pipe between $x_{1}$ and $x_{2}$ can only change by dye entering or leaving through the points $x_{1}$ and $x_{2}$. This is the law of conservation of mass which is the first great pillar of classical physics.

Let $f(x,t)$ be the amount of mass flowing through a point $x$ at time $t$ pr.unit time. This is the mass flux. By convention $f$ positive means that mass is flowing to the right.

Mass conservation then gives the balance equation

[TABLE]

Fick’s law states that we have the following relation between mass flux and mass density.

[TABLE]

Thus mass flow from high density areas towards low density areas.

Fick’s law is only a phenomenological approximate relation based on the empirical observations that mass flux is for the most part a function of $(x,t)$, thorugh the density gradient, thus $f(x,t)=f(u_{x}(x,t))$. Assuming that the flux is a smooth function of the density gradient, we can represent the flux as a power series in the density gradient

[TABLE]

Since there is no mass flux if the density gradient is zero we must have $f_{0}=0$. Furthermore, since mass flows from high density domains to low density domains we must have $f_{1}<0$. We can therefore write $f_{1}=-k$, where $k>0$. Clearly, if the density gradient is large enough, the third term in the power series 27 will be significant and must be included. Like for the traffic model, this will lead to a nonlinear equation instead of the linear equation that will be the result of truncating the power series at the second term. In general we will assume that the parameter $k$ in Fick’s law depends on $x$. This corresponds to the assumption that the relation between flux and density gradient depends on position. This is not uncommon. Using Fick’s law we can write the mass balance equation as

[TABLE]

[TABLE]

[TABLE]

[TABLE]

[TABLE]

[TABLE]

[TABLE]

This is the 1D diffusion equation. By using similar modeling in 2D and 3D, we get the 2D and 3D diffusion equations

[TABLE]

If the mass flow coefficient, $k$, is constant in space, the diffusion equation takes the form

[TABLE]

These equations and their close cousins describe a huge range of physical, chemical, biological and economic phenomena.

If the mass flow is stationary(time independent) and the mass flow coefficient independent of space, we get the simplified equation

[TABLE]

This is the Laplace equation and occur in the description of many stationary phenomena.

3.3 Heat conduction

Let us consider a narrow bar of some material, say a metal. We assume that the temperature of the bar, $T$, only depends on the coordinate along the bar which we call $x$.

Let the temperature of the bar at a point $x$ at time $t$ be $T(x,t)$. Let $e(x,t)$ be the energy density at $x$ at time $t$. The total energy in the bar between $x_{1}$ and $x_{2}$ is

[TABLE]

The second great pillar of classical physics is the law of conservation of energy. Assuming no energy is gained or lost along the length of the bar, the conservation of energy gives us the balance equation

[TABLE]

where $f_{i}=f_{i}(x_{i},t)$ and $f(x,t)$ is the energy flux at $x$ at time $t$. Here we use here the convention that positive energy flux means that energy flow to the right. In order to get a model for the bar, we must relate $e(x,t)$ and $f(x,t)$ to $T(x,t)$.

For many materials we have the following approximate identity

[TABLE]

where $\rho$ is the mass density of the bar and $C$ is the heat capacity at constant volume. For (32) to hold the temperature must not vary too quickly. Energy is always observed to flow from high to low temperature regions. If the temperature gradient is not too large, Fourier’s law hold for the energy flux

[TABLE]

The energy balance equation then become

[TABLE]

[TABLE]

[TABLE]

[TABLE]

[TABLE]

[TABLE]

[TABLE]

If the heat conduction coefficient, $q$, does not depend on $x$, we get

[TABLE]

Again we get the diffusion equation. In this context the equation is called the equation of heat conduction. Energy flow in 2D and 3D similarly leads to 2D and 3D version of the heat conduction equations

[TABLE]

and for stationary heat flow, $T_{t}=0$, we get the Laplace equation

[TABLE]

The wave equation, diffusion equation and Laplace equations and their close cousins will be the main focus in this course.

In the next example we will show that the diffusion equation also arise from a seemingly different context. This example is the beginning of a huge application that have ramifications in all areas of science, technology, economy, etc. I am talking about stochastic processes.

3.4 A random walk

Let $x$ be the position of a particle restricted to moving along a line. We assume that the particle, if left alone, will not move. At random intervals we assume that the particle experience a collision that makes it shift its position by a small amount $\delta$. The collision is equally likely to occur from the left or the right. Thus the shifts in position caused by a collision is equally likely to be $-\delta$ or $+\delta$. The collisions occur at random times, but let us assume that during some small time period, $\tau$, a collision will occur for sure.

Let $P(x,t)$ be the probability of finding the particle at a point $x$ at time $t$. We must then have

[TABLE]

This is just saying that probability is conserved. From this point of view (37) is a balance equation derived from a conservation law just like equations (11), (38), (25) and (31). Using Taylor expansions

[TABLE]

we get from equation (37)

[TABLE]

Thus by truncating the expansions we get

[TABLE]

This is again the diffusion equation. A more general and careful derivation of the random walk equation can be found in section 1.1 in [1].

3.5 The vibrating string

Let us consider a piece of string of length $l$ that is fixed at $x=0$ and $x=l$. We make the following modeling assumption

1. i)

The string can only move in a single plane, and in this plane it can only move transversally.

Using i) the state of the string can be described by a single function $u(x,t)$ measuring the height of the deformed string over the $x$-axis at time $t$. Negative values of $u$ means that the string is below the $x$-axis.

1. ii)

Any real string has a finite thickness, and bending of the string result in a tension (negative pressure) on the outer edge of the string and a compression(pressure) on the inner edge of the bend.

We model the string as infinitely thin and thus assume no compression force. The tension, $\mathbf{T}$, in a deformed string is always tangential to the string and assumed to be of constant length and time independent

[TABLE]

1. iii)

The string has a mass density given by $\rho$ which is independent of position and also time independent.

Note:

ii) and iii) are not totally independent because a typical cause of nonuniform tension, $T$, is that the string has a variable material composition and this usually imply varying mass density.

Let us consider a piece of string above the interval $[x_{1},x_{2}]$. The total vertical momentum (horizontal momentum is by assumption i) equal to zero), $P(t)$, in this piece of string at time $t$ is

[TABLE]

Recall that momentum is conserved in classical (nonrelativistic) physics. This is the third great pillar of classical physics. Note that with the $y$-axis oriented like in figure 5, a piece of the string moving vertically up gives a positive contribution to the total momentum $P(t)$. Recall also that force is by definition change in momentum per unit time. Thus, if we assume there are no forces acting on the string inside the interval $[x_{1},x_{2}]$, then momentum inside the interval can only change through the action of the tension force at the bounding points $x_{1}$ and $x_{2}$. If we orient the coordinated axes as in figure 5, a positive vertical component, $f_{1}$, of the tension force at $x_{1}$, will lead to an increase the total momentum $p(t)$. A similar convention apply for $f_{2}$ at $x_{2}$.

Conservation of momentum for the string is then expressed by the following balance equation

[TABLE]

Our last modelling assumption is

1. iv)

The string perform small oscillations

[TABLE]

Applying iv) gives us

[TABLE]

Thus our fundamental momentum balance equation is well approximated by

[TABLE]

[TABLE]

[TABLE]

[TABLE]

[TABLE]

[TABLE]

[TABLE]

[TABLE]

[TABLE]

This is the wave equation. Similar modeling of small vibrations of membranes (2D) and solids (3D) give wave equations

[TABLE]

Wave equations describe small vibrations and wave phenomena in all sorts of physical systems, from stars to molecules.

4 Initial and boundary conditions

PDEs typically have a very large number of solutions (free functions, remember?). We usually select a unique solution by specifying extra conditions that the solutions need to satisfy. Picking the right conditions is part of the modeling and will be guided by the physical system that is under consideration.

There are typically two types of such conditions:

1. i)

initial conditions 2. ii)

boundary conditions

An initial condition consists of specifying the state of the system at some particular time $t_{0}$. For the diffusion equation, we would specify

[TABLE]

for some given function $\varphi$. This would for example describe the distribution of dye at time $t_{0}$. For the wave equation, we need two initial conditions since the equation is of second order in time (same as for ODEs)

[TABLE]

The modelling typically leads us to consider the PDE on some domain $D$ (bounded or unbounded).

For the vibrating string it is natural to assume that the string has a finite length, $l$, so $D=[0,l]$. It is then intuitively obvious that in order to get a unique behavior we must specify what the conditions on the string are at the endpoints $x=0$ and $x=l$. After all, a clamped string certainly behave different from a string that is loose at one or both ends! Here are a couple of natural conditions for the string:

1. i)

$u(0,t)=0,u(l,t)=0$ , 2. ii)

$u(0,t)=g(t),u_{x}(l,t)=b(t)$.

For case i) the string is clamped at both ends and for case ii) the string is moved in a specified way at $x=0$, and is acted upon by a force $Tb(t)$ at the right end, $x=l$.

The three types of boundary conditions that most often occur are

1. D)

$u$ specified at the boundary. 2. N)

$\partial_{\mathbf{n}}u$ specified at the boundary. 3. R)

$\partial_{\mathbf{n}}u+au$ specified at the boundary.

The first is a Dirichlet condition, the second a Neumann condition and the third a Robin condition. Here $\mathbf{n}$ is a normal to the curve (2D) or surface (3D) that defines the boundary of our domain and $a$ is a function defined on the boundary.

Let us consider the temperature profile of a piece of material covering a domain $D$. We want to model heat conduction in this material.

1. i)

Assume that the material is submerged in a large body of water kept at a temperature $T=0$. The problem we must solve is then

[TABLE]

where $\varphi(\mathbf{x})$ is the initial temperature profile. 2. ii)

If the material is perfectly isolated (no heat loss), we must solve the problem

[TABLE]

The general morale here is that the choice of initial and boundary conditions are part of the modelling process just like the equation themselves.

4.1 Well posed Problems

Let us assume that the modelling of some physical system has produced a PDE, initial conditions and/or boundary conditions.

How can we know whether or not the model is any good? After all, the model is always a caricature of the real system. There are always many simplifying assumptions and approximations involved.

We could of course try to solve it and compare solutions to actual observation of the system. However, if this should even make sense to try, there are some basic requirements that must be satisfied.

1. i)

Existence: There is at least one solution that satisfy all the requirements of the model. 2. ii)

Uniqueness: There exists no more than one solution that satisfy all the requirements of the model. 3. iii)

Stability. The unique solution depends in a stable manner on the data of the problem: Small changes in data (boundary conditions, initial conditions, parameters in equations, etc, etc) lead to small changes in the solution.

Producing a formal mathematical proof that a given model is well posed can be a hard problem and belongs to the area of pure mathematics. In applied mathematics we derive the models using our best physical intuition and sound methods of approximation and assume as a “working hypothesis” that they are well posed. If they are not, the models will eventually break down and in this way let us know. The way they break down is often a hint we can use to modify them such that the modified models are well posed.

We thus use an “empirical” approach to the question of well-posedness in applied mathematics. We will return to this topic later in the class.

5 Numerical methods

In most realistic modelling, the equations and/or the geometry of the domains are too complicated for solution “by hand” to be possible. In the “old days” this would mean that the problems were unsolvable. Today, the ever increasing memory and power of computers offer a way to solve previously unsolvable problems. Both exact methods and approximate numerical methods are available, often combined in easy to use systems, like Matlab, Mathematica, Maple etc.

In this section we will discuss numerical methods. The use of numerical methods to solve PDEs does not make more theoretical methods obsolete. On the contrary, heavy use of mathematics is often required to reformulate the PDEs and write them in a form that is well suited for the architecture of the machine we will be running it on, which could be serial, shared memory, cluster etc.

Solving PDEs numerically has two sources of error.

Truncation error: This is an error that arises when the partial differential equation is approximated by some (multidimensional) difference equation. Continuous space and time can not be represented in a computer and as a consequence there will always be truncation errors. The challenge is to make them as small as possible and at the same time stay within the available computational resources. 2. 2)

Rounding errors: Any computer memory is finite and therefore only a finite set of real numbers can be represented. After any mathematical operation on the available numbers we must select the “closest possible” available number to represent the result. In this way we accumulate errors after every operation. In order for the accumulated errors not to swamp the result we are trying to compute, we must be careful about how the mathematical operations are performed and how many operations we perform.

1) and 2) apply to all possible numerical methods. Heavy use of mathematical methods are often necessary to build confidence that some numerical solution is an accurate representation of the true mathematical solution. The accuracy of numerical solutions are often tested by running them against known exact solutions to special initial and/or boundary conditions and to “stripped down” versions of the equations.

The simples truncation method to describe is the method of finite differences. This is the only truncation method we will discuss in this class. We will introduce the method by applying it to the 1D heat equation, the 1D wave equation and the 2D Laplace equation

5.1 The finite difference method for the heat equation

Consider the following initial/boundary value problem for the heat equation

[TABLE]

We introduce a uniform spacetime grid $(x_{i},t_{n})$, where

[TABLE]

where $x_{0}=0,x_{N}=1$ and thus $h=\frac{1}{N}$. The solution will only be computed at the grid points.

Evaluating the equation and the initial/boundary conditions at the grid points we find

[TABLE]

where $u_{i}^{n}=u(x_{i},t_{n})$.

We now must find expressions for $u_{t}$ and $u_{xx}$, using only function values at the grid points. Our main tool for this is Taylors theorem.

[TABLE]

Using the Taylor expansions we easily get the following expressions

[TABLE]

Since $t_{n}\pm k=t_{n\pm 1}$, we get the following three approximations to the first derivative with respect time

[TABLE]

The second derivative in space is also approximated using Taylor theorem

[TABLE]

From these two expressions we get

[TABLE]

Thus we get the following approximation to the second derivative in space

[TABLE]

Using the forward difference in time and center difference in space in (41) gives

[TABLE]

Define $s=\frac{k}{h^{2}}$. Then the problem we must solve is

[TABLE]

This is now a (partial) difference equation.

Let us use the initial condition

[TABLE]

The exact solution is found by separation of variables which we will discuss later in this class.

[TABLE]

In order to compute the numerical solution, we must choose values for the spacetime discretization parameters $h$ and $k$. In the figures 13 and 14, we plot the initial condition, the exact solution and the numerical solution for $h=0.05$, thus $N=50$, and $k=h^{2}s=0.0004s$ for $s=0.49$ and $0.51$. We use 100 iterations in $n$.

It is clear from the figure 14 that the numerical solution for $s=0.51$ is very bad. This is an example of a numerical instability. Numerical instability is a constant danger when we compute numerical solutions to PDEs.

For a simple equation like the heat equation we can understand this instability, why it appears and how to protect against it.

Let us try to solve the difference equation in (44) by separating variables. We thus seek a solution of the form

[TABLE]

[TABLE]

[TABLE]

[TABLE]

[TABLE]

Both sides must be equal to the same constant $\xi$. This gives us two difference equations coupled through the separation constant $\xi$.

[TABLE]

The solution of equation 1) is simply

[TABLE]

The second equation is a boundary value problem

[TABLE]

[TABLE]

We look for solutions of the form $X_{j}=z^{j}$. This implies that

[TABLE]

Thus (45) implies that

[TABLE]

[TABLE]

[TABLE]

This equation has two solutions $z_{1}$ and $z_{2}$ which we in the following assume to be different. (Show that the case when $z_{1}=z_{2}$ does not lead to a solution of the boundary value problem (45).)

The general solution of (45) is thus

[TABLE]

The general solution must satisfy the boundary conditions. This implies that

[TABLE]

Nontrivial solutions exists only if

[TABLE]

This condition can not be satisfied if $z_{1}\neq z_{2}$ are real. Since equation (46) has real coefficient, we must have

[TABLE]

But equation (46) $\Rightarrow$ $z_{1}z_{2}=1$ $\Leftrightarrow$ $\rho^{2}=1$ $\Leftrightarrow$ $\rho=1$. Thus $z_{1}=e^{i\theta}$, $z_{2}=e^{-i\theta}$. The condition (47) then implies that

[TABLE]

The corresponding solutions of the boundary value problem (45) are

[TABLE]

A basis of solutions can be choosen to be

[TABLE]

It is easy to verify that

[TABLE]

[TABLE]

[TABLE]

The vectors $\{X_{j}^{(r)}\}_{r=1}^{N-1}$ are thus linearly independent and indeed form a basis because $\mathbb{R}^{N-1}$ has dimension $N-1$.

Each vector $X_{j}^{(r)}$ gives a separated solution.

[TABLE]

to the boundary value problem

[TABLE]

where we have defined

[TABLE]

These solutions give exponential growth if $|\xi_{r}|>1$ for some $r$. Thus the requirement for stability is

[TABLE]

If this holds, the numerical scheme is stable.

Now

[TABLE]

[TABLE]

[TABLE]

Therefore $-1\leqslant\xi^{r}\leqslant 1$ if

[TABLE]

[TABLE]

[TABLE]

This is the stability condition indicated by the numerical test illustrated in figures (13) and (14).

5.2 The Crank-Nicholson scheme

It is possible to avoid the stability condition by modifying the simple scheme discussed earlier. Define

[TABLE]

and let $0\leqslant Q\leqslant 1$. For each such number, we define a “$Q$-scheme”

[TABLE]

When $Q=0$, we get the explicit scheme we have discussed. For $Q>0$ the scheme is implicit; At each step we must solve a linear system.

We analyze the stability of the scheme by inserting a separated solution of the form

[TABLE]

1. Note:

This separated solution will not satisfy Dirichlet boundary conditions but periodic boundary conditions, $u_{j+N}^{n}=u_{j}^{n}$ In this type of stability analysis we disregard the actual boundary conditions of the problem, using instead periodic boundary conditions. A numerical scheme is Von Neumann Stable if

[TABLE]

There are other notions of stability where the actual boundary conditions are taken into account, but we will not discuss them in this course. Usually Von Neumann Stability gives a good indicator of stability, also when the actual boundary conditions are used.

If we insert (49) into the $Q$-scheme we get

[TABLE]

Evidently $\xi_{r}\leqslant 1$. The condition $\xi_{r}\eqslantgtr-1$ requires that

[TABLE]

This is always true for $1-2Q\leqslant 0$. Thus if $\frac{1}{2}\leqslant Q\leqslant 1$ the numerical scheme is stable for all choices of grid parameters $k$ and $h$. The price we pay is that the scheme is implicit.

For $\frac{1}{2}\leqslant Q\leqslant 1$ we say that the $Q$-scheme is unconditionally stable. For $Q=\frac{1}{2}$, the $Q$-scheme is called the Crank-Nicolson scheme.

Observe that when $N$ is large, the angles

[TABLE]

form a very dense set of points in the angular interval $[0,2\pi]$. Thus, in the argument above we might for simplicity take $\theta$ to be a continuous variable in the interval $[0,2\pi]$. We will therefore in the following investigate Von Neuman stability by inserting

[TABLE]

where $\eta=e^{i\theta}$, into the numerical scheme of interest. Verify that using this simplified approach leads to the same stability conditions for the two schemes (44) and (48).

5.2.1 Finite difference method for the heat equation: Neumann boundary conditions

Let us consider the heat equation as before but now with Neumann boundary conditions.

[TABLE]

One possible discretization for the boundary conditions is to use forward difference at $x=0$ and backward difference at $x=1$.

[TABLE]

These approximations are only of order $h$. The space derivative in the equation was approximated to second order $h^{2}$. We get second order also at the boundary by using center difference. However we then need to extend our grid by “ghost points” $x_{-1}$ and $x_{N+1}$

[TABLE]

Thus our numerical scheme is

[TABLE]

Introduction of ghost points is the standard way to treat derivative boundary conditions for 1D,2D and 3D problems.

5.3 Finite difference method for the wave equation

We consider the following initial and boundary value problem for the wave equation

[TABLE]

on the real axis. We introduce a uniform spacetime grid $(x_{j},t_{n})$ where

[TABLE]

where $x_{0}=0,x_{N}=L$ and thus $h=\frac{L}{N}$.

The wave equation is required to hold at the grid points only.

[TABLE]

The partial derivatives are approximated by centered differences (order $h^{2}$, $k^{2}$).

[TABLE]

[TABLE]

[TABLE]

where $s=\frac{k^{2}}{h^{2}}$. The boundary conditions evaluated on the grid simply become $u_{0}^{n}=0,u_{N}^{n}=0$.

Note that the difference equation (50) is of order 2 in $n$ and $j$. In order to solve the difference equation and compute the solution at time level $n+1$, we need to know the values on time level $n$ and $n-1$. Thus to compute $u_{j}^{2}$, we need $u_{j}^{1}$ and $u_{j}^{0}$. These values are provided by the initial conditions

[TABLE]

Evaluated on the grid these are

[TABLE]

We now need a difference approximation for $u_{t}$ at timelevel $n=0$. Since the scheme (50) is of order 2 in $t$ and $x$, we should use a 2-order approximation for $u_{t}$. If we use a 1-order approximation for $u_{t}$ in the initial condition the whole numerical scheme is of order 1. Using a centered difference in time, we get

[TABLE]

We have introduced a “ghost point” $t_{-1}$ in time with associated $u$ value $u_{j}^{-1}$. If we insert $n=0$ in (50), we get a second identity involving $u_{j}^{-1}$

[TABLE]

We can use (51) and (52) to eliminate the ghost point and find in conjunction with $u_{j}^{0}=\varphi_{j}$ that

[TABLE]

These are the initial values for the scheme (50). We now test the numerical scheme using initial conditions

[TABLE]

The exact solution

[TABLE]

to this problem is found using a method that we will discuss later in this class. For the test we use $L=20$ and chose the parameter $a$ so large that $\varphi(x)$ to a good approximation satisfy the boundary conditions for the wave equation. In figures (17) and (18) we plot the initial condition, the exact solution and the numerical solution on the interval $(0,L)$. The space-time grid parameters are $h=0.1$, thus $N=200$, and $k=h\sqrt{s}=0.0004\sqrt{s}$ for $s=0.9$ in figure 17 and $s=1.1$ in figure 18. We use 60 iterations in $n$.

It appears as if the scheme is stable for $s\lesssim 1$ and unstable for $s>1$. Let us prove this using the Von Neumann Approach.

[TABLE]

Inserting this into the scheme we get

[TABLE]

[TABLE]

[TABLE]

Define $p=s(\cos\theta-1)$. Then the equation can be written as

[TABLE]

The roots are

[TABLE]

Clearly $p\leqslant 0$. If $p<-2$ then $p^{2}+2p=p(p+2)>0$ and there are two real roots. One of these must be $<-1\Rightarrow|\xi|>1\Rightarrow$ scheme unstable. If $p>-2$ then $p^{2}+2p<0$ and we have two complex conjugate roots

[TABLE]

These two roots satisfy

[TABLE]

we therefore have stability if $p>-2$. It is easy to verify that the scheme is also stable for $p=-2$. Taking into account the definition of $p$ we thus have stability if

[TABLE]

[TABLE]

This holds for all $\theta$ if $s\leqslant 1$ because clearly $\frac{2}{1-\cos\theta}\eqslantgtr 1$. Thus we have stability if

[TABLE]

which is what our numerical test indicated.

5.4 Finite difference method for the Laplace equation

Let us consider the problem

[TABLE]

We solve the problem by introducing a grid

[TABLE]

where $h=\frac{1}{N},p=\frac{1}{M}$.

Using center difference for $u_{xx}$ and $u_{yy}$, we get the scheme

[TABLE]

This scheme is 2-order in $h$ and $p$. Because of the symmetry there is no reason to choose $h$ and $p$ differently. For the case $h=p$, which corresponds to $M=N$, we get

[TABLE]

Observe that the value in the center is simply the mean of the value of all closest neighbors. The stencil of the scheme i shown in figure (20).

There is no iteration process in the scheme (53) (There is no time!). It is a linear system of variables for the $(N-1)(N-1)$ interior gridpoints. We are thus looking at a linear system of equations

[TABLE]

where $\mathbf{b}$ comes from the boundary values. Observe however, that the linear system (53) is not actually written in the form (54). In order to apply linear algebra to solve (53), we must write (53) in the form (54). The problem that arise at this point is that the interior gridpoints in (53) are not ordered, whereas the vector components in (54) are ordered. We must thus choose an ordering of the gridpoints. Each choice will lead to a different structure for the matrix A (for instance which elements are zero and which are nonzero). Since the efficiency of iterative methods depends on the structure of the matrix, this is a point where we can be clever (or stupid!).

For the problem at hand it is natural to order the interiour gridpoints lexiographical.

Thus we define the vector $\mathbf{x}$ as

[TABLE]

The $\mathbf{b}$ vector is obtained from the values of $u$ at the boundary. For example when $j=1$, we have

[TABLE]

The boundary condition gives

[TABLE]

$i=1,$

[TABLE]

[TABLE]

[TABLE]

$1<i<N-1,$

[TABLE]

[TABLE]

[TABLE]

$i=N-1,$

[TABLE]

[TABLE]

[TABLE]

Thus the matrix for the system has a structure that starts off as

[TABLE]

By writing down a larger part of the linear system we eventually realize that the matric of the system has the form of a $(N-1)\times(N-1)$ block matrix

[TABLE]

where $I$ is the $(N-1)(N-1)$ identity and where $B$ is the following $(N-1)\times(N-1)$ tri-diagonal matrix

[TABLE]

We observe that most elements in the matrices $A$ and $B$ are zero, we say the matrices are sparse. The structural placing of the nonzero elements represents the sparsity structure of the matrix. Finite difference methods for PDEs usually leads to sparse matrises with a simple sparsity structure.

When $N$ is large, as it often is, the linear system (54) represents an enormous computational challenge.

In 3D with 1000 points in each direction, we get a linear system with $10^{9}$ independent variables and the matrix consists of $10^{18}$ entries!

There are two broad classes of methods available for solving large linear systems of $P$ equations in $P$ unknowns.

Direct methods: These are variants of elementary Gaussian elimination. If the matrix is sparse (most elements are zero) direct methods can solve large linear systems fast. For dense matrices, the complexity of Gaussian elimination $(O(P^{3}))$ makes it impossible to solve more than few thousand equations in reasonable time. 2. 2)

Iterative methods: These methods try to find approximate solutions to (53) by defining an iteration

[TABLE]

that converges to a solution of the linear system of interest. Since each multiplication cost at most $O(P^{2})$ operations, fast convergence($n\ll P$) will give an approximate solution to (54) in $O(P^{2})$ operations, even for a dense matrix.

The area of iterative solution methods for linear systems is vast and we can not possibly give an overview of the available methods and issues in a introductory class on PDEs. However, it is appropriate to discuss some of the central ideas.

5.4.1 Iterations in $\mathbb{R}^{n}$

Let $\mathbf{u}^{0}\in\mathbb{R}^{n}$ be a given vector in $\mathbb{R}^{n}$ and consider the following iteration

[TABLE]

where $M$ is a $n\times n$ matrix and $\mathbf{c}$ is a fixed vector. Let us assume that $\mathbf{u}^{l}\to\mathbf{u}^{*}$ when $l\to\infty$. Then we have

[TABLE]

[TABLE]

[TABLE]

[TABLE]

[TABLE]

Thus $\mathbf{u}^{*}$ solves a linear system. The idea behind iterative solution of linear system can be described as follows:

Find a matrix $M$ and a vector $\mathbf{c}$ such that

1. i)

$A\mathbf{u}=f,$

[TABLE]

[TABLE] 2. ii)

The iteration

[TABLE]

converges.

The condition i) can be realized in many ways and one seeks choises $M$ and $\mathbf{c}$ such that the convergence is fast.

In order get anywhere, it is at this point evident that we need to be more precise about convergence in $\mathbb{R}^{n}$. This we do by introducing a vector norm on $\mathbb{R}^{n}$ that measure distance between vectors. We will also need a corresponding norm on $n\times n$ matrices.

A $n\times n$ matrix is a linear operator on $\mathbb{R}^{n}$. Let us assume that $\mathbb{R}^{n}$ comes with a norm $\lVert\cdot\rVert$ defined. There are many such norms in used for $\mathbb{R}^{n}$. The following three norms are frequently used

[TABLE]

Given a norm, $\lVert\cdot\rVert$, on $\mathbb{R}^{n}$ there is a corresponding norm on linear operators (matrices) defined by

[TABLE]

For matrix norms corresponding to vector norms we have the following two important inequalities.

[TABLE]

[TABLE]

For $\lVert\cdot\rVert_{1}$ and $\lVert\cdot\rVert_{\infty}$, one can show that

[TABLE]

For any matrix, define the spectral radius $\rho(A)$ by

[TABLE]

where $\lambda_{s}$ are the eigenvalues of $A$. For the $\lVert\cdot\rVert_{2}$ norm we have the following identity

[TABLE]

where $A^{*}$ is the adjoint of $A$. For the special case when $A$ is selfadjoint, $A=A^{*}$, we have

Property 1**.**

Let $A=A^{*}$. Then

[TABLE]

Proof.

Let $\lambda$ be an eigenvalue of $A$ with associated eigenvector $u$. Then we have

[TABLE]

[TABLE]

[TABLE]

[TABLE]

[TABLE]

Let us now consider the iteration

[TABLE]

and let us use the $\lVert\cdot\rVert_{2}$ norm on $\mathbb{R}^{n}$ with the associated norm $\lVert\cdot\rVert_{2}$ on $n\times n$ matrices. The iteration will then converge to $\mathbf{u}^{*}$ in $\mathbb{R}^{n}$ iff

[TABLE]

But

[TABLE]

[TABLE]

By induction we get

[TABLE]

Let us assume that $M$ is symmetric, then

[TABLE]

This gives us the following sufficient condition for convergence.

Theorem 1**.**

Let $M$ be symmetric and assume that for all eigenvalues $\lambda$ for $M$ we have

[TABLE]

Then the iteration

[TABLE]

converges to the unique solution of

[TABLE]

for all initial values $\mathbf{u}^{0}$.

Note that even though all initial values converges to the solution, the number of iterations we need to get a predetermined accuracy will depend on the choice of $\mathbf{u}^{0}$.

We are now ready to describe the two simplest iteration methods for solving $A\mathbf{u}=\mathbf{b}$.

5.4.2 Gauss-Jacobi

Any matrix can be written in the form

[TABLE]

where $L$ is lower triangle, $U$ is upper triangular and $D$ is diagonal. Using this decomposition we have

[TABLE]

[TABLE]

[TABLE]

[TABLE]

[TABLE]

Assume that $D$ is invertible (no zeroes on the diagonal). Then we get

[TABLE]

where $M_{J}=-D^{-1}(L+U)$, $\mathbf{c}=D^{-1}\mathbf{b}$. This gives the Gauss-Jacobi iteration.

[TABLE]

In component form the Gauss-Jacobi iteration can be written in the form

[TABLE]

Observe that

[TABLE]

[TABLE]

[TABLE]

[TABLE]

[TABLE]

[TABLE]

[TABLE]

Let $\lambda$ be an eigenvalue of $A$ with corresponding eigenvector $v$. Then we have

[TABLE]

This means that

[TABLE]

In order to find the spectrum of $A$, we will start from a property of the spectrum for block-matrices.

Let $C$ be a matrix on the form

[TABLE]

where $C_{ij}$ are $n\times n$ matrices. Thus $C$ is a $(nm)\times(nm)$ matrix. Let us assume that there is a common set of eigenvectors for the matrices $C_{ij}$. Denote this set by $\{v^{k}\}^{n}_{k=1}$. Thus we have

[TABLE]

where $\lambda_{ij}^{k}$ is the eigenvalue of the matrix $C_{ij}$ corresponding to the eigenvector $v^{k}$.

Let $\alpha_{1},\dots,\alpha_{m}$ be real numbers. We will seek eigenvectors of $C$ on the form

[TABLE]

Then $w^{k}$ is an eigenvector of $C$ with associated eigenvalue $\mu^{k}$ only if

[TABLE]

[TABLE]

[TABLE]

The eigenvalues of $C$ are thus given by the eigenvalues of the matrix

[TABLE]

Let us now return to the convergence problem for Gauss-Jacobi.

$A$ is a $(N-1)\times(N-1)$ block matrix given by

[TABLE]

where $B$ is a $(N-1)\times(N-1)$ tri-diagonal matrix given by

[TABLE]

$B$ obviously commute with $I$ and the zero matrix and they therefore has common eigenvectors. Since all vectors in $\mathbb{R}^{N-1}$ are eigenvectors of the identity matrix and the zero matrix, with eigenvalues 1 and zero, respectively, in order to set up the $\Lambda$ matrix we only need to find the eigenvalues of $B$.

Let $D$ be a $n\times n$ tri-diagonal matrix of the form

[TABLE]

Let $\lambda$ be an eigenvalue with corresponding eigenvector $\mathbf{v}=(v_{1},\dots,v_{n})$ Thus

[TABLE]

[TABLE]

[TABLE]

and $v_{0}=v_{n+1}=0$. This is a boundary-value problem for a difference equation of order 2. We write the equation on the form

[TABLE]

and look for solutions of the form

[TABLE]

Thus we get the algebraic condition

[TABLE]

[TABLE]

[TABLE]

If $r_{+}\neq r_{-}$, the general solution to the difference equation is

[TABLE]

the boundary conditions give

[TABLE]

nontrival solutions exits only if

[TABLE]

But  $r_{+}r_{-}=1\Rightarrow r_{-}=\frac{1}{r_{+}}$ and therefore from equation (5.4.2)

[TABLE]

Furthermore, from equation (56) we have

[TABLE]

The eigenvalues of the tri-diagonal matrix $D$ are thus

[TABLE]

Returning to our analysis of the matrix $A$, we can now assemble our matrix $\Lambda^{k}$

[TABLE]

where $\lambda_{k}=-4+2\cos\left(\frac{\pi k}{N}\right),\quad k=1,\dots,N-1$. $\Lambda^{k}$ is also tri-diagonal and the previous theory for tri-diagonal matrices give us the eigenvalues

[TABLE]

[TABLE]

[TABLE]

or using standard trigonometric identities

[TABLE]

and finally the spectrum of $M_{J}$ in the Gauss-Jacobi iteration method is from (55)

[TABLE]

Using Theorem 1, the stability condition for Gauss-Jacobi is

[TABLE]

[TABLE]

[TABLE]

Since both $x_{k},y_{l}\in(0,1)$, it is evident that the condition holds and that we have convergence. We get the largest eigenvalue by choosing $k=l=1$. The spectral radius can then for large $N$ be approximated as

[TABLE]

where we have recall that $h=\frac{1}{N}$ grid parameter for the grid defined in figure (19). The relative error is defined as

[TABLE]

The discretization error of our scheme is of order $h^{2}$. It is natural to ask how many iterations we need in order for the relative error to be of this size. There is no point in iterating more than this.

[TABLE]

[TABLE]

[TABLE]

[TABLE]

[TABLE]

Gauss elimination requires in general $O(N^{3})$ and since for large $N$ we have $N^{2}\ln N\ll N^{3}$ we observe that Gauss-Jacobi iteration is much faster that Gauss elimination.

5.4.3 Gauss-Seidel

We use the same decomposition of $A$ as for the Gauss-Jacobi method, $A=L+D+U$

[TABLE]

[TABLE]

[TABLE]

[TABLE]

[TABLE]

where $M_{S}=-(L+D)^{-1}U$, $\mathbf{c}=(L+D)^{-1}\mathbf{b}$. In component form we have the iteration

[TABLE]

We will now see that Gauss-Seidel converges faster than Gauss-Jacobi. It also use the memory in a more effective way. For the stability of this method we must find the spectrum of the matrix

[TABLE]

Note that the iteration matrix for the Gauss-Seidel method can be rewritten into the form

[TABLE]

Using straight forward matrix manipulations it is now easy to verify that $\lambda_{s}$ is an eigenvalue to the matrix $M_{S}$ iff

[TABLE]

easy to show that Let now $\alpha$ be any number, and define a matrix $E$ by

[TABLE]

Doing the matrix multiplication, it is easy to verify that we have the identity

[TABLE]

from which we immediately can conclude

[TABLE]

This derivation only require $D$ to be a diagonal matrix, $L$ to be a lower triangular matrix and $U$ to be a upper triangular matrix. Using (59) and (60), we have

[TABLE]

Thus, we have proved that the nonzero eigenvalues $\lambda_{s}$ of the Gauss-Seidel iteration matrix are determined by the eigenvalues $\mu$ of the Jacobi iteration matrix through the formula

[TABLE]

Since all eigenvalues of the Jacobi iteration matrix $M_{J}$ are inside the unit circle, this formula implies that the same is true for the eigenvalues of the Gauss-Seidel iteration matrix $M_{S}$. Thus Gauss-Seidel also converges for all initial conditions. As for the rate of convergence, using the results from the analysis of the convergence of the Gauss-Jacobi method, we find that

[TABLE]

and therefore

[TABLE]

For a fixed accuracy the Gauss-Seidel only needs half as many iterations as Gauss-Jacobi.

These two examples have given a taste of what one needs to do in order to define iteration methods for linear systems and study their convergence.

As previously mentioned there is an extensive set of iteration methods available for solving large linear systems of equation, and many has been implemented in numerical libraries available for most common programming languages. The choice of method depends to a large extent on the structure of your matrix. Here are two common choices:

[TABLE]

The conjugate gradient method has also been extended to more general matrices. A standard reference for iterative methods is “Iterative methods for sparse linear systems” by Yoosef Saad.

6 First order PDEs

So far we have only discussed numerical methods for solving PDEs. It is now time to introduce the so-called analytical methods. The aim here is to produce some kind of formula(s) for the solution of PDEs. What counts as a formula is not very precisely defined, but in general it involves infinite series and integrals that can not be solved by quadrature. Thus, if the ultimate aim is to get a numerical answer, and it usually is, numerical methods must be used to approximate the said infinite series and integrals. There is thus no clear-cut line separating analytical and numerical methods.

We started our exposition of the finite difference method by applying it to second order equations. This is because the simplest finite difference methods work best for such equations.

For our exposition of analytical methods we start with first order PDEs because the simplest such methods work best for these types of equations.

6.1 Higher order equations and first order systems

It turns out that higher order PDEs and first order systems of PDEs are not really separate things. Usually one can be converted into the other. Whether this is useful or not depends on the context.

Let us recall the analog situation from the theory of ODEs. The ODE for the Harmonic Oscillator is

[TABLE]

Define two new functions by

[TABLE]

Then

[TABLE]

Thus we get a first order system of ODEs.

[TABLE]

In the theory of ODEs this is a very useful device that can be used to transform any ODE of order $n$ into a first order system of ODEs in $\mathbb{R}^{n}$.

Note that the transition from one high order ODE, to a system of first order ODEs is not by any means unique. If we for the Harmonic Oscillator rather define

[TABLE]

we get the first order system

[TABLE]

Sometimes it is also possible to convert a first order system into a higher order scalar equation. This is however not possible in general but can always be done for linear first order systems of both ODEs and PDEs.

Let us consider the wave equation

[TABLE]

Observe that

[TABLE]

Inspired by this let us define a function

[TABLE]

Then the wave equation implies

[TABLE]

Thus the wave equation can be converted into the system

[TABLE]

Observe that the second equation is decoupled from the first. We will later in this class see that this can be used to solve the wave equation.

Writing a higher order PDE as a first order system is not always useful. Observe that

[TABLE]

Thus we might write the Laplace equation as the system

[TABLE]

However this is not a very useful system since we most of the time are looking for real solutions to the Laplace equation.

If we rather put

[TABLE]

Then using the Laplace equation and assuming equality of mixed partial derivatives we have

[TABLE]

Thus we get the first order system

[TABLE]

These are the famous Cauchy-Riemann equations. If you have studied complex function theory you know that (62) is the key to the whole subject; they determine the complex differentiability of the function

[TABLE]

(62) is a coupled system and does not on its own make it any simpler to find solutions to the Laplace equation. However, by using the link between complex function theory and the Laplace equation, great insight into it’s set of solutions can be gained.

Going from first order systems to higher order equations is not always possible but sometimes very useful.

A famous case is the Maxwell equations for electrodynamics, which in a region containing no sources of currents and charge are

[TABLE]

where $\mathbf{E}=(E_{1},E_{2},E_{3})$ and $\mathbf{B}=(B_{1},B_{2},B_{3})$ are the electric and magnetic fields. We can write this first order system for $\mathbf{E}$ and $\mathbf{B}$ as a second order system for $\mathbf{E}$ alone. We have, using (63), that

[TABLE]

[TABLE]

[TABLE]

[TABLE]

[TABLE]

and

[TABLE]

[TABLE]

[TABLE]

This is the 3D wave equation and it predicts, as we will see later in this course, that there are solutions describing electromagnetic disturbances traveling at a speed $v=\frac{1}{\sqrt{\varepsilon_{0}\mu_{0}}}$ in source free regions (for example in a vacuum). The parameters $\varepsilon_{0},\mu_{0}$ are constants of nature that was known to Maxwell, and when he calculated $v$, he found that the speed of these disturbances was very close to the known speed of light.

He thus conjectured that light is an electromagnetic disturbance. This was a monumental discovery whose ramifications we are still exploring today.

6.2 First order linear scalar PDEs

A general first order linear scalar PDE with two independent variables, $(x,t)$, is of the form

[TABLE]

where in general, $a,b,c$ and $d$ depends on $x$ and $t$. The solution method for (64) is based on a geometrical interpretation of the equation.

Observe that the two functions $a(x,t)$ and $b(x,t)$ can be thought of as the two components of a vector field $\mathbf{f}(x,t)$ in the plane.

[TABLE]

Let $\boldsymbol{\gamma}(s)=(x(s),t(s))$ be a parametrized curve in the plane whose velocity at each point $(x,t)$ is equal to $\mathbf{f}(x,t)$.

Thus we have

[TABLE]

Let $v(x,t)$ be a solution of (64) and define

[TABLE]

Then

[TABLE]

Thus $v(s)$ satisfy an ODE!

[TABLE]

From (65) and (66) we get a coupled system of ODEs

[TABLE]

The general theory of ODEs tells us that through each point $(x_{0},t_{0},v_{0})$, there goes exactly one solution curve of (67) if certain regularity conditions on $a,b,c$ and $d$ are satisfied.

The family of curves that solves (67) are called the characteristic curves for the PDE (64).

Observe that the system (67) is partly decoupled, the two first equations can be solved independently of the third. The plane curves solving

[TABLE]

are called characterictic base curves. They are also, somewhat imprecisely, called characteristic curves.

Most of the time we are not interested in the general solution of (64), but rather in special solutions that in addition to (64) has given values on some curve $C$ in the $(x,t)$ plane.

Any such problem is called an initial value problem for our equation (64).

This is a generalization of the notion of initial value problem for a PDE from what we have discussed earlier. This earlier notion consisted of giving the solution at $t=0$ (or $t=t_{0}$)

[TABLE]

where $\varphi(x)$ is some given function of $x$.

Let $(x(\tau),y(\tau))$ be a parametrization of $C$. Then our generalized notion of initial value problem for (64) consists of imposing the condition

[TABLE]

where $\varphi(\tau)$ is a given function. Our previous notion of initial value problem corresponds to the choice.

[TABLE]

Since this is a parametrization of the x-axis in the $(x,t)$ plane. Thus in summary, the (generalized) initial value problem for (64) consists of finding functions

[TABLE]

such that

[TABLE]

Where $(x(\tau),t(\tau))$ parametrize the initial curve and where $\varphi(\tau)$ are the given values of $v$ on the initial curve. The reader might at this point object and ask how we can possibly consider (68), solving (69), to be a solution to our original problem (64). Where is my $v(x,t)$ you might ask? This is after all what we started out to find.

In fact, the function $v(x,t)$, solving our equation (64), is implicitly contained in (68).

We can in principle use the two first equations in (68) to express $s$ and $\tau$ as functions of $x$ and $t$.

[TABLE]

It is known that this can be done locally close to a point $(s_{0},\tau_{0})$ if the determinant of the Jacobian matrix at $(s_{0},\tau_{0})$, this is what is known as the Jacobian in calculus, is nonzero.

[TABLE]

Thus if the Jacobian is nonzero for all points on the initial curve the initial value problem for (64) has a unique solution for any specification of $\varphi$ on the initial curve. This solution is

[TABLE]

Actually finding a formula for the solution $v(x,t)$ using this method, which is called the method of characteristics, can be hard. This is because in general it is hard to find explicit solutions to the coupled ODEs (69) and even if you succeed with this, inverting the system (70) analytically might easily be impossible.

However, even if we can not find an explicit solution, the method of characteristics can still give great insight into what the PDE is telling us. Additionally, the method of characteristic is the starting point for a numerical approach to solving PDEs of the type (64) and of related types as well.

In our traffic model, we encountered the equation

[TABLE]

This equation is of the type (64) and can be solved using the method of characteristic. Let the initial data be given on the x-axis.

[TABLE]

Our system of ODEs determining the base characteristic curves for the traffic equation (71) is thus

[TABLE]

Observe that the last ODE tell us that the function $n$ is constant along the base characteristic curves.

[TABLE]

Now we must express $s$ and $\tau$ in terms of $x$ and $t$

[TABLE]

and our solution to the initial value problem is

[TABLE]

This is the unique solution to our initial value problem. The base characteristic curves are given in parametrized form as

[TABLE]

or in unparametrized form as

[TABLE]

In either case they are a family of straight lines that is parametrized by $\tau$.

The solution in general depends on both the choice of initial curve and which values we specify for the solution on this curse. In fact, if we are not careful here, there might be many solutions to the initial value problem or none at all.

Since $n$ is constant along a base characteristic curve, its value is determined by the initial value on the point of the initial curve where the base characteristic curve and the initial curve coincide. If this occurs at more than one point there is in general no solution to the initial value problem.

Let for example $C$ be given through a parametrization

[TABLE]

The initial values is some function $\varphi(\tau)$ defined on $C$. The initial curve is thus a parabola and the base characteristic curves form a family of parallel lines

[TABLE]

Since $x-c_{0}t=A$ intersect the initial curve at two points corresponding to $x=\tau_{1}$ and $x=\tau_{2}$ and since $n=\varphi(\tau)$ is constant on the base characteristic curves, we must have

[TABLE]

This only holds for very special initial data. At a certain point $x=\tau_{3},\,t=\tau_{3}(1-\tau_{3})$, the characteristic base curve is tangent to the initial curve.

At this point the equations

[TABLE]

can not be inverted, even locally.

In order to see this, we observe that the base characteristic curves for the traffic equation (71), with initial data prescribed on the parabola (72), are given by

[TABLE]

The Jacobian is

[TABLE]

and at the point of tangency we have

[TABLE]

[TABLE]

[TABLE]

[TABLE]

[TABLE]

A point where the base characteristic curves and the initial curve are tangent is called a characteristic point. The initial curve is said to be characteristic at this point. Thus the initial curve $C$ is characteristic at the point $(x,t)$ where

[TABLE]

If the initial curve is characteristic at all of its points, we say that the initial curve is characteristic. It in fact is a base characteristic curve. The initial value problem then in general has no solution, and if it does have a solution, the solution is not unique. A characteristic initial value problem for our traffic equation (71) would be

[TABLE]

[TABLE]

General solution of the equation is

[TABLE]

The initial values require that

[TABLE]

[TABLE]

[TABLE]

In order for a solution to exist, $\varphi(t)=A$ must be a constant. And if $\varphi(t)=A$ then (74), with any function $f(y)$ where $f(0)=A$, would solve the initial value problem. Thus we have nonuniqueness.

Let us next consider the initial value problem

[TABLE]

The initial curce is the horizontal line $t=1$ which can be parametrized using

[TABLE]

Applying the method of characteristics, we get the following system of ODEs

[TABLE]

and

[TABLE]

We now solve for $s$ and $\tau$ in terms of $x$ and $t$

[TABLE]

The first equation only have solutions for $t>0$

[TABLE]

[TABLE]

[TABLE]

This solution is constant along the base characteristic curves

[TABLE]

[TABLE]

[TABLE]

The solution we have found exists only for $t>0$. This is its domain of definition

As a final example, let us consider the initial value problem for the wave equation

[TABLE]

We have previously seen that this can be written as a first order system

[TABLE]

[TABLE]

The initial conditions are

[TABLE]

We start by solving the homogeneous equation (78). We will use the method of characteristics on the problem

[TABLE]

[TABLE]

The corresponding system of ODEs is

[TABLE]

Solving we find

[TABLE]

We now express $\tau$ and $s$ as function of $x$ and $t$

[TABLE]

from which it follows that the function $u(x,t)=u(s(x,t),\tau(x,t))$ is given by

[TABLE]

Let us next apply the method of characteristic to the initial value problem for equation (77). This problem is

[TABLE]

[TABLE]

The corresponding system of ODEs is

[TABLE]

The first two equations give

[TABLE]

The last equation is solved by

[TABLE]

and using the initial data we have

[TABLE]

[TABLE]

[TABLE]

Observing that

[TABLE]

the formula for $v$ become

[TABLE]

Changing variable in the integral using $y=-2\gamma\sigma+\tau$, we get

[TABLE]

and $s=t,\,\tau=x+\gamma t$ from which it follows that

[TABLE]

Therefore, the function $v(x,t)=v(s(x,t),\tau(x,t))$ is given by the formula

[TABLE]

[TABLE]

[TABLE]

This is the famous D’Alembert formula for the solution of the initial value problem for the wave equation on an infinite domain. It was first derived by the French mathematician D’Alembert in 1747.

Let $G(\xi)=\int g(y)\mathrm{d}y$. Then the formula can be rewritten as

[TABLE]

The solution is clearly a linear superposition of a left moving and a right moving wave.

6.3 First order Quasilinear scalar PDEs

A first order, scalar, quasilinear PDE is an equation of the general form

[TABLE]

where

[TABLE]

This equation can also be solved through a geometrical interpretation of the equation.

Let $u(x,t)$ be a solution of (79).The graph of $u$ is a surface in $(x,t,u)$ space defined by the equation.

[TABLE]

This surface is called an integral surface for the equation.

From calculus we know that the gradient of $F$

[TABLE]

is normal to the integral surface $F$. Observe that the equation (79) can be rewritten as

[TABLE]

[TABLE]

[TABLE]

Thus the vector $(a,b,c)$ must always be in the tangent plane to the integral surface corresponding to a solution $u(x,t)$.

At each point $(x,t,u)$ the vector $(a,b,c)$ determines a direction that is the characteristic direction at this point. Thus the vector $(a,b,c)$ determines a characteristic direction field in $(x,t,u)$ space.

The characteristic curves are the curves that has $(a,b,c)$ as velocity at each of its points. Thus, if the curve is parametrized by $(x(s),t(s),u(s))$, we must have

[TABLE]

Observe that for this quasilinear case the two first equations does not decouple and the 3 ODEs has to be solved simultaneously. This can obviously be challenging.

For sufficiently smooth $a,b$ and $c$, there is a solution curve, or characteristic curve, passing through any point $(x_{0},t_{0},u_{0})$ in $(x,t,u)$ space.

The initial value problem is specified by giving the value of $u$ on a curve in $(x,t)$ space. This curve together with the initial data determines a curve, $C$, in $(x,t,u)$ space that we call the initial curve. To solve the initial value problem we pass a characteristic curve through every point of the initial curve. If these curves generate a smooth surface, the function, whose graph is equal to the surface, is the solution of the initial value problem.

Let the initial curve $C$ be parametrized by

[TABLE]

The characteristic system of ODEs with initial conditions imposed is

[TABLE]

The theory of ODEs gives a unique solution

[TABLE]

Observe that

[TABLE]

so the Jacobian is

[TABLE]

Let us assume that the Jacobian is nonzero along a smooth initial curve. We also assume that $a,b$ and $c$ are smooth. Then $\Delta(s,\tau)$ is continuous in $s$ and $\tau$, and therefore nonzero in an open set, surrounding the initial curve.

By the inverse function theorem, the equations

[TABLE]

can be inverted in an open set surrounding the initial curve. Thus we have

[TABLE]

Then

[TABLE]

is the unique solution to the initial value problem

[TABLE]

In general, going through this program and finding the solution $u(x,t)$ can be very hard. However deploying the method can give great insight into a given PDE.

Let us consider the following quasilinear equation

[TABLE]

This equation occurs as a first approximation to many physical phenomena. We have seen how it appears in the modelling of traffic.

It also describe waves on the surface of the ocean. If we assume that the waves moves in one direction only, then (81) appears as a description of the deviation of the ocean surface from a calm, flat state.

It is called the (inviscid) Burger’s equation. Let us solve Burger’s equation using the method of characteristics. Our initial value problem is

[TABLE]

The initial curve is parametrized by

[TABLE]

Here

[TABLE]

Let us compute the Jacobian along the initial curve. If the characteristic curves are

[TABLE]

where

[TABLE]

Then

[TABLE]

[TABLE]

[TABLE]

so $\Delta\neq 0$ along the entire initial curve. Let us next solve (82). From the last equation in (82) we get

[TABLE]

Thus $u$ is a constant independent of $s$. Since $u$ is a constant, the second equation in (82) give us

[TABLE]

and

[TABLE]

[TABLE]

[TABLE]

Finally for $t$ we get in the usual way

[TABLE]

If the conditions for the inverse function theorem holds, the equations

[TABLE]

can be inverted to yield

[TABLE]

and the solution to our problem is

[TABLE]

Without specifying $\varphi$, we can not invert the system and find a solution to our initial value problem. We can however write down an implicit solution. We have

[TABLE]

Even without precisely specifying $\varphi$ we can still gain great insight by considering the way the characteristics depends on the initial data. Projection of the characteristics in $(x,t,u)$ space down into $(x,t)$ space give us base characteristics of the form

[TABLE]

This is a family of straight lines in $(x,t)$ space parametrized by the parameter $\tau$. Observe that the slope of the lines depends on $\tau$, thus the lines in the family are not parallel in general. This is highly significant. The slope of the base characteristics is

[TABLE]

and is evidently varying as a function of $\tau$. This implies the possibility that the base characteristics will intersect. If they do, we are in trouble because the solution, $u$, is constant along the base characteristics. Thus if the base characteristics corresponding to parameters $\tau_{1}$ and $\tau_{2}$ intersect at a point $(x^{*},t^{*})$ in $(x,t)$ space then at this point $u$ will have two values

[TABLE]

that are different unless it happened to be the case that

[TABLE]

But in this case the two characteristics have the same sloop and would not have intersected. The conclusion is that base characteristics can intersect and when that happens $u$ is no longer a valid solution.

There is another way of looking at this phenomenon. Recall that the implicit solution is

[TABLE]

We know that a formula

[TABLE]

describe a wave moving at speed $c$ to the right. The formula (83) can then be interpreted to say that each point on the wave has a “local” speed $u(x,t)$ that vary along the wave. Thus the wave does not move as one unit but can deform as it moves.

If $u(x,t)>0$, the point $(x,u(x,t))$ on the wave move to the right and if $u(x,t)<0$, it moves to the left. If $u(x,t)$ is increasing, points on the wave move at greater and greater speed as we move to the right. The wave will stretch out. If $u(x,t)$ is decreasing the wave will compress. In a wave that has the shape of a localized “lump”,

there are points to the left that move faster than points further to the right. Thus points to the left can catch up with points to the right and lead to a steepening of the wave

Eventually the steepening can create a wave with a vertical edge. At this point $|u_{x}|=\infty$ and the model no longer describe what it is supposed to model. We thus have a breakdown.

It is evidently important to find out if and when such a breakdown occurs. There are two ways to work this out. The first approach use the implicit solution.

[TABLE]

We ask when the slope $u_{x}$ become infinite. From (84) we get

[TABLE]

[TABLE]

[TABLE]

Clearly $|u_{x}|$ become infinite if the denominator become zero

[TABLE]

[TABLE]

[TABLE]

where we recall that $\tau=x-ut$. What is important is to find the smallest time for which infinite sloop occurs. The solution is not valid for times larger than this. Thus we must determine

[TABLE]

When initial data is given, $\varphi(\tau)$ is known and $t^{*}$ can be found.

The second method involve looking at the geometry of the characteristics. It leads to the same breakdown time $t^{*}$.

We can immediately observe that if $\varphi$ is increasing for all $x$ so that

[TABLE]

then $t^{*}<0$ and no breakdown occur. For breakdown to occur we must have $\varphi^{\prime}(\tau)<0$ for some $\tau$. Let us look at some specific choices of initial conditions.

1. i)

$u(x,0)=\varphi(x)=-x$

[TABLE]

[TABLE]

[TABLE]

Thus we have breakdown at $t=1$. The solution is

[TABLE]

The whole graph becomes vertical at the breakdown time $t^{*}=1$. This is illustrated in figure 32. Let us describe this same case in terms of base characteristics. They are

[TABLE]

They all intersect the t-axis at $t=1$ and the base characteristics corresponding to $\tau$ intersect the x-axis at $x=\tau$. This case is illustrated in figure 33.

1. ii)

$u(x,0)=\varphi(x)=1-x^{2}$.

[TABLE]

[TABLE]

[TABLE]

Thus we have breakdown already at $t=0$! The implicit formula for the solution now is

[TABLE]

[TABLE]

[TABLE]

[TABLE]

[TABLE]

[TABLE]

[TABLE]

[TABLE]

[TABLE]

The initial condition requires that

[TABLE]

and for small $t$ we have

[TABLE]

Using the ‘ - ’ solution, we have for small $t$

[TABLE]

Thus

[TABLE]

Our unique solution of the initial value problem is therefore

[TABLE]

We have

[TABLE]

and thus infinite sloop occurs at all space-time points $(x,t)$ determined by

[TABLE]

[TABLE]

[TABLE]

This curve is illustarted in figure (34). There clearly are space-time points with arbitrarily small times on this curve and we therefore conclude that the breakdown time is $t^{*}=0$.

The actual wave propagation is illustrated in figure 35.

The steepening is evident.

1. iii)

$u(x,0)=\varphi(x)=\sin x.$

The implicit solution is now

[TABLE]

and thus

[TABLE]

From this formula it is evident that

[TABLE]

and this minimum breakdown time is achieved for each $\tau_{n}=(2n+1)\pi$.

Recall that the characteristic curves are parametrized by $\tau$ through

[TABLE]

Thus we can conclude that the breakdown occurs at the space-time points $(t^{*},x_{n}^{*})$ where

[TABLE]

The breakdown is illustrated in figure 36.

7 Classification of equations: Characteristics

We will start our investigations by considering second order equations. More specifically we will consider scalar linear second order PDEs with two independent variables. The most general such equation is of the form

[TABLE]

where $u=u(x,y),A=A(x,y)$ etc. The wave equation, diffusion equation and Laplace equation are clearly included in this class.

What separates one equation in this class from another are determined entirely by the functions $A\rightarrow G$. Thus the exact form of these functions determines all properties of the equation and the exact form of the set of solutions (or space of solutions as we usually say). In fact if $G=0$ the space of solutions is a vectorspace. This is because the differential operator

[TABLE]

is a linear operator. (Prove it!) Thus even if we at this point knows nothing about the functions forming the space of solutions of the given equation, we do know that it has the abstract property of being a vectorspace.

Studying abstract properties of the space of solutions of PDEs is a central theme in the modern theory of PDEs. Very rarely will we be able to get an analytical description of the solution space for a PDE.

Knowing that a given PDE has a solution space that has a set of abstract properties is often of paramount importance. This knowledge can both inform us about expected qualitative properties of solutions and act as a guide for selecting which numerical or analytical approximation methods are appropriate for the given equation.

The aim of this section is to describe three abstract properties that the equation (85) might have or not. These are the properties of being hyperbolic, parabolic or elliptic.

Let us then begin: We first assume that $A(x,y)$ is nonzero in the region of interest. This is not a very restrictive assumption and can be achieved for most reasonable $A$. We can then divide (85) by $A$

[TABLE]

We will now focus on the terms in (86) that are of second order in the derivatives. These terms form the principal part of the equation. The differential operator defining the principal part is

[TABLE]

We now try to factorize the differential operator by writing it in the form

[TABLE]

where $a,b,c$ are some unknown functions. We find these functions by expanding out the derivatives and comparing with the expression (87)

[TABLE]

Comparing with (87) we get

[TABLE]

From the first two equations we find that both $a$ and $b$ solves the equation

[TABLE]

Let us call these solutions $\omega^{\pm}$ so that

[TABLE]

We thus have

[TABLE]

and then

[TABLE]

We are in general interested in real solutions to our equations and aim to use this factorization to rewrite and simplify our equation. Thus it is of importance to know how many real solutions equation (89) has. From formula (90) we see that this is answered by studying the sign of $B^{2}(x,y)-A(x,y)C(x,y)$ in the domain of interest. We will here assume that this quantity is of one sign in the region of interest. This can always be achieved by restricting the domain, if necessary.

The equation (85) is called

1. i)

Hyperbolic if $B^{2}-AC>0,$ 2. ii)

Parabolic     if $B^{2}-AC=0,$ 3. iii)

Elliptic      if $B^{2}-AC<0.$

For the wave equation

[TABLE]

we have $A=1,B=0,C=-c^{2}$ ($x\to t,y\to x$ in the general form of the equation (85))

[TABLE]

Thus the wave equation is hyperbolic.

For the diffusion equation we have

[TABLE]

we have $A=0,B=0,C=-D$

[TABLE]

Thus the diffusion equation is parabolic

For the Laplace equation we have

[TABLE]

we have $A=1,B=0,C=1$

[TABLE]

Thus the Laplace equation is elliptic.

We will now start deriving consequences of the classification into hyperbolic, parabolic and elliptic.

7.1 Canonical forms for equations of hyperbolic type

For this type of equations $\omega^{+}\neq\omega^{-}$ and both $\omega^{+}$ and $\omega^{-}$ are real. Take $\omega^{+}$ and write down the following ODE

[TABLE]

Let the solution curves of this equation be parametrized by the parameter $\xi$. Thus the curves can be written as

[TABLE]

In a similar way we can write the solution curves for the ODE

[TABLE]

in the form

[TABLE]

where $\eta$ is the parameter parameterizing the solution curves for (92). We now introduce $(\xi,\eta)$ as new coordinates in the plane. The change of coordinates is thus

[TABLE]

For this to really be a coordinate system the Jacobian has to be nonzero.

[TABLE]

The solution curves, or characteristic curve, of the equation (85) are by definition determined by

[TABLE]

Differentiate implicit with respect to $x$. This gives, upon using (91) and (92),

[TABLE]

Thus

[TABLE]

because, by assumption, $\omega^{+}\neq\omega^{-}$ in our domain. Thus the proposed change of coordinates (93) does in fact define a change of coordinates.

Recall that $z=\omega^{+}$ is a solutions to the quadratic equation (89), thus

[TABLE]

Inserting $\omega^{+}=\frac{\xi_{x}}{\xi_{y}}$ from equation (94) we get, after multiplication by $\xi_{y}^{2}$

[TABLE]

In a similar way the root $z=\omega^{-}$ show that the function $\eta(x,y)$ satisfy the equation

[TABLE]

Thus both $\xi(x,y)$ and $\eta(x,y)$, that determines the characteristic coordinate system are solutions to the characteristic PDE corresponding to our original PDE (85)

[TABLE]

Using the chain rule we have

[TABLE]

[TABLE]

because

[TABLE]

Since

[TABLE]

describe the solution curves to the equation (91). In a similar way we find

[TABLE]

Therefore, using equations (94), we have

[TABLE]

and both operators are nonzero because $\omega^{+}\neq\omega^{-}$. But then we have

[TABLE]

Substituting this into (85) and expressing all first order derivatives with respect to $x$ and $y$ in terms of $\partial_{\xi}$ and $\partial_{\eta}$, we get an equation of the form

[TABLE]

where we have divided by the nonzero quantity $\xi_{y}\eta_{y}(w^{+}-w^{-})$ and where now $u=u(\xi,\eta)$, $a=a(\xi,\eta)$ etc. For any particular equation $a,b,c$ and $d$ will be determined explicitely in terms of $A,B,C$ and $D$.

Since (100) results from (85) by a change of coordinates, (100) and are equivalent PDEs. Their spaces of solutions will consist of functions that are related through the change of coordinates (93).

(100) is called the canonical form for equations of hyperbolic type.

There is a second common form for hyperbolic equation. We find this by doing the additional transformation

[TABLE]

[TABLE]

[TABLE]

for some new functions $\widetilde{a}=\widetilde{a}(\alpha,\beta)$ etc. This is also called the canonical form for hyperbolic equations. We observe that the wave equation

[TABLE]

is exactly of this form.

7.2 Canonical forms for equations of parabolic type

For parabolic equations we have

[TABLE]

in the domain of interest. For this case we have

[TABLE]

and there exists only one family of characteristic curves

[TABLE]

where the curves

[TABLE]

by definition are solutions to the ODE

[TABLE]

Let

[TABLE]

be another family of curves chosen such that

[TABLE]

defines a new coordinate system $(\xi,\eta)$ in the domain of interest. For example if $A,B$ and $C$ are constants we have

[TABLE]

and we can let

[TABLE]

where $c\neq\omega^{2}$. Then (102) is a coordinate system because

[TABLE]

We now express all derivatives $u_{x},u_{xx},u_{xy},u_{y},u_{yy}$ in terms of $u_{\xi},u_{\eta},u_{\xi\eta},u_{\xi\xi},u_{\eta\eta}$. using the chain rule. For example we have

[TABLE]

etc.

Inserting these expressions into the equation we get a principal part of the form

[TABLE]

Since

[TABLE]

we have

[TABLE]

Thus, the term $u_{\xi\xi}$ drops out of the equation. Also recall that for the parabolic case we have $\omega=-\frac{B}{A}$. Therefore we have

[TABLE]

Therefore the term $u_{\xi\eta}$ also drops out of the equation. Also, observe that the term

[TABLE]

must be nonzero. This is because if it was zero, $\eta=\eta(x,y)$ would define a second family of characteristic curves different from $\xi=\xi(x,y)$ and we know that no such family exists for the parabolic case.

We can therefore divide the equation by the term (104) and we get

[TABLE]

This is the canonical form for equations of parabolic type. We observe that the diffusion equation

[TABLE]

is of this form.

7.3 Canonical forms for equations of elliptic type

For elliptic equations we have

[TABLE]

so both $w^{+}$ and $w^{-}$ are complex and $\xi$ and $\eta$ become complex coordinates. In [1] all functions are extended into the complex domain, and an argument like the one for the hyperbolic case lead to the canonical form

[TABLE]

where $a=a(\alpha,\beta)$ etc. The same canonical form can be found using only real variables. The key step in the reduction to canonical form is to simplify the principal part of the equation using a change of variables. Let us therefore do a general change of variables

[TABLE]

From equation (7.2) we have the following expression for the principal part of the transformed equation

[TABLE]

The functions $\xi(x,y)$ and $\eta(x,y)$ are at this point only constrained by the Jacobi condition

[TABLE]

In the elliptic case, both solutions $\xi,\eta$ to the characteristic equation are complex. Thus, for the elliptic case the coefficients in front of $u_{\xi\xi}$ and $u_{\eta\eta}$ can not be made to vanish. However if we let $\eta(x,y)$ be arbitrary and let $\xi=\xi(x,y)$ be the solution family of the system

[TABLE]

Then along the solution curves $\xi(x,y)=\text{const},$ we have

[TABLE]

[TABLE]

[TABLE]

Thus

[TABLE]

because the characteristic equation has no real solutions.

Therefore

[TABLE]

is a change of coordinates.

Furthermore

[TABLE]

Thus, the coefficient in front of $u_{\xi\eta}$ is zero. Since the quantity

[TABLE]

the coefficients in front of $u_{\xi\xi}$ and $u_{\eta\eta}$ are nonzero and of the same sign. Dividing, we find that the principal part takes the form

[TABLE]

where the function $h(x,y)$ is a positive function.

Let us at this point simplify our considerations by assuming that the coefficients of the equation (85) are in fact constant. If that is the case, we choose the function $\eta(x,y)$ to be linear in $x$ and $y$. Then the derivatives of $\eta$ are constant and therefore the function $\xi(x,y)$ defined in (107) is also a linear function. For this particular case, the function $h(\xi,\eta)$ is a positive constant which we denote by $h_{0}$. Thus equation (108) takes the form

[TABLE]

Introduce a change of variables defined by

[TABLE]

Using the chain rule we find that the principal part of an elliptic equation with constant coefficients, takes the form

[TABLE]

Thus the canonical form of elliptic equations with constant coefficients is

[TABLE]

One can show that also elliptic equations with variable coefficients can be put into canonical form (112) by a suitable change of coordinates. This argument however either relies on complex analysis, or some rather intricate real variable manipulations.

7.4 Characteristic curves

One consequence of the classification of equations into hyperbolic, parabolic and elliptic, is the possibility to reduce each class into a canonical form where the principal part has the form

[TABLE]

In these reductions, certain special curve families played an important role. We will now see that these special curve families reappear when we consider the possibility of solutions that contain singularities. These are solutions whose existence,or not, is a consequence of which class our equation belong to. This will give us a new window into what the abstract classification of equations into the three classes hyperbolic, parabolic and elliptic actually means.

Let our domain, $D$, be separated into two parts $D_{1}$ and $D_{2}$ by a curve.

Let $u_{1}$ and $u_{2}$ be solutions of the equation

[TABLE]

in $D_{1}$ and $D_{2}$. Define a function $U$ on the whole domain $D$ by

[TABLE]

We will assume that $u,u_{x},u_{y}$ are continuous on $D$ but that $u_{xx},u_{xy},u_{yy}$ can be discontinuous across the curve separating $D_{1}$ and $D_{2}$. Thus $u$ is not a solution in the usual sense to (113) on the whole domain $D$.

The curve is defined by

[TABLE]

Introduce a curve family

[TABLE]

Our curve $y=h(x)$ corresponds to $\xi=0$. Introduce a second curve family determined $\eta=\eta(x,y)$ in such a way that

[TABLE]

define a change of coordinates. For this to be the case the Jacobian must be nonzero everywhere.

We now transform our equation into the new coordinate system in the same way as on page 79, and get

[TABLE]

where we only has displayed the principal part of the equation. What we precisely mean when we say that $u_{xx},u_{xy}$ and $u_{yy}$ can be discontinous across $y=h(x)$ is that

[TABLE]

are continuous but that $u_{\xi\xi}$ can be discontinuous. Evaluating (114) above and below the curve $\xi=0\Leftrightarrow y=h(x)$, taking the difference and going to the limit as the points approach $y=h(x)$ from above and below we get

[TABLE]

where $\left.\left[u_{\xi\xi}\right]\right|_{\xi=0}$ is the jump in $u_{\xi\xi}$ across $y=h(x)$. Since $\left.\left[u_{\xi\xi}\right]\right|_{\xi=0}\neq 0$, we must have

[TABLE]

Thus $\xi(x,y)$ must satisfy the characteristic PDE and consequently

[TABLE]

must be a characteristic curve for the PDE (113).

Thus, discontinuities in the second derivatives can only occur on characteristic curves.

Let us see how this works out for the wave equation, diffusion equation and Laplace equation. We have seen that these three equations are respectively of hyperbolic type, parabolic type and elliptic type.

The wave equation is

[TABLE]

[TABLE]

The characteristic PDE is

[TABLE]

and the characteristic curves are

1. i)

$\varphi_{t}-\gamma\varphi_{x}=0\Leftarrow\varphi(x,t)=x+\gamma t,$ 2. ii)

$\varphi_{t}+\gamma\varphi_{x}=0\Leftarrow\varphi(x,t)=x-\gamma t.$

The two families of characteristic curves are

[TABLE]

Singularities, that is discontinuities in $u_{xx},u_{xy},u_{yy}$, can only occur along these curves. Singularities given in initial data at $t=0$ is propagated along the characteristic curves and persists forever. This behavior is typical for all hyperbolic equations.

The diffusion equation is

[TABLE]

[TABLE]

So the characteristic PDE is

[TABLE]

and characteristic curves $\varphi(x,t)=\text{const}$ must be of the form

[TABLE]

Singularities in the initial data given on the x-axis can not leave the x-axis and will not be seen at any point $(x,t)$ for $t>0$. Thus, all singularities in initial data are washed out instantly for $t>0$. This behaviour is typical for all parabolic equations.

Finally, the Laplace equation is

[TABLE]

and the characteristic equation is

[TABLE]

There are no characteristic curves that can support singularities. The solutions of elliptic equations are very smooth. This is typical for all elliptic equations.

7.5 Characteristic curves again: A different point of view

We first discussed characteristic curves in the context of linear and quasilinear first order scalar PDEs in two independent variables. There they played the star role in a method for solving the initial value problem for such equations. This was the method of characteristics. Next, we met them in the context of second order linear scalar PDEs in two independent variables. Here we have seen that they played a dual role. On the one hand they were instrumental in the reduction of hyperbolic and parabolic equations to canonical form, and on the other hand they where curves along which solutions to the PDEs could have discontinuities in their second derivatives. This second role is also played by the characteristic of first order linear scalar PDEs.

We will now discuss a third way in which characteristic curves play a role.

Let us start by considering the case of first order linear scalar PDEs in two independent variables. The initial value problem in assumed to be given in terms of an initial curve that is the graph of a function $h(x)$

[TABLE]

On this curve the solution $u$ is assumed to be given by some function $f(x)$. Thus we have

[TABLE]

[TABLE]

Differentiating (116) we get

[TABLE]

and evaluating (115) along the initial curve we get

[TABLE]

(117) and (118) is a linear system of equations for determining $u_{x}$ and $u_{y}$

[TABLE]

For points close to the initial curve we can express the solution of the initial value problem (115),(116) as a Taylor series.

[TABLE]

Clearly if the solution $u(x,y)$ is to exist and be unique, then at least $u_{x}(x_{0},y_{0})$ and $u_{y}(x_{0},y_{0})$ must exist as a unique solution of the system (119).

This is ensured if

[TABLE]

Even if this is true it is not certain that the initial value problem is uniquely solvable, there is the question of the coefficients of the rest of the Taylor series and of course there is the matter of convergence for such series ….

However, if on the other hand,

[TABLE]

$u_{x}(x_{0},y_{0}),u_{y}(x_{0},y_{0})$ either does not exist and if they do exist they are not unique. This follows from the elementary theory of linear systems of equations. So if (120) holds we can know for sure that we do not have existence and uniqueness for the initial value problem (115),(116).

Condition (120) can be written as

[TABLE]

[TABLE]

[TABLE]

This equation tells us that $y=h(x)$ must be a characteristic curve! This follow, because the family of base characteristic curves $(x(s),y(s))$ are determined by being solutions to

[TABLE]

If the curves in question are graphs of functions we can parametrise them using $x$. So that they are given by a function $y(x)$. For this function (122) together with the chain rule gives

[TABLE]

and (121) say that $h(x)$ solves this equation so that $y=h(x)$ is a base characteristic curve.

Let us next consider the case of second order linear scalar PDEs in two independent variables.

The characteristic curves have preciously been found to satisfy

[TABLE]

The curves on which the PDE could have singularities, $y=h(x)$, were related to $\varphi(x,y)$ through

[TABLE]

and (123) then give the equation

[TABLE]

The solutions of this nonlinear ODE gives the characteristic curves for PDEs of the form

[TABLE]

where

[TABLE]

Recall that all parameters are functions of $x$ and $y$.

Let $y=h(x)$ be a initial curve where we have given $u,u_{x}$ and $u_{y}$. Thus

[TABLE]

For points close to the initial curve the solution can be expand as a Taylor series.

[TABLE]

Differentiating the last two equations in the initial data (7.5) with respect to $x$, we get

[TABLE]

[TABLE]

Evaluating (125) along $y=h(x)$ we get the equation

[TABLE]

Equations (127),(127) and (129) define a linear system of equations for $u_{xx},u_{xy}$ and $u_{yy}$

[TABLE]

The Taylor expansion show that the the problem (125) with initial data (7.5) can have a unique solution only if the determinant of the matrix in the linear system (130) is nonzero.

If, on the other hand, the determinant is zero, the problem is sure to either have no solution or to have a nonunique solution. This condition is

[TABLE]

and this is exactly the equation for the characteristic curves (124).

Thus characteristic curves are curves for which the (generalized) initial value problem either can not be solved or for which the solution is non-unique.

There is something awkward about the presentation given and equation (132). Since we are considering the general initial value problem as giving $u,u_{x}$ and $u_{y}$ along some curve in the $(x,y)$ plane, there is really no point in restricting to curves that is the graph of a function $y=h(s)$. Why not consider curves that are graph of a function $x=k(y)$ or indeed a curve that is neither a graph of a function of $x$ or a function of $y$?

Let us consider case iii). Let $\varphi(x,y)$ be a function such that $C$ is a level curve of $\varphi$

[TABLE]

Pick any point $(x_{0},y_{0})$ on $C$. Then the piece of the curve close to $(x_{0},y_{0})$ will either be the graph of a function of $x$ or the graph of a function of $y$

Let us assume, without loss of generality, that the piece of $C$ close to $(x_{0},y_{0})$ is the graph of a function of x, $y=h(x)$. But then we must have

[TABLE]

[TABLE]

[TABLE]

[TABLE]

[TABLE]

Inserting this into equation (132) gives

[TABLE]

[TABLE]

[TABLE]

The conclusion of this is that curves, $C$, for which the generalized initial value problem $(\ref{eq7-34}),(\ref{eq7-37})$ either has no solution or has a nonunique solution, must be level curves

[TABLE]

where $\varphi$ is a solution of (133). And all solutions of (133) determines such special initial curves.

We have met equation (133) before on page 9 of these notes where it was called the characteristic equation for $(\ref{eq7-34})$.

Thus the special initial curves discussed in this section are in fact characteristic curves as defined previously.

The characteristic curves are seen to play three different roles for PDEs of type (125).

1. i)

Level curves for special coordinate system used to reduce the PDE to canonical form. 2. ii)

Curves on which solutions of the PDE can have singularities. 3. iii)

Curves for which the initial value problem for the PDE either have no solution or have a non-unique solution.

So far the statements i), ii) and iii) only applies to linear second order scalar PDEs in two independent variables. Can we extend these results to more general equations?

The answer is yes we can! Many of these ideas can be extended far beyond the class of equations we have discussed so far. In [1] classification and associated characteristic curves and surfaces for large classes of equations are discussed on pages 137-151.

In this course we will not discuss all these various cases but will rather restrict to linear scalar second order equations in $n$ independent variables where $n$ can be any positive integer.

Equations of this type takes the form

[TABLE]

where all parameters $\{a_{ij}\}$ and $\{b_{i}\},c,d$ are functions of $(x_{1},\dots,x_{n})$

Remark on notation: Here and in many other areas of applied mathematics and theoretical physics it is common to use the so-called Einstein summation convention.

In this convention repeated indices in any expression are implicitly summed over. Thus the expressions

[TABLE]

[TABLE]

[TABLE]

means

[TABLE]

[TABLE]

[TABLE]

Using this convention, equation (134) is written as

[TABLE]

Observe that

[TABLE]

where $\widetilde{a}_{ij}=\frac{1}{2}(a_{ij}+a_{ji})$.

$\widetilde{a}_{ij}$ are evidently the components of a symmetric matrix

[TABLE]

Thus, we can without loss of generality assume that the matrix $A=(a_{ij})$ in equation (135) is a symmetrix matrix.

Let us first assume that all $a_{ij}$ are constants. Let $\lambda_{k}$ be the eigenvalues of $A$ and $R=(r_{ki})$ be the orthogonal matrix whose rows are the correponding eigenvectors. In terms of matrices we have

[TABLE]

This can be written, using the Einstein convention, in the form

[TABLE]

where $\Lambda_{pq}$ is defined by

[TABLE]

Introduce new coordinates $\xi_{1},\dots,\xi_{n}$ through

[TABLE]

This is evidently a linear change of coordinates. Observe that

[TABLE]

[TABLE]

[TABLE]

Thus the second order mixed partial derivatives all vanish. This is the canonical form for (135). If we return to the $n=2$ case we observe that

[TABLE]

Based on this we now make the following classification of equations of type (135).

[TABLE]

For the elliptic case we can make the further change of coordinates

[TABLE]

leading to a principal part of the form

[TABLE]

Assuming $\lambda_{1}>0$ and $\lambda_{2}\dots\lambda_{n}<0$, the hyperbolic case is reduced, using (136), to the form

[TABLE]

Similar reductions are possible for the parabolic case. Some $n=3,4$ basic examples of elliptic, hyperbolic and parabolic equations are

[TABLE]

The classification of equations of type (135) into elliptic, hyperbolic and parabolic based on eigenvalues of the matrix $A$ is still possible when $a_{ij}$ depends on $(x_{1},\dots,x_{n})$ and are thus not constant. The eigenvalues and the classification will however in general depend on which point we are at. Also, it is in general not possible to reduce it to a canonical form where all terms with mixed partial derivatives vanish. This is seen to be true by a simple counting argument.

Let

[TABLE]

be a change of coordinates in $\mathrm{R}^{n}$. Using the chain rule the principal part of equation (135) can be written as

[TABLE]

where $F_{ij}$ is some expression containing derivatives of the functions $\xi_{j}$. Let us try to eliminate all mixed partials. Then for this to happened, we must have

[TABLE]

There are exactly $\frac{1}{2}n(n-1)$ terms with mixed partial derivatives. Thus there are $\frac{1}{2}n(n-1)$ equations in (141). On the other hand we only have $n$ unknown functions $\xi_{1},\dots,\xi_{n}$. If

[TABLE]

the system (141) in general has no solution. This happened for $n\eqslantgtr 4$.

For equations of type (135), surfaces supporting singularities of solutions, and on which given initial data does not give a unique solution for the initial value problem, are given as solutions to

[TABLE]

These surfaces are called characteristic surfaces and equation (142) is called the characteristic equation for equation (135). For example the 3D wave equation

[TABLE]

has a characteristic equation

[TABLE]

One important type of solution is

[TABLE]

where $\omega^{2}=k_{1}^{2}+k_{2}^{2}+k_{3}^{2}$. For each value of $t$ the characteristic surface

[TABLE]

[TABLE]

[TABLE]

is a plane normal to the vector $\boldsymbol{k}=(k_{1},k_{2},k_{3})$. Such surfaces play a pivotal role in the theory of wave equations where they appear through special solutions called plane waves.

7.6 Stability theory, energy conservation and dispersion

So far we have used the highest order terms, the so-called principal part of the equation, to classify our equations. The lower order terms however also play an important role.

In this section of the notes we will investigate the effects of lower order terms by considering the class of equations of the form

[TABLE]

where $A,B,C,D,E$ and $F$ are real constants. Certain aspects of our conclusions for this class can be extended to a wider class of equations[1].

We start our investigations by looking for special solutions to (143) of exponential form

[TABLE]

This is a complex solution but both the real and imaginary part of $u$ are real solutions of (143). (Verify this). We will in this section interpret the variable $t$ as time and study the evolution of the solutions (144) for $t>0$. The special solutions (144) are called normal modes for (143). By the principle of linear super-position we can add together several normal modes corresponding to different values of $k$ to create more general solutions. Taking an infinite set of discrete values $\{k_{j}\}$ we get a solution of the form

[TABLE]

This is a Fourier series, and the sense in which such series solves (143) will be discussed extensively later in the class. There we will also discuss a much more general notion of normal mode that applies to a much wider class of equations. We will see that the normal modes (144) can be combined in a “continuous sum” by which we mean an integral

[TABLE]

to generate even more general solutions to (143). We recognize (146) as a Fourier transform. More general transforms are associated with other types of normal modes. These are discussed in chapter 5 in [1].

At this point, the only thing to note is that through the general representations (145) and (146) the behavior of the normal modes (144) will determine what any solution to equation (143) will do for $t>0$.

So, let us now start analyzing the normal modes. From (144) we have

[TABLE]

Inserting (7.6) into (143) we get

[TABLE]

A solution of this type thus only exists if

[TABLE]

This equation has two complex solutions. For each of these we can write

[TABLE]

and the corresponding normal mode can be written as

[TABLE]

The magnitude of $u(x,t)$ is

[TABLE]

We will assume that the expression $\left|a(k)\right|$ is bounded for all $k$. The growth of $|u(x,t)|$ as a function of $k$ for $t>0$ is then determined by $Re[\lambda(k)]$.

There are now two distinct possibilities. Either $Re[\lambda(k)]$ is bounded above for all $k\in\mathbb{R}$, or it is not.

Define

[TABLE]

If $Re[\lambda(k)]$ is unbounded we let $\Omega=+\infty$. Let us first consider the case $\Omega=+\infty$. Previously in this class we introduced the notion of being well-posed for initial and or boundary value problems for PDEs. In order to be well-posed, three criteria had to be satisfied

The problem must have at least one solution. 2. 2)

The problem must have at most one solution. 3. 3)

The solution must depend in a continuous way on the data of the problem: Small data must produce a small solution.

Data here means parameters occurring in the PDEs, description of the boundary curves or surfaces, if any, and functions determining boundary and/or initial conditions

For the case $\Omega=+\infty$, let $a(k)=\frac{1}{\lambda(k)^{2}}$ and consider the initial value problem for the normal mode with data at $t=0$. The data are from (144)

[TABLE]

Since $Re[\lambda(k)]$ are unbounded $(\Omega=+\infty)$ there are values of $k$ such that $|u(x,0)|$ and $|u_{t}(x,0)|$ are as small as we like. Thus the data in this initial value problem can be made arbitrarily small. However

[TABLE]

grows exponentially as a function of $k$ and can for any fixed $t>0$ be made arbitrarily large. This is because an algebraic decay $\frac{1}{\lambda^{2}(k)}$ is always overwhelmed by an exponential growth $e^{Re[\lambda(k)]t}$ (Basic calculus). Thus arbitrarily small data can be made arbitrarily large for any $t>0$. Condition 3) on page 52 is not satisfied and consequently the problem is not well posed. If the case $\Omega=+\infty$ is realized, it is back to the drawing board, something is seriously wrong with the model! If $\Omega<+\infty$ one can show that the problem is well posed.

Let us look at some specific examples.

The equation

[TABLE]

where $\rho$ is constant, is an elliptic equation. Our equation for $\lambda(k)$ for this case is

[TABLE]

[TABLE]

[TABLE]

$\lambda(k)$ is clearly unbounded when $k$ ranges over $\mathbb{R}$. Thus $\Omega=+\infty$ and the Cauchy problem for the equation with data at $t=0$ is not well posed.

The equation

[TABLE]

where $\rho$ is constant is a parabolic equation. In fact, we can also write it as

[TABLE]

For $\rho<0$ it is the heat equation. For $\rho>0$ it is called the reverse heat equation. We get the reverse heat equation from the heat equation by making the transformation $t^{\prime}=-t$. Increasing values of $t^{\prime}$ then corresponds to decreasing values of $t$. Solving the reverse heat equation thus corresponds to solving the heat equation backward in time.

We now get the following equation for $\lambda(k)$

[TABLE]

Clearly $\lambda(k)$ is unbounded for $\rho>0\Rightarrow\Omega=+\infty$ which implies that the reverse heat equation is not well posed. This say something profound about the process of heat diffusion or any other process described by the diffusion equation. Such processes are irreversible, they have a built-in direction of time. It is a well known fact that no fundamental physical theory is irreversible, they have no direction of time. So how come the diffusion equation has one? Where and how did the direction of time enter into the logical chain connecting microscopic reversible physics to macroscopic irreversible diffusion processes? This question has been discussed for more than a hundred years and no consensus has been reached. Nobody really knows.

For $\rho<0$, $\lambda(k)$ is bounded above by zero $\Rightarrow\Omega=0$ and the Cauchy problem for the heat equation is well-posed.

The equation

[TABLE]

where $\rho$ is a constant, is a hyperbolic equation. The equation for $\lambda(k)$ is

[TABLE]

[TABLE]

[TABLE]

For $\rho<0$ the real part of $\lambda$ vanish for $|k|$ large enough and for $\rho>0$ the real part vanish for all $k$. Thus in either case $\Omega<+\infty$ and the Cauchy problem for this equation is well-posed.

The constant $\Omega$ is called the stability index for the equation. If $\Omega>0$ there is some set or range of $k$ values whose corresponding normal modes experience exponential growth. If this is the case the equation is said to be unstable. If $\Omega<+\infty$ the Cauchy problem is well posed, but the unbounded growth can still invalidate the PDE as a mathematical model for some $t>0$. The reason for this is that PDE models are usually derived using truncated Taylor expansions or expansions of some other sort. In order for these truncations to remain valid the solutions must remain small. This assumption will break down in the case of unbounded growth. If $\Omega<0$ all normal modes decay exponentially. If this is the case the PDE is strictly stable. If $\Omega=0$ there may be some modes that does not decay and will continue to play a role for all $t>0$. If this is the case the PDE is neutrally stable.

Consider the parabolic equation

[TABLE]

For this equation the equation for $\lambda(k)$ is

[TABLE]

[TABLE]

[TABLE]

[TABLE]

Clearly $\Omega=-b$. It is stable if $b>0$ and unstable if $b<0$ and neutrally stable for $b=0$.

The equation

[TABLE]

is parabolic and has an equation for $\lambda$ given by

[TABLE]

[TABLE]

[TABLE]

[TABLE]

[TABLE]

and thus $\Omega=0$. The equation is therefore neutrally stable. It has a normal mode solution

[TABLE]

It is straight forward to verify that it also have a solution

[TABLE]

By linear superposition it has a solution

[TABLE]

The solution to the Cauchy problem with start data

[TABLE]

is then easily shown to be

[TABLE]

Thus, even if the data for this case is bounded, the solution grows without bound as $t$ approach infinity. However, the growth is algebraic rather than exponential, which was the case for an unstable equation. This kind of behavior is generally true for the neutrally stable case; if there is an instability, then the solution grows algebraically rather than exponentially.

A special situation occur when $Re[\lambda(k)]=0,\,\,\forall k$. Then

[TABLE]

The amplitude of the normal mode is preserved for all time. In many cases of importance, $|u|^{2}$ is a measure of the energy of the normal mode and we have that the energy in the system is conserved in time.

In general, an equation where $Re[\lambda(k)]=0$, is called conservative. For the conservative case we define the positive quantity

[TABLE]

to be the frequency of the mode. Given this definition we have

[TABLE]

depending on the sign of the imaginary part of the complex root $\lambda(k)$ of equation (152).

The normal mode then takes the form

[TABLE]

The equation

[TABLE]

is called the dispersion relation for the equation. If

[TABLE]

and thus $\omega(k)$ is not a linear function of $k$, the equation is said to be of dispersive type. The normal mode can, for the conservative case, be written as

[TABLE]

where $\theta=kx\pm\omega(k)t$ is called the phase of the normal mode. It is clear that for $k>0$, (156) represents a wave moving from right to left(+) or from left to right(-), with speed

[TABLE]

The quantity $v_{f}$ is called the phase speed of the wave. If $\omega^{\prime\prime}(k)=0$ then

[TABLE]

where $c$ is some constant. For this case the phase speed is independent of $k$

[TABLE]

Thus all normal modes moves at the same speed.

For the case when $w^{\prime\prime}(k)\neq 0$, each normal mode moves at its own speed, and the linear superpositions (145) and (146) represents spreading or dispersing waves. For such solutions the phase speed is not the relevant quantity to consider. We will see later that for such waves the group velocity.

[TABLE]

is the relevant quantity.

If the stability index $\Omega\leqslant 0$ and $Re[\lambda(k)]<0$ for all, except a finite number of $k$-values, the equation is said to be of dissipative type.

Let us consider the so-called Telegrapher’s equation

[TABLE]

This is a hyperbolic equation and the polynomial determining $\lambda(k)$ is

[TABLE]

[TABLE]

[TABLE]

Since $\widehat{\lambda}>0,Re[\lambda(k)]<0,\,\forall k\neq 0$. For $k=0$, either $\lambda(0)=0$ or $\lambda(0)=-2\widehat{\lambda}$ so that

[TABLE]

The telegraphers equation is thus neutrally stable. It is also of dissipative type. When $\widehat{\lambda}=0$ we get the regular wave equation which is hyperbolic and conservation. For this reason the telegrapher’s equation is also called the damped wave equation.

The Klein-Gordon equation takes the form

[TABLE]

The equation for $\lambda(k)$ is

[TABLE]

[TABLE]

[TABLE]

The real part of $\lambda(k)$ is zero so the equation is of conservative type, also $\Omega=0$ so it is neutrally stable. The dispersion relation is

[TABLE]

and clearly $\omega^{\prime\prime}(k)=0$ so the Klein-Gorden equation is of dispersive type.

If $\gamma=0$ we get

[TABLE]

and for this case the equation is not of dispersive type.

The Klein-Gordon equation was proposed as a relativistic extension of the Schrodinger equation of quantum mechanics. It was first found by Schrodinger but was rediscovered by Klein and Gordon. In both cases it was rejected as an equation for a relativistic electron on physical grounds. It predicted negative probability! Later it however reappeared in quantum field theory and is now one of the fundamental equations of nature. It for example describes the Higgs Boson.

8 Adjoint differential operators

Let us consider the differential operator

[TABLE]

For any pair of functions $u,w$ we have

[TABLE]

Thus $uLw-wLu$ is a total derivative. Identities like (159) are of outmost importance in the theory of PDEs. Here is another example. Let

[TABLE]

and define a operator $L^{*}$ by

[TABLE]

Then

[TABLE]

Thus we get a total derivative. Next consider an operator

[TABLE]

Then from vector calculus we have

[TABLE]

Thus we get a total derivative. In genera,l if $L$ is a differential operator and if we have, for some vector field $\mathbf{F}$, the identity

[TABLE]

then $L^{*}$ is the formal adjoint of $L$. If $L=L^{*}$ the operator is formally self-adjoint. The operators (159), (163) are formally self-adjoint. The operator $L^{*}=-\frac{\mathrm{d}}{\mathrm{d}x}$ is the formal adjoint of $L=\frac{\mathrm{d}}{\mathrm{d}x}$.

Let $D$ be some domain in $\mathbb{R}^{n}$ and let $L$ be a differential operator defined for functions whose domain include $D$. Let $L^{*}$ be the formal adjoint of $L$. Then by the divergence theorem we have

[TABLE]

It is the possibility of reducing a integral over a domain to a integral over the boundary of the domain that is the reason why the formal adjoint of differential operators are important. We can use this possibility in many important ways. We will return to this topic many times in this and the next class (Mat-3202) in the spring.

For now observe the following:

Let $L=\nabla^{2}\Rightarrow L^{*}=\nabla^{2}$. Then (164) is

[TABLE]

where $\partial_{\mathbf{n}}=\mathbf{n}\cdot\nabla$ is the normal derivative on the surface $S$. Let us assume that the domain of $L$ consists only of functions that satisfy the boundary conditions

[TABLE]

Then (165) implies that

[TABLE]

which can be rewritten into the form

[TABLE]

Let $V$ be the vector space of smooth functions defined on $D$ and that satisfy the boundary conditions (166) (Prove that this is a vector space). On the vector space $V$, we can define and inner product using the formula

[TABLE]

Then (168) show that

[TABLE]

So, $L$ is a self-adjoint operator on $V$. If we choose different boundary conditions so that the right hand side of (165) does not vanish then (169) will not hold and the operator is not self-adjoint in the sense of linear algebra. It is however still formally self-adjoint according to the definition on the previous page.

We know from linear algebra that self-adjoint operators have many good properties, we will discuss many of these in the next section of these lecture notes.

An important class of differential operators in the theory of PDEs consists of operators of the form

[TABLE]

where we use the Einstein convention.

The formal adjoint of this operator is

[TABLE]

To see why this is so take the single term $wa_{11}\partial_{x_{1}}\partial_{x_{1}}u$. We have

[TABLE]

The other terms are treated in the same way and (171) appears.

By carrying out all differentiations in (171) and comparing with (170) it is easy to verify that $L$ is formally self-adjoint only if

[TABLE]

If this is the case we have

[TABLE]

where we have used the fact that $a_{ij}$ is a symmetric matrix. From (172) we see that a self-adjoint operator of the form(170) can be rewritten as

[TABLE]

If the coefficients are constants (173) imply that a self-adjoint operator can not have any first derivative terms.

Let us look at some examples.

The operator

[TABLE]

is an elliptic operator (Show this).

Observe that

[TABLE]

so $L$ is formally self-adjoint.

The operator

[TABLE]

is a hyperbolic operator (Show this).

Assuming that $\rho$ is independent of time we observe that

[TABLE]

where $\widetilde{\nabla}=(\nabla\cdot,\partial_{t})$ is a 4D divergence operator. Thus $\widetilde{L}$ is formally self-adjoint.

Finally define the operator

[TABLE]

where we still assume that $\rho$ is independent of time. This operator is parabolic (Show this). Define another parabolic operator

[TABLE]

Then

[TABLE]

Thus $\widehat{L}^{*}$ is the formal adjoint of $\widehat{L}$. Since $\widehat{L}^{*}\neq\widehat{L}$ the operator $\widehat{L}$ is not formally self-adjoint.

The notion of a formal adjoint can be extended to a much larger class of operators[1] than (170).

9 Initial and boundary value problems in bounded domains

In this section we will introduce the standard approach for solving PDEs in bounded domains. This approach consists of expanding solutions into infinite series of eigenfunctions. The eigenfunctions are determined by the PDE and the boundary values. Eigenfunctions are generalizations of the normal modes introduced in the previous section. In fact eigenfunctions are often called normal modes.

This approach is very general and work for a large range of linear PDEs. We will see, at the end of this section, that they are also very useful when we investigate nonlinear PDEs.

In these lecture notes, we choose a certain level of generality.

Our first restriction is to consider only scalar PDEs. Within the class of scalar PDEs we chose a subclass that is quite general and cover situations that occur in many applications in science and technology. In chapter 8 of [1] a much larger class of PDEs is investigated using similar methods to the ones we discuss here.

Let $L$ be a second order differential operator defined by

[TABLE]

where $\mathbf{x}\in\mathbb{R}^{2}$ or $\mathbb{R}^{3}$. For one space-dimension we define $L$ by

[TABLE]

We will in general assume that $p>0$ and $q\eqslantgtr 0$. The type of equations we want to consider are

Elliptic case:

[TABLE]

Parabolic case:

[TABLE]

Hyperbolic case:

[TABLE]

It is an important point here that the functions $p,\rho$ and $q$ does not depend on time.

These equations, of elliptic, parabolic and hyperbolic type appear as the first approximation to a large set of mathematical modelling problems in science and technology.

Let us consider the case of mass flow. Let $D$ be some domain with boundary $S$.

In $D$ we have a concentration of some material substance given by $u(\mathbf{x},t)$. The total mass of the material inside $D$ is

[TABLE]

where $\rho$ have the interpretation of mass density.

[TABLE]

This requirement explains the assumption $\rho>0$ in equations (176),(177) and (178).

We have the balance equation for mass

[TABLE]

where $\mathbf{j}$ is the mass flux and $H$ is a source of mass. In general, $\mathbf{j}=\mathbf{j}(\mathbf{x},u,\nabla u)$, but for low concentration and small gradients we can Taylor expand and find

[TABLE]

Here a constant term and a term only depending on $u$ has been dropped on physical grounds: All mass flux is driven by gradients, and we have introduced a minus sign and assumed $p>0$ in order to accommondate the fact that mass always flow from high to low concentration. In general $H=H(\mathbf{x},t,u)$ but for low concentration we can also here Taylor expand and find

[TABLE]

For notational convenience we introduce $F$ by

[TABLE]

The minus sign in the first term is introduced to model a situation where high concentration leads to high loss. Imagine a sink hole where the amount leaving increase in proportion to the mass present.

Using (181), (182) and (183) in (180) we get

[TABLE]

[TABLE]

[TABLE]

This is assumed to hold for all domains $D$. Therefore we must have

[TABLE]

where $L$ is the operator defined in (174). In a similar way, using conservation of momentum, equation (178) can be shown to model vibrations in homogeneous solids, membranes and bars.

Our basic model equations must be subject to boundary conditions since the domains in this section are assumed to be bounded. In two and three spatial dimensions we pose the conditions

[TABLE]

where $\alpha,\beta$ and $B$ are given functions defined on the boundary, $\partial G$, of $G$.

In the elliptic case, (176), there is no time dependence and $B=B(\mathbf{x})$. For the parabolic and hyperbolic case $G$ is a space-time domain of the form $G=G_{\mathbf{x}}\times[t_{0},t_{1}]$

For most common situations occurring in applications, both $\alpha$ and $\beta$ are positive. They can not both be zero at the same point since then the boundary condition (184) would not actually be a constraint!

There are 3 special cases of (184).

1. i)

$\alpha\neq 0,\beta=0$ – Diricelet condition or condition of the first kind. 2. ii)

$\alpha=0,\beta\neq 0$ – Neuman condition or condition of the second kind. 3. iii)

$\alpha\neq 0,\beta\neq 0$ – Robin condition or condition of the third kind.

It is not uncommon that conditions of the first, second and third kind are specified at different parts of the boundary.

For heat conduction problems, part of the boundary could have a given temprature (first kind) whereas other parts could be insulated (second kind).

For the one dimensional case, $\partial G$ consists of the two end points of the interval $G=(0,l)$. So $\partial G=\{0,l\}$. For the hyperbolic and parabolic case the general conditions (184) are

[TABLE]

The minus sign at $x=0$ occur because the unit normal at the point $x=0$ is given by the number (-1).

If $F$ and $B$ in (176), (177), (178) and (184) are zero the corresponding boundary value problem is homogeneous.

9.1 Separation of variables

We will now apply the method of separation of variables to the equations and boundary conditions introduced in the previous section.

For the hyperbolic case we consider the homogeneous PDE

[TABLE]

and homogeneous boundary conditions

[TABLE]

in 2D and 3D.

In 1D, the homogeneous condition is

[TABLE]

In order to get a completely determined problem we also must specify initial conditions

[TABLE]

For the parabolic case we consider the homogeneous equation

[TABLE]

with the boundary conditions (187) or (188) and the single initial condition

[TABLE]

In order to give a uniform treatment of the hyperbolic, parabolic and elliptic case, we introduce for the elliptic case a function $u=u(\mathbf{x},y)$ where now $\mathbf{x}$ is in $\mathbb{R}$ or $\mathbb{R}^{2}$. The variable $y$ is scalar and restricted to the domain $0<y<\widehat{l}$. The elliptic PDE can in terms of these variables be written as

[TABLE]

where $\rho,p,q$ are functions of $\mathbf{x}$ alone. The equation (192) is supplied with the boundary conditions (187) for the 2D case or (188) for the 1D case. The analog of the initial conditions in the time dependent case are the following conditions at $y=0$ and $y=l$.

[TABLE]

This assumed form for the elliptic equation certainly does not cover all possibilities that occur in applications. We observe for example that the domain always is a generalized cylinder

However, the assumed form makes a unified presentation for the hyperbolic, parabolic and elliptic case possible.

After all this preparatory work, let us start separating variables!

The method of separation of variables starts by looking for special solutions

[TABLE]

for (186) and (187) and

[TABLE]

for (192). Inserting (194) and (195) into (186), (190) and (192) gives, upon division by $NM$

[TABLE]

Since the two sides of equations (196) $\rightarrow$ (198) depends on different variables, both sides must be equal to the same constant. This constant we write, by convention, in the form $-\lambda$. This give us in all three cases the following equation for $M$.

[TABLE]

and

[TABLE]

The function $M(\mathbf{x})$ in equation (199) is subject to boundary conditions (187) or (188). Since both conditions are homogeneous, as are equation (199), $M(\mathbf{x})=0$ is a solution. This solution, also called the trivial solution, gives $u=0$ and thus plays no role.

Defining $A=\frac{1}{\rho}L$, we observe that the structure of equation (199) is

[TABLE]

Thus $M$ is an eigenvector of $A$ and $\lambda$ the corresponding eigenvalue. The vector space is the set of smooth functions defined on $G$ and satisfying the homogeneous boundary conditions (187) or (188) on the boundary of $G$. (Show that this is a vector space!) In this vector space, vectors are functions, so solutions of (199) are often called eigenfunctions. The big difference between the eigenvalue problems that are studied in elementary linear algebra and the current one, is that here, the vector space is infinite dimensional. This makes the problem much harder, much more interesting and also much more useful.

9.2 Self-adjoint and positive operators

The operator

[TABLE]

occurring in equation (199) and determining the functions $M(\mathbf{x})$ in the separation

[TABLE]

for the hyperbolic case (186), the parabolic case (190) and the elliptic case (192), is very special. We will now determine some of the properties of the spectrum of $L$. The spectrum of a linear operator is the collection of eigenvalues and corresponding eigenspaces of the operator. Recall that the eigenspace corresponding to a given eigenvalue consists of all eigenvectors of the operator that has the same eigenvalue.

Using vector calculus we observe that

[TABLE]

[TABLE]

[TABLE]

[TABLE]

Thus $L$ is formally self-adjoint.

Let $G$ be a domain in $\mathbb{R}^{2}$ or $\mathbb{R}^{3}$. Integrating (204) over $G$ gives

[TABLE]

Let us assume that both $u$ and $w$ satisfy the boundary conditions (184) with $B=0$. Thus

[TABLE]

For a given $u,w$ that satisfy (206) we can think of (206) as a linear system for $\alpha$ and $\beta$

[TABLE]

But our assumptions on the functions $\alpha$ and $\beta$ imply that they can not both be zero at any point. Thus, the system (207) must have a nonzero solution. This imply that the determinant of the matrix defining the linear system (207), must be zero or

[TABLE]

Thus, the right hand side of (205) is zero

[TABLE]

On the vector space, $V$, consisting of real valued functions defined on $G$ and satisfying the homogeneous boundary conditions (206), we define an inner product

[TABLE]

Given this inner product we define the norm of functions in $V$ as

[TABLE]

Assuming $f$ is continuous we have

[TABLE]

Thus, what we have defined is a norm in the usual sense of linear algebra. Using this inner product (209) can be written as

[TABLE]

Thus the operator $\tilde{L}=\frac{1}{\rho}L$ is self-adjoint in the sense of linear algebra. In the most common cases when $\rho=1$, $L$ itself is self-adjoint.

Let now $M_{m}$, $M_{l}$ be eigenfunctions of $L$ corresponding to different eigenvalues, $\lambda_{k}\neq\lambda_{l}$. Since both $M_{k}$ and $M_{l}$ satisfy the homogeneous boundary conditions, (206), we have

[TABLE]

but

[TABLE]

Thus

[TABLE]

But $\lambda_{k}\neq\lambda_{l}\Rightarrow\lambda_{k}-\lambda_{l}\neq 0$, and therefore we must have $(M_{k},M_{l})=0$.

"Eigenfunctions corresponding to different eigenvalues are orthogonal."

Actually we have in our choice of inner product assumed that the eigenfunctions are real valued and the spectrum real. This is in fact true:

Let $\varphi$ be a complex valued function. Then, redoing the argument leading up to (205) using $w=u^{*}$ (the complex conjugate of $u$), we get

[TABLE]

and (208) becomes

[TABLE]

[TABLE]

[TABLE]

Let now $\lambda$ be a complex eigenvalue and $M$ the corresponding complex valued eigenfunction. Then (9.2) imply that

[TABLE]

But

[TABLE]

so we must have

[TABLE]

which means that the eigenvalues are real. The eigenfunctions can then always be chosen to be real.

Let now $u$ be a (real valued) function that satisfies the boundary condition (184) with $B=0$. Thus $u\in V$, as defined on page 108. Then we have

[TABLE]

Split the boundary $\partial G$ into three pieces $S_{1},S_{2}$ and $S_{3}$, where on $S_{i}$ the boundary conditions are of kind $i$, where $i=1,2,3$ (see page 7). Then on $S_{1}\cup S_{2},pu\partial_{\mathbf{n}}u=0$ and on $S_{3}$, $pu\partial_{\mathbf{n}}u=-p\frac{\alpha}{\beta}u^{2}$

[TABLE]

The assumptions we have put on $p,\alpha,\beta$ and $q$ earlier implies that the right hand side of (214) is positive. Thus, using the inner product defined on page 16, we have found that

[TABLE]

This means, by definition, that $\frac{1}{\rho}L$ is a positive operator.

Let now $M_{k}$ be an eigenfunction corresponding to the eigenvalue $\lambda_{k}$. From (215) we have

[TABLE]

and

[TABLE]

Thus we must have

[TABLE]

Therefore, the spectrum of $L$ is not only real but non-negative.

One can prove[1] that there are countably infinite many eigenvalues $\{\lambda_{k}\}_{k=1}^{\infty}$ and that $\infty$ is the only limit point of the sequence of eigenvalues. Recall that a limit point for a sequence $\{a_{n}\}$ is a point $a^{*}\in\mathbb{R}$ such that $\forall\epsilon>0$ there exists infitely many $n$ such that

[TABLE]

Thus the eigenvalues only "accumulate" around $\lambda=\infty$.

Let now $\{\lambda_{k}\}$ be the countably infinite set of eigenvalues. One can also prove that there are at most a finite set of linearly independent eigenvectors corresponding to each $\lambda_{k}$. We can assume that these has been orthogonalized using the Gram-Schmidt orthogonalization process.

Let us in the following for simplicity only consider the case when each eigenspace is of dimension one and where all $k$, $\lambda_{k}>0$. Thus $\lambda=0$ is not an eigenvalue. Then the equations (200-202) for the function $N$ gives upon solution the following:

[TABLE]

Here, $a_{k}$, $b_{k}$ are arbitrary constants, one for each eigenvalue $\lambda_{1},\lambda_{2},\dots$. We now have infinitely many solutions

[TABLE]

that solves equations (186), (190) or (192) and the boundary conditions (187) or (188). These functions will however not in general satisfy the initial conditions. In order to also satisfy the initial conditions (183) for the hyperbolic case, (191) for the parabolic case and (193) for the elliptic case we consider formal superposition

[TABLE]

For the hyperbolic case we get the conditions

[TABLE]

[TABLE]

(220) and (221) require us to find eigenfunction expansions for the the functions $f$ and $g$ in terms of the eigenfunctions $\{M_{k}\}^{\infty}_{k=1}$. In what sense $f$ and $g$ can be represented by their eigenfunction expansions is an important, and in general very hard, question that will be discussed more fully only for the 1D case in this class.

You might recall from elementary calculus that for infinite series of functions, several different notions of convergence are available. The most important ones are pointwise convergence, uniform convergence and mean-square convergence. We will see that the last type of convergence will hold for very mild conditions on $f$ and $g$. However, mean-square convergence is not enough to ensure that (219) is an actual, honest to God, solution to our initial-boundary value problem (186),(187),(189). For this to be the case we need a stronger convergence than mean-square. This stronger form of convergence will allow us to differentiate the formal superposition series (219) term by term. If this is possible, we have a classical solution. If it is not possible, (219) is a generalized solution that might still carry important information about the problem.

In order to express the coefficients $\{a_{k}\}$ and $\{b_{k}\}$ in terms of $f$ and $g$, we use the fact that the eigenfunctions $\{M_{k}\}$ are orthogonal for $k\neq k^{\prime}$

[TABLE]

[TABLE]

In the same way we find

[TABLE]

The formal solution to the hyperbolic initial-boundary value problem (186-188) is then

[TABLE]

where $\{a_{k}\}$ and $\{b_{k}\}$ are given by (222) and (223).

In a similar way, the solution of the initial-boundary problem for the parabolic case is

[TABLE]

where $\{a_{k}\}$ is given by (222).

For the elliptic case we get the conditions

[TABLE]

Proceeding as for the hyperbolic case we get

[TABLE]

These equations can be solved to give unique values for $\{a_{k}\}$ and $\{b_{k}\}$. Our formal solution for the elliptic case is then

[TABLE]

We will now take a closer look at the formal solutions (224), (225) and (227) for the 1D case.

10 The Sturm-Liouville problem and Fourier series

The 1D version of the eigenvalue problem (199) + (187) is

[TABLE]

for $0<x<l$ and boundary conditions

[TABLE]

The problem defined by (228) and (10) is known as the Sturm-Liouville problem. If we require the following additional conditions

[TABLE]

to hold, (228), (10) define a regular Sturm-Liouville problem. If we relax one or more condition on the coefficient functions $p,q$ and $\rho$ we get more generalized versions of the Sturm-Liouville problem. For example, one can allow the functions $p$ and/or $\rho$ to vanish at the endpoints $x=0,x=l$. We then get a singular Sturm-Liouville problem. This generalization do occur in important applications of this theory.

Before we start developing the theory, let us introduce some important constructions. These constructions have already been introduced in the more general setting discussed earlier, but will be repeated here in the 1D context.

For real valued functions on $0\leq x\leq l$ which are bounded and integrable, we introduce the inner product

[TABLE]

Observe that (231) is evidently bilinear and symmetric. For continuous functions we have

[TABLE]

Thus, (231) really defines an inner product in the sense of linear algebra on the linear vector space consisting of continuous functions on $0\leq x\leq l$.

Introduce the norm of a function by

[TABLE]

A function such that $\|\varphi\|<\infty$ is said to be square-integrable. Furthermore, a square-integrable function is normalized to unity if

[TABLE]

If $\varphi$ is any square-integrable function we can normalize this function by defining a new function

[TABLE]

Two functions $\varphi$ and $\psi$ are orthogonal on the interval $0<x<l$ if

[TABLE]

A set of functions $\{\varphi_{k}(x)\},\,k=1,2,\dots$ is said to be an orthogonal set if

[TABLE]

If also $\|\varphi_{k}\|=1$ for all $k$, the set is orthonormal. An orthogonal set can always be turned into an orthonormal set by using the approach described above.

If we allow the functions to assume complex values, we replace the inner product (231) by

[TABLE]

where for any complex number $z$, $\overline{z}$ is its complex conjugate. Observe that we have bilinearity and

[TABLE]

Furthermore for any $\varphi$

[TABLE]

and if $\varphi$ is, say, continuous then

[TABLE]

Thus (236) defines a Hermitian inner product in the sense of linear algebra.

Let us return to the real valued case. Let an orthonormal set of functions $\{\varphi_{k}\}$ be given and let $\varphi$ be a square-integrable function on $0<x<l$. Then the infinite set of numbers $(\varphi,\varphi_{k})$ are called the Fourier coefficients of $\varphi$ with respect to $\{\varphi_{k}\}$. The formal series

[TABLE]

is called the Fourier series of $\varphi$ with respect to the set $\{\varphi_{k}\}$. Note that we use the summation sign in (240) in a formal sense - we have at this point not assumed anything about the convergence of the infinite series of functions (240).

The name, Fourier series, is taken from the special case when $\{\varphi_{k}\}$ is a set of trigonometric functions. Often in the literature the term Fourier series refer only to this trigonometric case, but we will use the term Fourier series for all formal series (240).

So, in what sense does (240) converge? And if it converge, does it converge to the function $\varphi$ that we used to construct the formal series?

Let us start by considering a finite sum,

[TABLE]

For this sum we have

[TABLE]

But the norm of any function is non-negative. Thus we must have

[TABLE]

This inequality says that all partial sums of the numerical series $\sum\limits_{k=1}^{N}(\varphi,\varphi_{k})^{2}$, consisting entirely of positive terms, are bounded by $\|\varphi\|^{2}$. It then follows, from the fact that all increasing bounded sequences are convergent that the series $\sum_{k=1}^{N}(\varphi,\varphi_{k})^{2}$ is convergent and that

[TABLE]

holds. (241) is known as Bessel’s inequality.

A sequence of square-integrable functions $\{S_{N}\}_{N=1}^{\infty}$ is said to converge to a function $\varphi(x)$ in the mean if

[TABLE]

This kind of convergence does not in general imply pointwise convergence. It is called mean square convergence. Observe that if we have equality in Bessel’s inequality (241),

[TABLE]

then we have

[TABLE]

(242) is called Parsevals’s equality. The argument we just did shows that if Parseval’s equality holds, for a given $\varphi$ and given orthonormal set $\{\varphi_{k}\}$, then the Fourier series (240) converges to $\varphi$ in the mean.

An orthonormal set $\{\varphi_{k}\}_{k=1}^{\infty}$ is said to be complete if the Parseval equality (242) holds for all square integrable functions. We have proved that if an orthonormal set of functions is complete then the formal Fourier series (240) for any square integrable functions $\varphi$, converge in the mean square square sense to the function $\varphi$.

10.1 Properties of eigenvalues and eigenfunctions

The properties of the eigenvalues and eigenspaces are the same as for the 2D and 3D case discussed earlier. After all, they follow from structural properties of $L$ and the boundary conditions that are completely analogous in 1D, 2D and 3D. Let us just record these properties here:

1. i)

Eigenfunctions corresponding to different eigenvalues are orthogonal. 2. ii)

The eigenvalues are real and non-negative and the eigenfunctions can be chosen to be real. 3. iii)

Each eigenvalue is simple. This means that each eigenspace is of dimension one. 4. iv)

There is a countable infinity of eigenvalues having a limit point at infinity. The set of eigenvalues can be arranged as

[TABLE]

with $\lim\limits_{k\rightarrow\infty}\lambda_{k}=\infty$. 5. v)

The set of eigenfunctions $\{u_{k}(x)\}_{k=1}^{\infty}$ form a complete orthonormal set of square integrable functions. Thus, for any square integrable function on $0<x<l$, the Fourier series

[TABLE]

converge to $u$ in the mean. 6. vi)

If $u(x)$ is continuous on $0\leq x\leq l$ and has a piecewise continuous first derivative on $0\leq x\leq l$ and $u(x)$ satisfies the boundary conditions (10) then the Fourier series (240) converge absolutely and uniformly to $u(x)$ i $0\leq x\leq l$. 7. vii)

If $u(x)$ has some jump discontinuity at a point $x_{0}$ in the interior of the interval $0\leq x\leq l$ then (240) converge to $\frac{1}{2}(u(x_{0}^{-})+u(x_{0}^{+}))$ at this point.

Here

[TABLE]

so the only thing that remains to do is to actually find the eigenvalues and eigenfunctions for particular cases of interest. We will do so for a few important cases.

All regular Sturm-Liouville problems can be simplified if we use the following uniform approach:

Let $v(x;\lambda),\,w(x;\lambda)$ be solutions of the initial value problem

[TABLE]

with initial conditions

[TABLE]

Then it is easy to verify that the function

[TABLE]

is a solution to (243) and satisfies the boundary condition

[TABLE]

In order to also satisfy the boundary condition at $x=l$ we must have

[TABLE]

The eigenvalues are solutions of equation (246) and the corresponding eigenfunctions are determined by (245).

We will now apply the uniform approach to several cases leading to various trigonometric eigenfunctions and then turn to a couple of important examples of singular Sturm-Liouville problems.

10.1.1 Trigonometric eigenfunctions

For this case we choose $p(x)=\rho(x)=1$ and $q(x)=0$. Equation (9.2) in the Sturm-Liouville problem for this case reduces to

[TABLE]

It is easy to verify that the functions

[TABLE]

satisfy (247) and the initial values (10.1). Here we assume $\lambda>0$. We leave the case $\lambda=0$ as an exercise for the reader. The eigenfunction for $\lambda>0$ is found from (245) and is

[TABLE]

The eigenvalues are determined from the equation (245), which for this case simplifies into

[TABLE]

We will get different sets of eigenvalues and eigenfunctions depending on which boundary conditions (227) we assume for the Sturm-Liouville problem (247).

Fourier Sine series

Assume $\alpha_{1}=\alpha_{2}=1,\beta_{1}=\beta_{2}=0$. Thus the boundary conditions are

[TABLE]

Our equation (251) then reduces to

[TABLE]

which gives us the following nonzero eigenvalues

[TABLE]

and the corresponding normalized eigenfunctions are from (250)

[TABLE]

For this case $\lambda=0$ is not an eigenvalue because for $\lambda=0$ equation (247) reduces to

[TABLE]

The expansion of a function in a series

[TABLE]

is called the Fourier Sine series. Using explicit expressions for the eigenfunction and the inner product we get the following formulas for the Fourier Sine series

[TABLE]

where

[TABLE]

Fourier Cosine series

Assume $\alpha_{1}=\alpha_{2}=0,$ $\beta_{1}=\beta_{2}=1$. Thus the boundary conditions are

[TABLE]

Our equation (251) now reduces to

[TABLE]

$\lambda_{0}=0$ is clearly a solution. For $\lambda_{0}=0$ our equation reduces to

[TABLE]

Thus, $\lambda=0$ is an eigenvalue for this case. The corresponding normalized eigenfunction is

[TABLE]

For $\lambda>0$ we can divide the eigenvalue equation by $\lambda$ and get

[TABLE]

The corresponding normalized eigenfunctions are

[TABLE]

The series expansion

[TABLE]

is the Fourier Cosine series. Using explicit expression for the eigenfunctions and the inner product we get

[TABLE]

where

[TABLE]

Note that all the general properties of eigenvalues and eigenfunctions for Sturm-Liouville problems can be verified to hold for the eigenvalues and eigenfunctions we have found for the Fourier Sine and Fourier Cosine cases.

Fourier series

In this case we consider a slightly generalized Sturm-Liouville problem. We assume that the domain is $(-l,l)$ instead of $(0,l)$. We still make the assumptions $p(x)=\rho(x)=1$, $q(x)=0$, so the eigenvalue equation is

[TABLE]

Instead of using the boundary conditions (227) we use periodic boundary conditions

[TABLE]

The general solution of (261) is for $\lambda>0$

[TABLE]

and

[TABLE]

[TABLE]

The eigenvalues $\lambda>0$ are then clearly

[TABLE]

For each such eigenvalue $a$ and $b$ are arbitrary. Thus (264) gives us an eigenspace of dimension 2 spanned by $\cos\sqrt{\lambda_{k}}x$ and $\sin\sqrt{\lambda_{k}}x.$ Introducing an inner product and corresponding norm

[TABLE]

we can normalize the eigenfunctions and find

[TABLE]

One can show that there are no negative eigenvalues and that $\lambda_{0}=0$ is an eigenvalue with corresponding normalized eigenfunction

[TABLE]

Using these functions we get a Fourier series

[TABLE]

Using experssions for the eigenfunctions and inner product we get

[TABLE]

where

[TABLE]

If $u$ is continuous, has a piecewise continous derivative in the interval $[-l,l]$ and satisfies the boundary conditions (262), then the series (270) converges absolutely and uniformly to $u(x)$.

10.1.2 Bessel eigenfunctions

We are now going to consider a singular Sturm-Liouville problem that arise from the eigenvalue problem for Bessel’s equation of order $n$

[TABLE]

where $n\geq 0$ is an integer. This is a Sturm-Liouville problem with $p(x)=\rho(x)=x$, $q(x)=\frac{n^{2}}{x}$. The boundary conditions of interest for (271) are

[TABLE]

Requiring that $u(x)$ is bounded at $x=0$ is a significant condition here because (271) is a singular Sturm-Liouville problem. We observe that $p,\,\rho$ and $q$ are positive on $0<x<l$ but that $p$ and $\rho$ vanish at $x=0$ and that $q$ is singular at $x=0$. We might therefore in general expect there to be solutions that are singular at $x=0$ and (272) eliminate such solutions as acceptable eigenfunctions. One might worry that all solutions are singular at $x=0$. Then we will get no eigenvalues and eigenfunctions at all! By using the method of Frobenious from the theory of ordinary differential equations one can show that (271) has solutions that are bounded at $x=0$. These solutions can be expressed in terms of Bessel functions of order $n$, $J_{n}$, in the following simple way

[TABLE]

We find the eigenvalues by solving the equation

[TABLE]

One can prove that this equation has infinitely many solutions $\lambda_{kn},\,k=1,2,3,\dots$ They are all positive and $\lambda_{kn}\rightarrow\infty$ as $k\rightarrow\infty$. A large number of zeros for the function $J_{n}$ has been found using numerical methods. Let these zeros be $\alpha_{kn},\,k=1,2,3,\dots$ In terms of these zeros the eigenvalues are

[TABLE]

The corresponding eigenfunctions are

[TABLE]

Using the inner product

[TABLE]

the eigenfunctions can be normalized

[TABLE]

where we have used the formulas

[TABLE]

As you have seen in calculus there are a large number of formulas linking elementary functions like $\sin(x)$, $\cos(x)$, $\ln(x)$, $e^{x}$ etc. and their derivatives and anti-derivatives. The Bessel functions $J_{n}$ is an extension of this machinery. There are a large set of formulas linking the $J_{n}$ and their derivatives and anti-derivatives to each other and other functions like $\sin(x)$, $\cos(x)$, $\ln(x)$ etc. The functions in calculus ($\sin(x)$, $\cos(x),\dots$) are called elementary functions. The function $J_{n}$ is our first example of a much larger machinery of functions called special functions. There are many handbooks that list some of the important relations and formulas involving special functions. Symbolic systems like mathematica, maple, etc. embody most of these formulas and often produce answers to problems formulated in terms of such functions.

Note that the boundary conditions (272) are not of the type (227) discussed in the general theory of Sturm-Liouville theory. Thus, there is no guarantee at this point that eigenfunctions corresponding to different eigenvalues are orthogonal. This is typical for singular Sturm-Liouville problems; the properties of the eigenvalues and eigenfunctions must be investigated in each separate case.

Observe that, in general, with $L$ defined as in (9.2) we have

[TABLE]

[TABLE]

Recall that for our current singular Sturm-Liouville problem, the space of functions we need to expand in terms of the eigenfunctions (277), consists of functions such that

[TABLE]

Thus on the right hand side of (279) the evaluation at $x=0$ must be interpreted as a limit as $x$ approach [math] from inside the interval $(0,l)$.

Let us apply (279) to a pair of eigenfunctions $u_{i},u_{j}$ corresponding to different eigenvalues $\lambda_{i}\neq\lambda_{j}$. Thus

[TABLE]

But then (279) gives us

[TABLE]

In order for eigenvectors corresponding to different eigenvalues to be orthogonal, the right hand side (281) must be zero. The evaluation at $x=l$ is zero since the eigenfunctions are zero at this point. In order for the right hand side of (281) to vanish, it is in general not enough to require that the eigenfunctions are bounded at $x=0$. The derivative of the eigenfunctions might still be unbounded at $x=0$ in such a way that the right hand side of (281) is nonzero and thus ruining the orthogonality of the eigenfunctions. However, if we also require that

[TABLE]

everything is okay and the eigenfunctions form an orthogonal set. Since $p(x)=x$ we thus have the requirement

[TABLE]

From the theory of Bessel functions we have

[TABLE]

Thus the functions are bounded at $x=0$. Also

[TABLE]

but then

[TABLE]

Thus, by the argument given, the eigenfunctions

[TABLE]

form an orthonormal set.

One can show that these eigenfunctions form a complete set. Any smooth function defined on the interval $0\leq x\leq l$ that satisfies the boundary condition $u(l)=0$ can be expanded in a convergent series

[TABLE]

10.1.3 Legendre polynomial eigenfunctions

Another singular Sturm-Liouville problem that occur frequently is

[TABLE]

with boundary conditions

[TABLE]

The function $p(x)=1-x^{2}$ vanishes at both endpoints $x=\pm 1$ so we indeed have a singular Sturm-Liouville problem.

One can use the method of Frobenius, alluded to earlier, to show that the only bounded solutions to (284) are in fact polynomials and that these exists only for

[TABLE]

The normalized eigenfunctions, $P_{k}$, are called the Legendre polynomials of degree $k$

[TABLE]

It is very easy to verify that

[TABLE]

so this singular Sturm-Liouville problem has an orthogonal set of eigenfunctions according to previous arguments. One can show that

[TABLE]

Thus we get an orthonormal set of eigenfunctions

[TABLE]

One can show that this set is complete and that any smooth function defined over the interval $-1\leq x\leq 1$ can be expanded in a convergent series

[TABLE]

10.2 Series solutions to initial boundary value problems

It is now time for returning to the solution of initial boundary value problems. We will consider a few specific examples and do detailed calculations for these cases. However, before we do that let us recapitulate the separation of variable strategy for the hyperbolic case. As we have seen, the strategy for parabolic and elliptic equations is very similar.

1. i)

Look for functions of the form

[TABLE]

that solves the PDE (186) and the boundary conditions. These could be of the form (187), but as we have seen in the previous section on singular Sturm-Liouville problems, other possibilities exist. 2. ii)

If the boundary condition has been choosen in an appropriate way we find that $M$ must be one of a complete set of eigenfunctions $\{M_{k}\}_{k=1}^{\infty}$ of a self-adjoint eigenvalue problem with eigenvalues $\{\lambda_{k}\}_{k=1}^{\infty}$ corresponding to the eigenfunctions $\{M_{k}\}_{k=1}^{\infty}$. 3. iii)

The function $N(t)$ is found to be one of

[TABLE]

where $\{a_{k}\},$ $\{b_{k}\}$ are constants that can be freely specified. 4. iv)

A formal solution to the PDE and the boundary conditions is

[TABLE] 5. v)

Let $a_{k}$ and $b_{k}$ be the Fourier coefficients

[TABLE]

where $f$ and $g$ are the initial data. Then

[TABLE]

and thus the formal solution (294) satisfies the PDE, the boundary conditions and the initial conditions. 6. vi)

Investigate the convergence of the series (294). If the initial data $f$ and $g$ are smooth enough the series will converges so fast that it can be differentiated twice termwise. In this case the formal series solution is a classical (honest to God!) solution to our initial/boundary value problem. For less smooth $f$ and $g$ the convergence is weaker and will lead to generalized (non-classical) solutions.

10.2.1 Vibrations of a string with fixed endpoints

We have seen earlier in this class that a vibrating string with fixed endpoints can in certain cases be modelled using the equation

[TABLE]

[TABLE]

where $c=\sqrt{\frac{T}{\rho}}$, and where $T$ and $\rho$ are the, assumed constant, tension and mass density of the string. In order to get a unique motion of the string, we must specify the displacement and movement of the string at $t=0$

[TABLE]

These are the initial values. We now separate variables, $u(x,t)=M(x)N(t)$, and get the following regular Sturm-Liouville problem for $M$

[TABLE]

This is the problem we investigated on page 122. We found eigenvalues and normalized eigenfunctions of the form

[TABLE]

For the function $N(t)$ we find

[TABLE]

and the corresponding separated solutions $u_{k}(x,t)$ are

[TABLE]

In order to satisfy the initial conditions we write down the formal solution

[TABLE]

Upon applying the initial conditions (299) we get

[TABLE]

These are formally satisfied if

[TABLE]

With these values for $a_{k}$ and $b_{k}$, (304), is the formal solution of the initial value problem for the vibrating string fixed at both ends.

The separated solutions (303) are called normal modes for the vibrating string. The formal solution (304) is a sum, or linear superposition, of normal modes. The notion of normal modes is a concept that goes way beyond their use for the vibrating string. Separated solutions satisfying a hyperbolic PDE and boundary conditions are in general called normal modes.

Normal modes is a physical concept that has very wide applicability and deep roots. Physicists in general are not very familiar with eigenvalue problems for linear differential operators and their associated eigenfunctions, but they all know what normal modes are!

In quantum theory applied to fields, normal modes is a key idea. In this setting, they are in fact the mathematical embodiment of elementary particles. The newly discovered Higgs Boson is a normal mode for a certain field equation closely related to the 3D wave equation. In the quantum theory of materials, so-called many-body theory, the interaction between material normal modes (phonons) and electronic normal modes is responsible for the phenomenon of superconductivity.

Now then, when can we expect the formal solution (304) to be an actual, honest to God, classical solution? In general, finding sharp conditions on the initial data $f(\mathbf{x})$ and $g(\mathbf{x})$ that will imply the existence of a classical solution is a hard mathematical problem that is best left to specialists in mathematical analysis. In fact, sometimes the problem is so hard that entirely new mathematical ideas have to be developed in order to solve this problem.

You might know that all of modern mathematics is formulated in terms of, and gets its consistency from, set theory. You might not know that set theory was invented by George Cantor in the 1880’s while studying the convergence of trigonometric Fourier series!

So what does mathematical analysis tell us about the convergence of the formal solution (304)? The following conditions are sufficient to ensure the existence of a classical solution for the vibrating string:

1. i)

$f$ must have two continuous derivatives and a third piecewise continuous derivative, 2. ii)

$f(0)=f(l)=0,\quad f^{\prime\prime}(0)=f^{\prime\prime}(l)=0$, 3. iii)

$g$ must have one continuous derivative and a second piecewise continuous derivative, 4. iv)

$g(0)=g(l)=0.$

By using standard trigonometric formulas, the normal modes (303) can be rewritten as

[TABLE]

where

[TABLE]

The formula (306) shows that for a normal mode each point $x_{0}$ perform a harmonic vibration with frequency

[TABLE]

and amplitude

[TABLE]

Solutions to the wave equation are in general called waves. The normal mode (303) then is a standing wave; each point on the string oscillate in place as $t$ varies.

Some points on the string remain fixed during the oscillation. These are the points

[TABLE]

which are called nodes for the standing wave. The standing wave has maximum amplitude at the points

[TABLE]

which are called antinodes.

Since each normal mode performs a harmonic oscillation, $u_{k}(x,t)$ is often called the $k^{th}$ harmonic. Since a vibrating string makes sound, the normal mode $u_{k}$ will produce a pure tone at frequency $\omega_{k}$. Thus the decomposition of the string into normal modes corresponds to the decomposition of sound into pure tones. If the tone-spectrum is pleasing we call it music!

Recall that the energy of a small piece of a vibrating string located between $x$ and $x+\mathrm{d}x$ is

[TABLE]

Thus, the total energy of the vibrating string is

[TABLE]

Writing

[TABLE]

where

[TABLE]

we observe that

[TABLE]

In a similar way we find that

[TABLE]

Therefore

[TABLE]

where

[TABLE]

where $m=\rho l=\text{ mass of string}$.(Do this calculation!)

Thus each $E_{k}$ is independent of $t$ and so is the total energy

[TABLE]

which is the sum of the energies, $E_{k}$, of each normal mode.

The decomposition of the total energy of a system into the energy of each normal mode plays a key role in the quantum theory of fields and in many-body theory.

Using trigonometric addition formulas one can show that

[TABLE]

Using these formulas on the expression (303) for the normal modes we find

[TABLE]

Thus the normal mode can be written as a sum of left and right traveling waves. This is not surprising since we know that all solutions of the wave equation can be written as a sum of left and right traveling waves. What is special for the normal modes is that their traveling wave components are arranged in such a way that we get a standing wave.

10.2.2 Heat conduction in a finite rod

As we have seen previously in this class, heat conduction in a finite rod, where there are no heat sources and where the endpoints of the rod are kept at zero temperature, is in certain cases modeled by the boundary value problem

[TABLE]

where $u(x,t)$ is the temperature of the rod as a function of time, and where $c^{2}$ is the coefficient of heat conductivity. In order to get a unique solution we must specify the initial temperature distribution

[TABLE]

The separation of variable approach leads to the same eigenvalue problem as for the string. The corresponding functions $N_{k}(t)$ now become

[TABLE]

and we get a formal solution

[TABLE]

The initial value (313) implies that the constants $a_{k}$ must be

[TABLE]

(317) with the constants given by (318) is a formal solution to our initial boundary value problem.

Because of the exponential term in the solution for $t>0$, the series (317) converges so strongly that the series can be differentiated termwise as many times as we like. Thus (317) satisfies the PDE and boundary conditions if we merely assume that $f(x)$ is a bounded function on $0\leq x\leq l$. However, as $t\to 0$ the exponential term vanishes, so in order to ensure that (317) also satisfies the initial condition, we need stronger conditions on $f$. Our friendly mathematical analyst tells us that the following conditions are sufficient:

1. i)

f must be continuous and f’ piecewise continuous in $0\leq x\leq l$ 2. ii)

f(0)=f(l)=0

For the wave equation we know that discontinuities in the initial data are preserved by the equation and propagated along the characteristic curves $x\pm ct=\text{const}$. For the heat conduction problem, the exponential factor will ensure that the solution is infinitely many times differentiable for any $t>0$ even if the initial data is discontinuous!

The heat conduction equation can not support any kind of discontinuities for $t>0$. We have seen this same result appear previously from the fact that the only characteristic curves for the heat conduction equation are of the form $t=\text{const}$.

10.2.3 The Laplace equation in a rectangle

Consider the boundary value problem

[TABLE]

The equation, which is the Laplace equation, and boundary conditions determine the steady state displacement of a stretched membrane on a rectangular domain. The functions $f$ and $g$ specify the stretching of the membrane on the boundary curves $y=0$ and $y=\hat{l}$.

The separation of variable method leads to the same eigenvalue problem as before but now the corresponding functions $N_{k}$ are

[TABLE]

The separated solutions then are

[TABLE]

where we have introduced hyperbolic functions instead of exponential ones.

The formal solution to our equation and homogeneous boundary conditions are then

[TABLE]

The boundary conditions at $y=0$ and $y=\hat{l}$ imply that

[TABLE]

These are formally satisfied if

[TABLE]

(322) with $a_{k}$ and $b_{k}$ given by (325) determines the formal solution to the problem (319).

Assume that

[TABLE]

Then

[TABLE]

and using the fact that for $z\geq 0$

[TABLE]

we get the estimate

[TABLE]

Thus, this term is bounded by an exponential as $k\to\infty,\medspace\forall\,0<y<\hat{l}.$ In a similar way $|a_{k}\sinh\left(\frac{\pi k}{l}y\right)|$ is bounded by a decaying exponential as $k\to\infty$. Therefore the terms in the formal solution (322) decay exponentially as $k\to\infty$ for all $0<x<l,0<y<\hat{l}$. The formal series (322) therefore converges and defines a function in the inside of the rectangle $[0,l]\times[0,\hat{l}]$ whose series representation can be differentiated termwise as many times as we like. Thus the formal series (322) defines a classical solution to the Laplace equation inside a closed rectangle $[0,l]\times[0,\hat{l}]$. It, by construction of the separated solutions, satisfies the boundary conditions on the edges $x=0,\,x=l,\,0\leq y\leq\hat{l}$. As we approach $y=0$ or $y=\hat{l}$, the exponential decay is weakened and in order to ensure that the formal solution satisfies the boundary conditions on $y=0,\,\hat{l},\,0\leq x\leq l$ our friendly mathematical analysis tells us that it is sufficient that:

1. i)

$f$, $g$ are continuous and $f^{\prime}$, $g^{\prime}$ piecewise continuous in the interval $0\leq x\leq l$ 2. ii)

$f(0)=f(l)=g(0)=g(l)=0$

The fact that solutions to the Laplace equations are smooth in the interior of the rectangle even if the boundary data has discontinuities is, as we have seen, also a consequence of having no real characteristics at all.

10.3 Inhomogenous equations - Duhamel’s principle

So far, we have only considered the separation of variable method for homogeneous equations. This means physically that we are modeling source-free problems. It is now time to consider the presence of sources. We will here focus on the hyperbolic and parabolic case:

[TABLE]

The equations (329) are supplied with homogeneous boundary conditions and also homogeneous initial conditions

[TABLE]

The initial/boundary value problem (329),(330) is solved using Duhamel’s principle. We proceed as follows:

Let $v(x,t;\tau)$ be a function that satisfies

[TABLE]

with initial conditions

[TABLE]

We assume this solution has been found using separation of variables or some other method.

Define a function $u(\mathbf{x},t)$ by

[TABLE]

Then we have by the rules of calculus

[TABLE]

and from the initial conditions

[TABLE]

Thus we have for the hyperbolic case

[TABLE]

(10.3) show that the $u(\mathbf{x},t)$ constructed in (332) solves the inhomogeneous problem for the hyperbolic case. In a similar way we can show that the function $v(\mathbf{x},t;\tau)$ for the parabolic case also solves the ihomogenous equation for that case. (332) is Duhamel’s principle.

10.3.1 The inhomogeneous wave equation

Let the following initial/boundary value problem be given

[TABLE]

According to Duhamel’s principle we should construct a function $v(x,t;\tau)$ which solves

[TABLE]

We solve (10.3.1) using the d’Alembert formula

[TABLE]

The solution to (336), which we denote by $u_{p}$, is then according to Duhamel’s principle

[TABLE]

Let $u_{h}(x,t)$ be the solution to the homogeneous wave equation with arbitrary initial data

[TABLE]

According to the d’Alembert formula the solution with initial data (340) is

[TABLE]

Let now

[TABLE]

Verify that $u(x,t)$ satisfies the inhomogeneous wave equation with initial data (340). Inserting the formulas for $u_{p}$ and $u_{h}$ we get

[TABLE]

From this formula it is clear that the value of $u$ at a point $(x_{0},t_{0}$), with $t_{0}>0$ depends only on the data at points in a characteristic triangle.

The characteristic triangle is called the domain of dependence of the solution at the point $(x_{0},t_{0})$.

Similarly, we define the domain of influence of a point $(x_{0},t_{0})$ to be the set of space-time points $(x,t)$ such that the domain of dependence of the solution at $(x,t)$ includes the point $(x_{0},t_{0})$.

The point $(x_{0},t_{0})$ will have no causal connection to points beyond its domain of influence. These notions play a fundamental role in the theory of relativity, where the domain of influence of a space-time point $(x,t)$ is called a light cone. All material particles and force mediating particles must move within the light cone. This is strictly true only for classical physics. In quantum theory, both material particles and force mediators can, in the form of virtual particles, move beyond the light cone.

10.3.2 The inhomogeneous heat equation

Let us consider the 1D heat equation on a finite interval, $0<x<l$ , driven by a heat source $g(x,t)$ located in the interval. Assuming the endpoints are kept at zero temperature, our problem is

[TABLE]

The function $u(x,t;\tau)$ appearing in Duhamel’s principle satisfies (344) with $g$ replaced by zero and the initial condition

[TABLE]

Separation of variables for this problem gives (page 72)

[TABLE]

where

[TABLE]

Duhamel’s principle now gives the solution to (344) in the form

[TABLE]

Let us consider the special case

[TABLE]

It is easy to verify that for this choice

[TABLE]

[TABLE]

Observe that

[TABLE]

when $t\to+\infty$. By direct differentiation we find that the limiting function, $u(x)$, is a solution to

[TABLE]

which shows that the function $u(x)$ is a time independent, or stationary solution to the heat equation.

In fact, if $g(x,t)=G(x)$ then from (348) we get

[TABLE]

For this case the limiting function $u(x)$ is

[TABLE]

It is easy to verify that $u(x)$ is in fact a stationary solution to the driven heat equation (344)

[TABLE]

10.4 The finite Fourier transform

Duhamel’s principle can, as we have seen, be used to construct solutions to inhomogeneous equations. However, the solutions are expressed as an integral over the parameter $\tau$ and are rarely explicit solutions. In order to get actual numbers out we typically have to do numerical integration. This can quickly become awkward and expensive to do in realistic situations.

There is however another approach that is more general and which is better behaved numerically. This is the finite Fourier transform. In this method we expand the solutions in terms of eigenfunctions of the spatial part of the differential operator defining the equation. The expansion coefficients are time dependent and the PDE is reduced to an infinite set of coupled ODEs. We have thus transformed the PDE into a system of ODEs. This process is called the finite Fourier transform. The infinite series of eigenfunctions expressing the solution of the PDE in terms of the time dependent coefficients is called the inverse finite Fourier transform.

We have thus the following solution strategy:

This transform strategy is an important approach to solving complex problems in applied mathematics and theoretical physics.

We will stick to the restricted class discussed previously in these notes. For that class we have an operator, $L$, given by

[TABLE]

We can give a unified presentation of the finite Fourier transform if we introduce an operator $K$ which is

[TABLE]

In fact, if we also allow $K=0$, the general 3D elliptic problem can be included. With this definition of $K$ our problem of interest is

[TABLE]

Let $\{M_{k}\}$ be the complete set of eigenfunctions for the eigenvalue problem

[TABLE]

where as usual we apply the spatial boundary conditions. We will assume that the eigenfunctions have been normalized so that

[TABLE]

We now expand the solutions of (361) in terms of the eigenfunctions $M_{k}$

[TABLE]

where

[TABLE]

The coefficients $N_{k}$ are functions of $t$ for the hyperbolic and parabolic case and a function of $y$ for the elliptic case with (generalized) cylinder symmetry. For the general elliptic case the $N_{k}$’s are constants.

The formula (365) defines the finite Fourier transform and (364) the inverse finite Fourier transform.

In order to derive an equation for the coefficients $N_{k}$ we multiply (361) by $M_{k}(\mathbf{x})$ and integrate over the domain for the variable $\mathbf{x}$. We let this domain be denoted by $G$.

This integration gives us the equation

[TABLE]

Using the definition of the finite Fourier transform and the inner product, each term in this equation can be reformulated into a more convenient form.

For the first term we have

[TABLE]

For the second term in(366) we have

[TABLE]

The second term in (10.4) can be further simplified using the boundary conditions.

Recall that the boundary conditions on $u$ and $M_{k}$ are

[TABLE]

where $u$ and $B$ are functions of $\mathbf{x}$ and $t$ for the hyperbolic and parabolic case, functions of $\mathbf{x}$ and $y$ for the restricted elliptic case and a function of $\mathbf{x}$ for the general elliptic case. As usual $\alpha=\alpha(\mathbf{x})$ and $\beta=\beta(\mathbf{x})$.

Using the classification of the boundary conditions into type $1,2,3$ introduced on page 104 we have for the second term in (10.4)

[TABLE]

Substituting (367), (10.4) and (10.4) into (366) gives us the system

[TABLE]

where

[TABLE]

For the general elliptic case $K=0$ so (371) are algebraic equations for the $N_{k}$. For the restricted elliptic case we have a system of ordinary differential equations

[TABLE]

This system of ODEs must be solved as a boundary value problem. The boundary conditions are

[TABLE]

where we have used the boundary conditions and where $u(0)\equiv u(\mathbf{x},t)$ etc.

For the parabolic case $K=\partial_{t}$ we have a system of ODEs

[TABLE]

subject to the initial condition

[TABLE]

Observe that (374) is an infinite system of ODEs, but they are uncoupled. This is a result of the linearity of the underlying PDEs and the fact that the $M_{k}$’s are eigenfunctions for the operator $L$ defining the equations. The finite Fourier transform can also be used on nonlinear PDEs, as we will see, and for this case the ODE’s for the $N_{k}(t)$ will be coupled. The ODEs (374) are of a particular simple type. They are linear first order equations. Such equations can be solved using an integrating factor. The solution is found to be

[TABLE]

This formula is only useful if we can get some exact or approximate analytical solution of the integral. If a numerical solution is required it might be more efficient to solve the ODEs (374) numerically, using for example Runge-Kutta.

For the hyperbolic case, $K=\partial_{tt}$ and we get a system of second order ODEs

[TABLE]

subject to the initial conditions

[TABLE]

This second order equation is also of a simple type which is solvable by the method of variation of parameters. We get the formula

[TABLE]

Just as for the parabolic case the formula (379) is only useful if some exact or approximate analytical solution can be found. If a numerical solution is required it is more efficient to solve the ODEs (377) using a numerical ODE solver.

Recall that when the $N_{k}$s have been computed, the solution to the PDE is constructed using the inverse finite Fourier transform

[TABLE]

The representation (380) is problematic when we approach the boundary since all the $M_{k}$s satisfy homogeneous boundary conditions, whereas $u$ satisfies inhomogeneous boundary conditions. If these conditions are of, say, type 1 and if (380) is a classical solution (series converges pointwise) then (380) implies that $u$ is zero on the boundary. But we know that $u$ in general satisfies nonzero boundary conditions. Thus the series (380) can not converge pointwise as we approach the boundary. Therefore, (380) is in general not a classical solution, but rather some form of generalized solution.

The series might still converge pointwise away from the boundary, but the convergence will typically be slow since the expansion coefficients, $N_{k}$, "know about" the function over the whole domain and this leads to slow convergence everywhere.

There is however a way, which we will discuss shortly, to avoid this boundary problem and get series that converge fast enough to get a classical solution.

Let us now look at a couple of examples.

10.4.1 Hyperbolic equations

We consider the following hyperbolic PDE

[TABLE]

where $M_{i}(\mathbf{x})$ are one of the eigenfunctions of $L$. We assume that initial and boundary conditions are homogeneous

[TABLE]

For this case we have

[TABLE]

and because the initial conditions for the PDE (381) are homogeneous, we have

[TABLE]

Thus our system of ODEs simplifies into

[TABLE]

with initial conditions

[TABLE]

We get $N_{k}=0$, $k\neq i$ and from (379)

[TABLE]

The inverse finite Fourier transform gives us

[TABLE]

This is now the exact solution to this problem. From (387) it is clear that as $\omega\to\sqrt{\lambda_{i}}$ the amplitude for the solution (387) becomes larger. We can not evaluate (387) at $\omega=\sqrt{\lambda_{i}}$, since the denominator is then zero, but if we take the limit of (387) as $\omega\to\sqrt{\lambda_{i}}$, using L’Hopitals rule on the way, we find

[TABLE]

The formulas (387) and (388) tell the following story: The general solution of the homogeneous equation

[TABLE]

is a sum of normal modes $u_{k}(\mathbf{x},t)$, whose oscillation frequencies are $\omega_{k}=\sqrt{\lambda_{k}}$. If we drive (386) with a source that oscillates at a frequency $\omega$ and whose spatial profile is an eigenfunction, $M_{i}$, the solution is a linear superposition of terms oscillating at $\omega$ and $\omega_{i}$, with a common spatial profile $M_{i}(\mathbf{x})$. As $\omega$ approaches $\omega_{i}$, the amplitude of the oscillation grows and in the limit $\omega\to\omega_{i}$ the system displays linear unbounded growth in time.

What we have here is the well known phenomena of resonance from particle mechanics. Driven continuum system modeled by PDEs can display resonance. Understanding, predicting and preventing resonance is of major importance in applied mathematics and physics. After all we don’t want to end up in a plane whose jet engines fall off because resonance has weakened the bolts tying the engine to the plane fuselage!

10.4.2 The Poisson equation in a disk

We consider the boundary value problem for Poisson’s equation

[TABLE]

For this geometry it makes sense to use polar coordinates

[TABLE]

Then

[TABLE]

Thus in polar coordinates (390) can be written as

[TABLE]

By finding the eigenfunctions, $M_{k}(r,\theta)$, of the 2D elliptic operator

[TABLE]

we can express the solution $u$ as a Fourier series

[TABLE]

This is done in chapter 8 of our text-book. Here we will solve this problem using the finite Fourier transform. Introduce the operators

[TABLE]

Then (392) can be written as

[TABLE]

The idea is now to introduce a finite Fourier transform based on the eigenfunctions of the operator $\hat{L}$. Observe that since $\theta$ is the polar angle, running from zero to $2\pi$, all functions are periodic in $\theta$ with period $2\pi$. Thus we are lead to the eigenvalue problem

[TABLE]

This problem is easiest to solve using complex exponentials. The general solution of the equation (398) is for $\lambda>0$

[TABLE]

and the periodicity condition (399) implies that

[TABLE]

The functions $e^{\pm i\sqrt{\lambda}\theta}$ are linearly independent (verify this), so we must conclude from (401) that

[TABLE]

If $\left(1-e^{-i\sqrt{\lambda}2\pi}\right)\neq 0$ both $A$ and $B$ must be zero. This gives $M(\theta)=0$, so such $\lambda$s are not eigenvalues. The other possibility is that

[TABLE]

and then $A$, $B$ are arbitrary. The equation (403) is only satisfied if

[TABLE]

These are the positive eigenvalues. The corresponding eigenspaces are two-dimensional and spanned by the normalized eigenfunctions

[TABLE]

In order to find out if $\lambda=0$ is an eigenvalue we put $\lambda=0$ in (398) and get the equation

[TABLE]

whose general solution is

[TABLE]

The boundary condition now gives

[TABLE]

and $B$ is arbitrary, so the eigenspace corresponding to $\lambda=0$ is one-dimensional and spanned by the normalized eigenfunction

[TABLE]

Verify that there are no negative eigenvalues.

Using (405) and (408) we define a finite Fourier transform for functions $u(r,\theta)$ by

[TABLE]

The corresponding inverse Finite Fourier transform is

[TABLE]

The ODEs for the coefficients $a_{k}$ and $b_{k}$ are now found by multiplying (397) by the eigenfunctions and integrating over the interval $(0,2\pi)$. Verify that we get

[TABLE]

where $A_{k}$ and $B_{k}$ are the Fourier coefficients of $F$,

[TABLE]

Note that only solutions of the singular equations (10.4.2) that are bounded as $r\to 0$ are acceptable. So our boundary condition at $r=0$ is

[TABLE]

At $r=R$ the boundary condition is

[TABLE]

Taking the Finite Fourier transform of the boundary data $f(\theta)$, we find the the boundary conditions for the functions $a_{k}(r),b_{k}(r)$ at $r=R$ aret

[TABLE]

Note that both of the equations in (10.4.2) are inhomogeneous versions of the homogeneous Euler equation

[TABLE]

A basis for the solution space of such equations can be found using standard methods. We have

[TABLE]

The solutions of (10.4.2) can now be found using variation of parameters. We have

[TABLE]

This completes the formal solution of (390).

Let us now return to the general scheme. As we have indicated, the solutions found using the finite Fourier transform will in general have slow convergence because the eigenfunctions satisfy the homogeneous boundary conditions, whereas in general the solution might satisfies inhomogeneous boundary conditions.

However, the problem of slow convergence can in general be solved using an approach which we will now explain.

Let us first consider the case when the inhomogeneous terms in the equation and boundary conditions are independent of time (we are here thinking about the hyperbolic and the parabolic cases).

Thus we have

[TABLE]

Let $v(\mathbf{x})$ be a time independent solution to (10.4.2). Thus

[TABLE]

We will assume that we can solve (423). Define $w(\mathbf{x},t)=u(\mathbf{x},t)-v(\mathbf{x})$. Then $w(\mathbf{x})$ solves

[TABLE]

and has boundary condition

[TABLE]

Thus $w(\mathbf{x},t)$ satisfies a homogeneous equation and homogeneous boundary conditions. The formal series for $w(\mathbf{x},t)$ will therefore have better convergence properties than the formal series for $u(\mathbf{x},t)$ . If $v(\mathbf{x})$ can be found without using eigenfunction expansions, the solution $u(\mathbf{x},t)$ will be represented by a Fourier series with good convergence properties.

Finding an expression for $v(\mathbf{x})$ might however not be easy.

There is however a related approach that can be used and that do not rely on finding a stationary solution.

Let us assume that we can find a function $v(\mathbf{x},t)$ such that

[TABLE]

The function $v$ must be sufficiently smooth so that it is possible to substitute into the equation

[TABLE]

Define a function $w(\mathbf{x},t)=u(\mathbf{x},t)-v(\mathbf{x},t)$. Then

[TABLE]

and

[TABLE]

where

[TABLE]

Thus $w(\mathbf{x},t)$ solves an inhomogeneous equation with homogeneous boundary conditions. The solution to this problem can be found using Duhamel’s principle or the finite Fourier transform. Convergence of the series should now be much faster and thus lead to more efficient numerical calculations.

Let us illustrate this method for a 1D inhomogeneous heat equation.

[TABLE]

It is easy to see, using linear interpolation that the function

[TABLE]

satisfy the inhomogeneous boundary conditions. Defining $w(x,t)=u(x,t)-v(x,t)$, we find immediately that

[TABLE]

for $t>0$ and $0<x<l$. The boundary condition are by construction homogeneous

[TABLE]

and the initial value is

[TABLE]

This problem can now be solved by, for example series expansion, and the convergence of the series should be good.

10.5 Nonlinear stability theory

The main focus of this class is the theory of linear partial differential equations. The reason for this is two-fold. Firstly, many systems of scientific and technological interest can be described by such equation. We have seen several modeling examples where this turns out to be the case. Secondly, there is no theory of nonlinear PDEs with a generality anywhere close to the case for linear PDEs. Nonlinear PDEs in general require an individual treatment. In fact some people spend their whole working career studying one or a few such nonlinear equations. This is, however, if one is looking for exact solutions, precise statements about solution spaces etc. If one is satisfied with approximate answers, analytic or numeric, a lot more can be achieved.

The reason for this is that in many situations the nonlinear terms in the PDEs are much smaller than the linear ones. For example, we saw at the beginning of this class that the PDE

[TABLE]

is a simple model of traffic flow on a one-lane road. Assuming $n$ is small, thus the density of cars is low, we can Taylor expand $a(n)$ and truncate at first order

[TABLE]

Inserting (439) into (438) and rearranging the terms we get the nonlinear equation

[TABLE]

By assumption (439), the nonlinear term in this equation is smaller that the linear ones.

The use of truncated expansions like (439) are very common when modeling physical systems by PDEs. We saw this for example for the vibrating string.

The refractive index, $n$, is a material property that determines the optical response of dielectrics and gases. For inhomogenous materials, like corrugated glass for example, it will depend on position $n=n(\mathbf{x})$.

An important equation in the theoretical description of light is

[TABLE]

This is a equation which model light under normal conditions well. It is clearly a linear equation. For light under more extreme conditions, for example light generated by a laser, one find that the refractive index depends on the light intensity $I=|E|^{2}$. We thus in general have have $n=n(\mathbf{x},I)$. This functions has been measured for many materials and one find that in many cases a good approximation is

[TABLE]

where both $n_{0}$ and $n_{2}$ in general depends on position. Inserting (442) into (441) and keeping only terms that are linear in the intensity $I$, we get the equation

[TABLE]

This is now a nonlinear equation. The factor $n_{2}$ is for most materials very small, so ligh of very high intensity $I=|E|^{2}$ is required in order for the nonlinearity to have any effect. From your basic physics class you might recall that light tend to bend towards areas of higher refractive index; this property is the basis for Snell’s law in optics. Thus the intensity dependent index (442) will tend to focus light into areas of high light intensity. This increases the intensity, creating an even larger index, which focus light even more strongly. This is a runaway effect that quickly creates a local light intensity high enough to destroy the material the light is traveling through.

The study of nonlinear properties of high-intensity light is a large and very active field of research which has given us amazing technologies, like the laser, optical fibers, perfect lenses of microscopic size and even invisibility cloaks! (Well almost anyway)

For some physical systems there are no good linear approximations. One such area is fluid dynamics. The most important equation in this field of science is the Navier-Stokes equation.

[TABLE]

where $\rho_{0}$ is the mass density, $\mathbf{u}(\mathbf{x},t)$ the fluid velocity field, $p(\mathbf{x},t)$ the pressure in the fluid and $\mu$ the viscosity (friction). The second term on the left hand side of (444) is clearly nonlinear and in many important applications not small. The essential nonlinearity of the Navier-Stokes equation gives rise to the phenomenon of turbulence and other physical phenomena that are extremely hard to describe and predict. In fact, it is not even known, after almost 200 years, if the Cauchy problem for (444) is well posed. There is a possibility that solutions of (444) develop singularities in finite time. If you can solve this problem, the Clay mathematics institute will pay you one million dollars!

For systems where the nonlinear terms are small, one would expect that methods introduced for solving linear PDEs should be useful, somehow. They are, and we will illustrate this using a nonlinear heat equation

[TABLE]

The parameter $\hat{\lambda}$ measures the strength of the nonlinearity. The equation can model situations where we have a nonlinear heat source. The equation (445) clearly has a solution

[TABLE]

This corresponds to a situation where the system has a uniform temperature equal to zero. We know that if the nonlinear heat source is not present, any initial disturbance will decay, and the system will return to a uniform zero temperature state. We have seen this happening using normal modes and also Fourier series solutions. The question we want to ask is if $u_{0}(x,t)=0$ is stable also in the presence of the nonlinear heat source. We will investigate this problem both for an infinite system, $-\infty<x<\infty$, and a finite one $(0,l)$. We will see that the stability condition is different in the two cases. This is a general fact: Finite and infinite systems display different stability behavior.

Let us first set the problem up. Since we will be concerned with solutions close to $u_{0}(x,t)=0$, we introduce a small parameter $0<\epsilon<1$ and write the initial condition as

[TABLE]

where $h(x)$ is of normal size, $h(x)=\mathcal{O}(1)$. We are similarly introducing a $\mathcal{O}(1)$ function $w(x,t)$ through

[TABLE]

The equation (445) now gives

[TABLE]

The first step is to do a linear stability analysis. For this we recognize that, if $w=\mathcal{O}(1)$, the last term in (449) is much smaller than the other terms because $\epsilon<<1$. We therefore drop the small term and are left with the equation

[TABLE]

(450) is an arbitrary good approximation to (449) if $\epsilon$ is small enough and $w=\mathcal{O}(1)$. Going from (449) to (450) is called linearization. Linearization is an idea of outmost importance in applied mathematics.

We will now investigate the linearized equation (450) for the infinite and finite case separately.

10.5.1 The infinite case

In the infinite case $-\infty<x<\infty$, we can analyze the stability of equation (450) using normal modes as discussed previously in this class. Inserting a normal mode solution

[TABLE]

into (450), we find the algebraic equation

[TABLE]

From this we see that the stability index for the equation (450) is

[TABLE]

The conclusion is then that

1. i)

$u_{0}$ is stable if $\hat{\lambda}<0$, 2. ii)

$u_{0}$ is unstable if $\hat{\lambda}>0$.

In the unstable case the normal mode (451) experience unbounded exponential growth. However, when the growth makes

[TABLE]

the nonlinear term is of the same size as the linear one and the linearized equation is not valid anymore. From (455) this happend when

[TABLE]

which gives us, using (451), a break-down time

[TABLE]

So what happens after the break-down time? Observe that the nonlinear equation (445) has another simple uniform solution

[TABLE]

We can now ask when this solution is stable. We introduce a small parameter $\epsilon$ and a $\mathcal{O}(1)$ function $v(x,t)$ through

[TABLE]

Inserting this into (445) we find

[TABLE]

Linearizing, we find the equation

[TABLE]

Inserting a normal mode now gives

[TABLE]

and we can conclude

1. i)

$u_{1}$ is stable if $\hat{\lambda}>0$, 2. ii)

$u_{1}$ is unstable if $\hat{\lambda}<0$.

Thus, the solution $u_{1}$ is stable when $u_{0}$ is unstable. A reasonable guess now is that a small perturbation to the solution $u_{0}$ in the unstable regime $(\hat{\lambda}>0)$ grows and eventually approaches the stable solution $u_{1}$. Solving equation (445) numerically with a broad selection of small initial data $u(x,0)=\epsilon h(x)$ is needed to support our conjecture. Such situations can not prove, in the mathematical sense, that our conjecture is true. However, in real-life applied mathematics, proofs are few and far between, and we must use the tools at our disposal to do as much as we can to dispel uncertainty as to the validity of our conclusions.

Another way to support our conjecture is to find solutions for special initial data. Let us consider the special initial data

[TABLE]

Let us look for a solution to (445) that is independent of $x$, $u=u(t)$. Such a solution must satisfy

[TABLE]

This first order ODE that can be solved exactly

[TABLE]

It is easy to verify that

[TABLE]

for the case $\hat{\lambda}>0$. Thus this special solution lends further support to our conjecture.

10.5.2 The finite case

Let us consider the nonlinear heat equation (445) on a finite interval $0<x<\pi$ with zero boundary conditions

[TABLE]

Recall that linearizing (442) around the simple solution $u_{0}(x,t)=0$ gives us an equation of the form

[TABLE]

We can solve this problem using the finite Fourier transform. Introduce the operator $L=-\partial_{xx}$. The relevant eigenfunctions satisfy

[TABLE]

We have analyzed this eigenvalue problem before. The eigenvalues are

[TABLE]

and the normalized eigenfunctions are

[TABLE]

The associated finite Fourier transform of a given function $f(x)$ is

[TABLE]

and the inverse finite Fourier transform is

[TABLE]

We now express an arbitrary solution to (467) in terms of an inverse finite Fourier transform

[TABLE]

The evolution equations for the $N_{k}(t)$ are found from (467) if we multiply by $M_{k}(x)$ and integrate

[TABLE]

For each of the three terms in this equation we have

[TABLE]

[TABLE]

[TABLE]

Inserting (466), (476) and (478) into (474) gives the following uncoupled system of ODEs for the unknown functions $N_{k}(t)$

[TABLE]

The general solution to (467) can thus be written as

[TABLE]

where the normal modes, $w_{k}(x,t)$, are

[TABLE]

We see that

1. i)

$\hat{\lambda}<1\,\Rightarrow$ all normal modes decay $\Rightarrow\,u_{0}(x,t)=0$ is stable. 2. ii)

$\hat{\lambda}>1\,\Rightarrow$ at least one normal mode grows exponentially $\Rightarrow\,u_{0}(x,t)$ is unstable.

For the infinite domain case the stability condition was $\hat{\lambda}<0$. Thus the stability condition is not the same for the case of infinite and finite domain. This is in general true.

For the unstable case, $\hat{\lambda}>1$, one or more normal modes experience exponential growth, which after a finite time invalidate the linearized equation (467). What is the system going to do after this time? This is a questing of nonlinear stability.

Nonlinear stability is a much harder problem than linear stability but powerful techniques exist in the form of asymptotic expansions which in many important cases can provide answer to the question of nonlinear stability.

Here we will not dive into this more sophisticated approach but will rather introduce an approach based on the finite Fourier transform.

Observe that the finite Fourier transform can also be applied to the nonlinear equation

[TABLE]

Representing the solution using the Finite Fourier transform we have

[TABLE]

The evolution equations for the $N_{k}(t)$s are found as for the linearized case. We have

[TABLE]

The only new term as compared with the linearized case is the second term on the right hand side of (484). For this term we have

[TABLE]

where

[TABLE]

Using (475), (476), (478) and (485) in (484) gives us

[TABLE]

(487) is an infinite system of nonlinear coupled ODEs. This is evidently a very hard problem to solve in any exact sense. (487) shows the complexity of nonlinear PDEs. The very compact notation from calculus makes the nonlinear PDE (482) to appear simple, but as (487) shows the problem of solving (482) is actually horrendously difficult.

In order to get anywhere, we must restrict to situations where some of the complexity of (487) is reduced. Finding such situations is what good applied mathematicians are trained to do.

Let us assume that $\hat{\lambda}>1$ so we are in the unstable regime, but that $\hat{\lambda}<4$ so that only the normal mode $u_{1}(x,t)$ is unstable. All the other ones decay exponentially. For the $N_{k}(t)$s this means that

[TABLE]

Thus, on the right hand side of (487) the term

[TABLE]

is much larger than the other terms. A good approximation to (487) therefore appears to be

[TABLE]

The coefficient $a^{111}_{1}$ can be calculated using (486)

[TABLE]

[TABLE]

If we multiply this equation by $N_{1}(t)$ we obtain

[TABLE]

where $y(t)=N_{1}^{2}(t)$. (493) is a so-called Ricatti equation and can be solved exactly. There is a general way to solve such equations.

A general Ricatti equation is of the form

[TABLE]

First introduce a function $v(x)$ by

[TABLE]

Then we get the following equation for $v(x)$

[TABLE]

where

[TABLE]

Next, introduce the function $u(x)$ by

[TABLE]

Then we get the following linear equation for $u$

[TABLE]

Using (498) and (495), any solution of (499) gives a solution of the Ricatti equation (494). Since the general solution to (499) can be found explicitly in some cases, the general solution (494) can also be found for those cases.

For our equation (493) we have

[TABLE]

We can find the general solution to this equation which leads to the general solution (493). With initial condition given as in (447) we get

[TABLE]

This solution tends to a finite stationary value as $t\to\infty$

[TABLE]

rather than growing unbounded as the linearized equation predicts.

Note that the current method can be extended to study more general nonlinear terms $\hat{\lambda}f(u)$.

11 Integral transforms

Integral transforms are used to solve PDEs on unbounded domains. They have very wide applicability and play a fundamental role in applied mathematics and theoretical physics.

In these notes I will concentrate on equations of the form

[TABLE]

where

[TABLE]

and where

[TABLE]

Recall that the basic idea in the separation of variables is to look for special solution to (503) of the form

[TABLE]

($t\to y$ in the restricted elliptic case). We found that a separation constant $\lambda$ appeared and that (506) is a solution to (503) only if

[TABLE]

and

[TABLE]

Here we have written the separation constant as $\lambda^{2}$ for convenience. Equation (507) must be supplied with boundary conditions. They must of course be inherited from the boundary conditions on equation (503). The most common ones are that the solution is bounded or that it vanishes as $\|x\|\to\infty$. Thus we will require that $M(\mathbf{x})$ is bounded or vanishes at infinity.

As for the case of bounded domains, (507) with boundary conditions will not have solutions for all $\lambda$. The $\lambda$s that give a solution are part of the spectrum of $L$. Things look very similar to the case of bounded domains, so far. There is however a big difference as compared with the bounded case.

In order to illustrate this difference, let us consider the boundary value problem

[TABLE]

The general solution to the equation is for $\lambda>0$

[TABLE]

The boundary conditions imply that

[TABLE]

[TABLE]

In order for (512) to have a nontrivial solution, and thus for $\lambda$ to be in spectrum of $L=-\partial_{xx}$, the matrix in (512) must have zero determinant

[TABLE]

There are no eigenvalues for $\lambda\leq 0$. Thus (514) defines the spectrum of $l$.

Observe that the distance $\Delta\lambda$ between two spectral points is

[TABLE]

When the size of the domain increases without bound we observe that the distance between spectral points go to zero

[TABLE]

Thus in this limit the spectrum of $L$ will form a continuum

[TABLE]

This is in general the case for unbounded domains; the spectrum of the operators form a continuum.

For the discrete case, with eigenvalues $\{\lambda_{k}\}_{k=1}^{\infty}$ and corresponding eigenfunctions $\{M_{k}\}_{k=1}^{\infty}$, the basic idea was to look for solutions to the PDE in the form of an infinite sum

[TABLE]

For the unbounded case we have eigenfunctions $M_{\lambda}$ and corresponding $N_{\lambda}$, where the $\lambda$s form a continuum that we denote by $D$. By analogy with (512) we seek solutions to the PDE in the form of "sums over the continuum". Such sums are of course integrals. Thus

[TABLE]

We recall that in the case of finite domains, the fitting of initial conditions to our solution consisted in finding constants $\{a_{k}\}_{k=1}^{\infty}$ such that

[TABLE]

Using orthonormality of the eigenfunctions we found that

[TABLE]

Fitting initial conditions for the case of unbounded domains will by analogy involve finding coefficients $a_{\lambda}$, $\lambda\in D$ such that

[TABLE]

In order to determine $a_{\lambda}$ we need a continuum analog of the orthonormality condition from the discrete case

[TABLE]

where we recall that $\delta_{kk^{\prime}}$ is the Kronecker delta.

If we had something like a continuum version of the Kronecker delta one would guess at

[TABLE]

This is in fact what Paul Dirac did early in the previous century. Dirac was an English genius that was one of the founding fathers of quantum theory. He came to this idea while trying to solve the Schrodinger equation from quantum mechanics. He found that $\delta_{\lambda\lambda^{\prime\prime}}$ only depended on $\lambda-\lambda^{\prime}$ and he wrote it as

[TABLE]

$\delta(\lambda)$ is called the Dirac delta function after him. In order to do its intended job he found that it would have to have the following two properties

[TABLE]

[TABLE]

It is a trivial exercise to prove that no function exists that have the two properties (526) and 527). Mathematicians were quick to point this out but Dirac didn’t care, he did not know much mathematics, and he did not have to know it either, because he was in fact a genius and people like him move to a different drum.

Eventually mathematicians did make sense of the Dirac delta function and in fact discovered a whole new mathematical universe called the theory of distributions. We will discuss elements of this theory next semester. In this semester we will use objects like the Dirac delta sparingly and when we do only in a purely formal way like Dirac himself.

There is an important moral to the story of the Dirac delta. Pure mathematics is important because it gives us abstract frameworks to help organize our mathematical methods. Pure mathematics can also give us more efficient ways of doing our calculations and make precise under which conditions our manipulations will keep us within the boundaries of the known mathematical universe.

However the framework and bounds imposed by pure mathematics should not stop us from sometimes disregarding the rules and play with the formulas, like Dirac did. In the final analysis, virtually all of mathematics originate from such rule-defying play.

It is however also appropriate to give a warning; disregarding the rules and playing with the formulas should not be a cover for sloppy thinking. This is unfortunately the case quite frequently in applied science.

We can’t all be a Dirac. Disregarding the rules imposed by pure mathematics will, by most people, most of the time, produce useless nonsense.

In these notes I will sometimes break the rules a la Dirac and I will try to remember telling you when I do. All my rule breaking can be made sense of within the universe of distributions.

All finite Fourier transforms will in the limit of unbounded domains give rise to integral transforms. There is thus a large set of such transforms available, each tailored to a particular class of problems. In these notes I will focus on the Sine-, Cosine-, Fourier- and Laplace transforms.

11.1 The 1D Fourier transform

Let us consider the eigenvalue problem

[TABLE]

The boundary condition is that $M(x)$ is bounded when $|x|\to\infty$. The general solution of the equation (528) is

[TABLE]

where $\lambda\in\mathbb{C}$. Writing $\lambda=\alpha+i\beta$ we have

[TABLE]

and the functions $e^{\mp\beta x}$ are unbounded on $-\infty<x<\infty$ as long as $\beta\neq 0$. Thus to satisfy the boundary conditions we must restrict $\lambda$ to the real axis $-\infty<x<\infty$. On the real axis impose the conditions

[TABLE]

then the boundedness of $M_{\lambda}$ follows from the triangle inequality

[TABLE]

Therefore the spectrum in this case consists of the whole real axis. We can now use these eigenfunctions to represent functions $f(x)$ in the following way

[TABLE]

where $c(\lambda)=a(-\lambda)+b(\lambda)$. Note that in (532), $f(x)$ is represented by $e^{i\lambda x}$ in terms of a standard improper integral. In the theory of integral transforms, other, less standard, improper integrals appear. We will discuss these when they arise.

The formula (532) is by definition the inverse Fourier transform. The key question is how to find the $c(\lambda)$ corresponding to a given function $f(x)$.

For this purpose, let us introduce a small number $\epsilon>0$ and define

[TABLE]

Observe that for all real $\lambda$ we have

[TABLE]

Thus

[TABLE]

This result makes the following formula plausible (I am playing with formulas now)

[TABLE]

In the theory of distributions it becomes clear that the constant of proportionality is $2\pi$. So we have the formula

[TABLE]

Since the Dirac delta function is zero for $\lambda\neq 0$ it appears that for all reasonable functions, $f(\lambda)$, we have

[TABLE]

(I am playing now.)

Let us now return to our original question of how to find the $c(\lambda)$ corresponding to a given function $f(x)$ in formula (532)

[TABLE]

Multiplying (542) by $e^{i\hat{\lambda}x}$ and integrate with respect to $x$ over $-\infty<x<\infty$ gives

[TABLE]

So the formula we are looking for is

[TABLE]

This formula defines the Fourier transform. Note that there are several conventions about where to place the factor of $2\pi$. A more symmetric choice[1], is

[TABLE]

[TABLE]

Note that not only are there different conventions on where to place the factor of $2\pi$ but also in which exponent to place the minus sign.

In this class (546) defines the Fourier transform and (547) defines the inverse Fourier transform.

Pure mathematicians have made extensive investigations into for which pair of functions $f(x)$ and $F(\lambda)$ the formulas (546) and (547) are valid. For example if $f(x)$ is piecewise continuously differentiable on each finite interval and if additionally

[TABLE]

formula (546) and (547) holds pointwise. This result and related ones are important but we will use (546) and (547) in a much wider context than this.

Let $f(x)=\delta(x)$. Then (546) gives

[TABLE]

and (547) then gives

[TABLE]

where we have used the formula (540) which we argued for previously. More generally, the pair of formulas (546) and (547) holds for any distribution, not just the Delta distribution. More about this next semester. In this semester we will basically use (546) and (547) to take Fourier transform and inverse Fourier transform of anything we want!

The Fourier transform has many powerful properties. You can find them in handbooks or on the web whenever you need them. For example if $F(\lambda)$ and $G(\lambda)$ are Fourier transforms of $f(x)$ and $g(x)$ then $aF(\lambda)+bG(\lambda)$ is the Fourier transform of $af(x)+bg(x)$. A less trivial, and extremely useful for us, property involves products of Fourier transforms $F(\lambda)G(\lambda)$ of functions $f(x)$ and $g(x)$. Such products often occur when we use Fourier transforms to solve ODEs and PDEs.

Playing with formulas we get the following

[TABLE]

The integral on the right hand side of (11.1) is known as a convolution integral. Canceling the factor $\frac{1}{\sqrt{2\pi}}$ we get the convolution integral theorem for Fourier transforms

[TABLE]

Evaluating at $x=0$ we get

[TABLE]

Now observe that

[TABLE]

where $\overline{F}(\lambda)$ is the complex conjugate of $F(\lambda)$.

Thus the fourier transform of $\overline{f}(-t)$ is $\overline{F}(\lambda)$. Putting $g(t)=\overline{f}(-t)$ into (553) we get

[TABLE]

The identity

[TABLE]

is called the Parseval equation and is clearly a generalization of the Parseval equation for generalized Fourier series to continuous spectrum.

Another important property of Fourier transforms involve derivatives. Observe that

[TABLE]

Thus, derivatives are converted into multiplication by the Fourier transform. (11.1) can be easily generalized to

[TABLE]

Formulas (11.1) and (558) are the main sources for the great utility of Fourier transforms in the theory of PDEs and ODEs.

11.1.1 A boundary value problem for an ODE

Let us consider the following boundary value problem

[TABLE]

where $k$ is a fixed constant and $f$ a given function. The boundary conditions attached to (559) are

[TABLE]

Let $F(\lambda)$ and $Y(\lambda)$ be the Fourier transforms of $f(x)$ and $y(x)$. Then the transform of the equation (559) is, upon using (558), given by

[TABLE]

This is a simple algebraic equation whose solution is

[TABLE]

Using the inverse Fourier transform and the convolution theorem we get the following representation of $y(x)$:

[TABLE]

where

[TABLE]

Using a table of transforms or residue calculus from complex function theory we have

[TABLE]

so that the solution $y(x)$ is

[TABLE]

Because of the exponentially decaying factor, $y(x)$ and $y^{\prime}(x)$ will both go to zero as $|x|\to\infty$ if $f(t)$ decays fast enough at infinity.

Actually, computing $y(x)$ for a given $f(t)$ using formula (566) is no simple matter; tables of integrals, complex function theory, numerical methods and asymptotic methods are all of help here.

For some simple choices for $f(t)$ the integral in (566) can be evaluated easily. Let us first choose $f(t)=1$. Formula (566) implies

[TABLE]

Observe that

[TABLE]

so $y(x)$ defined in (567) is a solution to the equation (559), but it does clearly not satisfy the boundary conditions (560). The reason why this occurs is that in deriving the formula (566) we assumed that $f(x)$ had a Fourier transform $F(\lambda)$. However $f(x)=1$ leads to the divergent integral

[TABLE]

for $F(\lambda)$. Using (550) we conclude that $F(\lambda)$ is not a function but a distribution

[TABLE]

Calculations like the one leading up to the "solution" (567) are very common whenever transforms are used to solve PDEs. They are ideal calculations that is meant to represent the limit of a sequence of problems, each producing an honest to God solution to the boundary value problem (559),(560). Such ideal calculations are often possible to do in closed form and does often lead to some insight into the general case, but in what sense the result of the calculation is a limit of a sequence of actual solution is for the most part left unsaid. This will be clarified when we discuss distributions next semester. However, let us preview next semester by saying exactly what we mean by a limit here.

Define a sequence of functions $\{f_{N}(x)\}$

[TABLE]

This function has a Fourier transform for each $N$ and our formula (566) should give a solution to the boundary value problem (559),(560)

[TABLE]

Let us first assume that $|x|\leq N$

[TABLE]

The calculations for $|x|>N$ is completed in a similar way. This gives us

[TABLE]

$y_{N}(x)$ is a solution to (559) and it is easy to verify that $y_{N}(x)$, $y^{\prime}_{N}(x)\to 0$ when $|x|\to\infty$, so we have an infinite sequence of solutions $\{y_{N}(x)\}$ to the boundary value problem defined by (559),(560).

For any $x$ there evidently exists a $N$ so large that $N\geq|x|$. Therefore for such an $x$ we have $f_{N}(x)=f(x)$ for all $N>|x|$ and consequently the sequence of functions $\{f_{N}\}_{N}^{\infty}$ converge pointwise to the functions $f(x)=1$. This convergence is however not uniform.

Observe that for any $x$ there obviously exists inifinitely many $N$ such that $N\geq|x|$. For such $N$ we have from (575) that

[TABLE]

when $N\to\infty$

Thus, this convergence, like the one for the source $f(x)$, is also pointwise but not uniform. This non-uniformity of the limit is why the limiting function, $y(x)$, does not satisfy the boundary condition (560) even though all $y_{N}(x)$ does.

As a second example of an ideal calculation let us choose $f(t)=\delta(t)$. Then formula (566) gives

[TABLE]

This function satisfies the boundary condition and the equation for all $x\neq 0$. It is however not differentiable at $x=0$.

Introduce a sequence of functions $\{g_{N}(x)\}$ by

[TABLE]

Using formula (566) with $f(x)=g_{N}(x)$, it is straight forward to show that the corresponding solution is

[TABLE]

For any $x$ there exist infinitely many $N$ such that $|x|\geq\frac{1}{N}$, and for such $N$ we have from (579) that

[TABLE]

when $N$ approach infinity. Furthermore, for $x=0$ we have from (579) that

[TABLE]

when $N$ approach infinity.

We have thus established that the infinite sequence of solutions (579) of equation (559) corresponding to the sequence $\{g_{N}\}$ converge pointwise to the solution $y(x)=\frac{1}{2k}e^{-k|x|}$ corresponding to the Dirac delta function $f(x)=\delta(x)$. The sequence (579) does not converge uniformly to the function $y(x)$ and this explains the fact that all functions from the sequence (579) are differentiable at $x=0$ whereas the limiting function $y(x)=\frac{1}{2k}e^{-|x|}$ is not.

In the previous example the sequence $\{f_{N}(x)\}$ converged pointwise to the limiting function $f(x)=1$ that generated the solution $y(x)=\frac{1}{k^{2}}$. Is it similarly true that in the current example the sequence $\{g_{N}(x)\}$ converge to the Dirac delta function that generated the solution $f(x)=\frac{1}{2k}e^{-k|x|}$?

In order to answer this question we start by observing that for all $N$

[TABLE]

so $g_{N}(x)$ satisfies the same integral identity (527) as the Dirac delta function (Here $D=(-\infty,\infty)$). Furthermore for any $x\neq 0$ there exists a $N_{x}$ such that

[TABLE]

Therefore for all $N\geq N_{x}$, we have $g_{N}(x)=0$. Therefore

[TABLE]

so the limit of the sequence $g_{N}$ when $N\to 0$ also satisfies the identity (526). So in fact

[TABLE]

I am playing with formulas here; the limit in (583) does not mean the usual thing. After all, $\delta(x)$ is not a function but something different, namely a distribution. It would perhaps be more fair to say that the limit of the sequence $\{g_{N}(x)\}$, when $N\to\infty$ does not exist, since there can be no function satisfying (526) and (527).

What you see here is an example of a very subtle and profound mathematical idea called completion. You have actually used this idea constantly since grade school, but you may not have been told.

When we learn children arithmetic we start with whole numbers and eventually move on to fractions or rational numbers as they are called. The universe of numbers are rational numbers at grade school level. At some point the children are confronted with the equation

[TABLE]

If the teacher is any good he will by trial and error convince the children that (584) has no solutions. So what is there to do? He will say that (584) actually has a solution that we denote by $\sqrt{2}$, and he will show, perhaps using a calculator, how to find a sequence of rational numbers that more and more closely approximate $\sqrt{2}$. In a similar way the teacher introduces $\sqrt{3}$, $\sqrt{5}$ etc. and he develops rules for calculations with such things.

But the children has been tricked. When the teacher argues that calculated sequences of rational numbers approximate $\sqrt{2}$ better and better and in the limit equals $\sqrt{2}$, he has been pulling wool over their eyes. What he has done is to complete the set of rational numbers, in the process creating the real numbers. There is no sequence of rational numbers converging to $\sqrt{2}$, in reality $\sqrt{2}$ is defined to be a certain set of sequences of rational numbers. The existence of $\sqrt{2}$ is postulated, pure and simple, through the completion process. After the completion process has lead to the set of real numbers, pure and applied mathematics can give them additional intuitive meaning by mapping them onto other physical and mathematical systems. In this way real numbers have come to be interpreted as geometrical points on a line etc.

The relation between functions and distributions is exactly the same as the one between rational and real numbers. For the case of distributions Dirac is the teacher and the rest of us the children. Next semester we will map distributions onto a system derived from linear algebra and in this way gain more insight into the nature of distributions.

The particular "solution" (577) found by choosing $f(x)=\delta(x)$ is an example of what is called a Green’s function. Writing

[TABLE]

We have seen that $g(x-y)$ is a solution to the equation

[TABLE]

where we are taking the derivative with respect to the variable $x$, the variable $y$ is a parameter here. Green’s functions are of outmost importance in applied mathematics and theoretical physics and will be discussed extensively next semester. As you can see from (586), Green’s functions are solutions to differential equations formulated in the universe of distributions. For this reason Green’s functions are not actually functions but distributions.

Let us now proceed to consider a slightly different ODE

[TABLE]

For this equation the boundary conditions can not be that $y(x)$ and $y^{\prime}(x)\to 0$ as $|x|\to\infty$. The reason is that (587) does not have solutions decaying at infinity. In fact for $f(x)=0$ the solution consists of $\sin kx$ and $\cos kx$, which are bounded but does not vanish. We will therefore only require that $y$ and $y^{\prime}$ are bounded at infinity. This implies, by an argument identical to the one on page 11, that the spectrum is real. We apply the Fourier transform like for the previous example and find

[TABLE]

We observe that, compared to the previous case, a problem has appeared: $Y(\lambda)$ is not defined at $\lambda=\pm k$. The proper way to handle this is to use complex function theory and for more complicated situations complex function theory becomes a necessity. Here we will proceed in a more heuristic way by observing that (587) becomes (559) if we let $k\to ik$ or $k\to-ik$. Making the same replacement in the solution formula (566) we get two solutions

[TABLE]

Both of $y_{\pm}(x)$ satisfy (587) and are bounded if $f(t)$ is absolutely integrable $\left(\int^{\infty}_{-\infty}\mathrm{d}t|f(t)|<\infty\right)$. Thus the boundary value problem (587) does not have a unique solution. Why is that?

Let us start with the wave equation

[TABLE]

and let the source, $f(x)$, be supported in some bounded domain $G$. Outside this domain $G$, (590) is the homogeneous wave equation whose general solution consists of left- and right going waves.

We assume that the source is the originator of all waves. Thus, no source, no waves. Then it is evident that there can be no left moving waves to the right of $G$ and no right-moving waves to the left of $G$. This is a physical condition satisfied by all known wave phenomena; waves move away from the source that creates them. This is in fact part of the definition of what a physical source is!

Let us separate variables in (590)

[TABLE]

Defining $k=\frac{\omega}{c}$ we find our equation (587). The homogeneous equation (592) has a solution proportional to $e^{\pm ikx}$. Taking these solutions together with (591) we observe that

[TABLE]

The physical boundary condition thus consists of killing terms containing $e^{-ikx}$ to the right of the support of $f$ and terms containing $e^{ikx}$ to the left of the support of $f$. This can be expressed by saying that

[TABLE]

where $\text{Supp(f)}=\{x;f(x)\neq 0\}$. Note that Supp(f) by definition includes its boundary points and is thus a closed set.

We have here assumed that the support of $f$ is bounded. However it is enough that $f(x)\to 0$, $|x|\to\infty$. For this more general case we pose the conditions

[TABLE]

(594) or more generally (595) is called the radiation condition at infinity, or outgoing waves at infinity and is the correct boundary condition for equations describing wavelike phenomena. (595) picks out the unique solution $y_{+}(x)$ to our boundary value problem (587),(595). Choosing $f(x)=\delta(x)$ gives us the following Green’s function for this case

[TABLE]

11.1.2 The Cauchy problem for the heat equation

Let us next apply the Fourier transform to the Cauchy problem for the heat equation

[TABLE]

In order to apply the Fourier transform we must assume that $f$ has a Fourier transform and that $u$, $u_{x}$ goes to zero at $x=\pm\infty$.

Taking the Fourier transform with respect to $x$ of the heat equation (597) we get

[TABLE]

where

[TABLE]

The equations (598) are uncoupled for different values of $\lambda$. It is the linearity of the heat equation that produce an uncoupled system of ODEs for the set of functions $U(\lambda,t)$ indexed by the spectral parameter $\lambda$. Note how similar this is to what we got when we applied the finite Fourier transform to linear PDEs on a finite domain.

The initial value problems (598) are easy to solve. We get

[TABLE]

Using the inverse Fourier transform we get the solution

[TABLE]

where we have changed order of integration. Defining

[TABLE]

the solution can be written in the form

[TABLE]

Note that we could have gotten this expression more directly by using the convolution theorem. Observe that

[TABLE]

This integral, because of the exponential factor present for $t>0$, converges very fast, so we can differentiate under the integral sign

[TABLE]

so

[TABLE]

This ODE is simple to solve. The solution is

[TABLE]

and

[TABLE]

(This is a Gaussian integral.) Thus

[TABLE]

The formula for $u(x,t)$ thus becomes

[TABLE]

If we formally use the initial condition

[TABLE]

we get the solution

[TABLE]

$G$ is yet another example of a Green’s function. In this context it is called the fundamental solution to the heat equation. The fundamental solution is the temperature profile generated by a point-like heat source located at $x=x^{\prime}$ for $t=0$.

Note that it is an actual classical solution (verify!) for all $t>0$ even if the initial condition is extremely singular. This is yet another example illustrating the extreme smoothening properties of the heat equation.

From the solution formula (611) it is evident that even if the initial heat profile is zero beyond a finite interval $[-L,L]$, the solution is nonzero for any point $x\in\mathbb{R}$ for any $t>0$, however small $t$ is.

The heat equation clearly does not respect causality.

Also note that for the fundamental solution we have

[TABLE]

Thus

[TABLE]

The fundamental solution is thus a sequence of smooth functions converging to the Dirac delta distribution. Which is exactly what it needs to do in order to satisfy the initial condition!

We can argue for the correctness of (614) as follows. Observe that as $t\to 0^{+}$ the function $G(x-\xi,t)$ becomes more and more localized around $x=\xi$. Thus for $t$ close enough to zero and for all reasonable functions $g(x)$ we have

[TABLE]

because

[TABLE]

Verify this!

11.1.3 The Cauchy problem for the wave equation

Let us consider the problem

[TABLE]

Taking the Fourier transform of this problem, proceeding like for the heat equation, we get

[TABLE]

The initial value problem (619) is entirely standard and the solution is

[TABLE]

Taking the inverse transform we get

[TABLE]

By the definition of the inverse Fourier transform we have

[TABLE]

But then

[TABLE]

where D is a constant. Given (623), (621) can be written as

[TABLE]

This is the d’Alembert formula, now derived using the Fourier transform.

There is however something rather fishy with what we did. Let us say that we did not realize that the identity (623) holds. Then, for given $F(\lambda)$, $G(\lambda)$, we would have to find $u(x,t)$ by solving the integrals in (621) analytically or numerically. We should however worry about the last two integrals. They are of the general form

[TABLE]

where in general $H(0)\neq 0$. The integral is singular at $\lambda=0$ and the question how the integral should be understood must be raised.

Away from $\lambda=0$ there are no problems, so the question we need to ask is if we can make sense of integrals of the type

[TABLE]

where $L$ is arbitrary but fixed. In calculus we have seen that at least some integrals with a singular integrand in the domain of integration can be made sense of by considering it to be an improper integral. Let us try to use this approach on (626):

[TABLE]

This limit does not exist. Thus (626) can not be given meaning as an improper integral. Observe however that if we in (627) let $\epsilon=\delta$ and take the difference before we take the limit we get zero, which is a finite answer. Thus what we are proposing is to understand (626) in the following way

[TABLE]

When we use the definition (628) we understand the integral as a Cauchy principal value integral. This is now a new way for you to calculate integrals with a singular integrand. In general for a function $f(x)$, with a singularity at $x=x_{0}$, the definition is

[TABLE]

where $PV$ signifies that we are calculating the principal value.

The principal value is finite only if the singularity in $f$ is of the right type. For example for $f(x)=\frac{1}{x^{2}}$ we get

[TABLE]

If we interpret the two integrals containing $\dfrac{1}{\lambda}$ in (615) as Cauchy principal value integrals, we can use the formula to calculate the solution to the wave equation.

However, there is still a weak fishy smell here; why is interpreting the singular integrals in (615) as Cauchy principal value integrals the right thing to do? How do we know that this will give the same solution as, say, the d’Alembert formula? Could one not imagine some other funny way of getting a finite number out of integrals like (626)? And would this not perhaps produce a different solution; which should not be possible because the Cauchy problem for the wave equation has a unique solution? What is going on?

Answering this question, which is an important one, would take us too far from the main focus of the course, so I will refrain from giving an answer. Here I will only ensure you that this is the right way to understand the integral (626) and it does give the same solution as the d’Alembert formula. It is also worth noting that questions like the one we have been considering frequently turn up when we use transforms to solve PDEs, and solving them can often be a deep and vexing problem.

11.1.4 The Laplace equation in a half-plane

Let us consider the Laplace equation in the upper halfplane:

[TABLE]

We introduce the Fourier transform in $x$

[TABLE]

Using this Fourier transform on (631) gives us the following infinite system of uncoupled ODEs:

[TABLE]

This problem is standard and the solution is

[TABLE]

The inverse transform gives us the solution

[TABLE]

Furthermore, observe that

[TABLE]

and therefore we have the following expression for the solution of the Laplace equation in a half-space

[TABLE]

By direct substitution it can be verified that (637) is a solution to the Laplace equation for $y>0$.

When $y\to 0$ the factor $\dfrac{y}{\pi}$ goes to zero, but the integrand approaches

[TABLE]

which is singular at $x^{\prime}=x$. A careful study of the limit shows that the limit of (637) when $y\to 0$ is actually $f(x)$. We will not do this here.

For simple choices of $f(x)$ we can solve the integral (637). For example if $f(x)=H(x)$ where $H(x)$ is the Heaviside step function

[TABLE]

we find

[TABLE]

Substitute $u=\dfrac{x^{\prime}-x}{y}$ and $\mathrm{d}u=\dfrac{1}{y}\mathrm{d}x^{\prime}$ in this integral. This gives us the following explicit expression for $u(x,y)$

[TABLE]

This function satisfies the Laplace equation and

[TABLE]

so $u(x,y)$ satisfies the boundary condition at $y=0$. Also

[TABLE]

so $u(x,y)$ is bounded at $y=+\infty$.

However, $u(x,y)$ does not go to zero as $|x|\to\infty$ and so does not actually have a Fourier transform, which we assumed. The reason for this is that our boundary function $f(x)=H(x)$ does not have a Fourier transform in the conventional sense. Observe that $u(x,y)$ is infinitely smooth for $y>0$ even if it is discontinuous in $x$ for $y=0$. We are rediscovering the fact that the Laplace equation can not support solutions with singularities. As we have seen, this is a general feature of elliptic equations.

Let us say that we were rather to solve the Laplace equation with Neumann conditions

[TABLE]

We can reduce this to a problem of the previous type using a trick, invented by an English mathematician, called Stoke’s rule.

Introduce $v(x,y)=u_{y}(x,y)$. Differentiating (637) with respect to $y$ gives

[TABLE]

The solution to this problem is given by (637) and since

[TABLE]

we get

[TABLE]

Observe that $u(x,y)$ defined by (11.1.4) is not the only solution to (644).

For any solution $v(x,y)$ to (645) we can add an arbitrary function $C(x)$ to the right hand side of (646) and in this way define a whole infinite familty of functions $u(x,y)$

[TABLE]

All such functions satisfy the boundary condition $u_{y}(x,0)=g(x)$, but in order to satisfy the Laplace equation (644) the added function $C(x)$ must be linear. Furthermore, in order for the solution to be bounded in $x$ this linear function must be a constant.

11.2 The Fourier sine- and cosine transform

Let us consider the following eigenvalue problem

[TABLE]

on a semi-infinite 1D domain. We will consider two sets of boundary conditions

[TABLE]

Let us first consider (649) with boundary conditions (650) .

We will only consider (649) for real $\lambda$. For $\lambda=0$ the general solution of (649) is

[TABLE]

and

[TABLE]

Thus $\lambda=0$ is not an eigenvalue. For $\lambda\neq 0$ the general solution to (649) is

[TABLE]

and

[TABLE]

This function is clearly bounded for all real $\lambda$ and therefore all such $\lambda$ are in the spectrum. Observe however that

[TABLE]

Therefore, the eigenfunctions corresponding to $\lambda$ and $-\lambda$ are linearly dependent. We therefore loose nothing by restricting to the range

[TABLE]

and eigenfunctions

[TABLE]

For (649) with the boundary conditions (651) we find in a similar way the spectrum

[TABLE]

and eigenfunctions

[TABLE]

In order to find the transform corresponding to (655) and (656) it is useful to consider the finite domain

[TABLE]

and then let $l$ approach infinity. On the bounded domain (659), the boundary condition at $x=+\infty$ is replaced by the condition $M(l)=0$. This gives us the following normalized eigenfunctions

[TABLE]

where

[TABLE]

A function $\tilde{f}(x)$ can now be expanded in a Fourier series with respect to (660) according to

[TABLE]

[TABLE]

Observe that the distance between the points $\lambda_{k}$ on the $\lambda$-axis is

[TABLE]

We use this and rewrite (662), (663) in the following way

[TABLE]

[TABLE]

where $\tilde{F}_{k}=\sqrt{l}F_{k}$. If we formally let $l\to\infty$ in (665) and (666) we get

[TABLE]

[TABLE]

In order to get a more symmetric and pleasing form of the transform (667), (668) we introduce scalings

[TABLE]

[TABLE]

[TABLE]

Choose $\alpha=\sqrt{\pi}$, $\beta=1$. Then we get the Fourier sine transform

[TABLE]

Our friendly analysis can tell us for which pair of nice functions the improper integrals converge and for which (672) and (673) holds.

However we will apply it to much more general things and will proceed formally and hope that in the end everything will be ok in the universe of distributions. Here it is only worth noting that even though $\sin\lambda x$ is zero for $x=0$, $f(x)$ in (672) might not approach zero as $x$ approach zero. This is because in general the convergence of the improper integral will not allow us to interchange integral and limit.

A similar derivation for the boundary condition (651) gives us the Fourier cosine transform

[TABLE]

The Fourier sine- and cosine transforms have several interesting and useful properties[1]. Here we will discuss how they behave with respect to derivatives. Observe that

[TABLE]

If we assume that $f(x)$ vanishes at $x=+\infty$ we get the identity

[TABLE]

In a similar way we find, assuming that also $f^{\prime}(x)$ vanishes at infinity,

[TABLE]

Thus the Fourier sine transform of $f^{\prime}$ is expressed in terms of the Fourier cosine transform. The Fourier sine transform of $f^{\prime\prime}$ is expressed in terms of the Fourier sine transform of $f$ and a term depending on $f(0)$.

For the Fourier cosine transform we get in a similar way

[TABLE]

[TABLE]

We observe that in order to apply the Fourier sine- and cosine transforms to solve differential equations on a semi infinite domain the boundary conditions must be of a particular type. For second order equations the Fourier sine transform requires Dirichlet condition and the cosine transform requires Neumann condition. Since one typically solve time dependent problems on a semi-infinite domain $t>0$, one might think that the Fourier Sine- and Cosine transforms should be useful for this, but they are in fact not.

If the equation contains a first derivative with respect to time we will find that the Fourier Sine transform of $u_{t}$ is expressed in terms of the Fourier Cosine transform for $u$ and vice versa. If the equation contains two time derivatives, $u_{tt}$, there are two data given at $t=0$, $u$ and $u_{t}$. But both the Fourier Sine and the Fourier Cosine transform can only include one condition at $t=0$, $u$ for the Sine transform and $u_{t}$ for the Cosine transform.

We will see that for initial value problems, the Laplace transform is the appropriate one.

Another useful property of the Fourier Sine- and Cosine transform is

[TABLE]

11.2.1 The heat equation on a semi-infinite interval

We consider the equation

[TABLE]

with initial condition

[TABLE]

and boundary condition

[TABLE]

Because of the Diricelet boundary condition at $x=0$ it is appropriate to use the Fourier sine transform on the $x$-variable.

Taking the Fourier Sine transform of (683) and make use properties (678) and (685), we get

[TABLE]

with initial condition

[TABLE]

where

[TABLE]

The initial value problem (686), (687) is standard with solution

[TABLE]

Taking the inverse Fourier Sine transform (610) we get the solution in the form

[TABLE]

Using the addition formulas for cosine we have

[TABLE]

But then

[TABLE]

where we have used quantities defined in (605) and (610). We also have

[TABLE]

so the solution can compactly be written as

[TABLE]

Let us look at a couple of simple cases.

For the first case we assume $g(t)=0$, so the boundary condition is homogeneous, and that the initial temperature profile is uniform $u(x,0)=f(x)=u_{0}$. Then the second term in (695) vanishes and

[TABLE]

where $\operatorname{erf}(\xi)$ is the error function

[TABLE]

We observe that

[TABLE]

Using the properties of the error function we find that

[TABLE]

[TABLE]

[TABLE]

Thus $u(x,t)$ satisfies all requirements of our problem even if $u(x,0)=u_{0}$ does not have a Fourier sine transform. One can show that the solution (697) is the limit of a sequence of problems that all have initial data for which the transform exist.

For the second case we don’t assume that $f(x)$ is constant, but merely that it is bounded

[TABLE]

We let the boundary data at $x=0$ be given by an arbitrary time dependent function $g(t)$. Observe that for this case, the first term in the solution formula (695) approach zero for large $t$

[TABLE]

Therefore as $t\to\infty$ the solution $u(x,t)$ approaches

[TABLE]

This solution is called a steady state. It is not time invariant but represent the solution for times so large that the initial data has been dissipated. For such large times the solution only depends on the boundary condition, or drive at the boundary. Recall that an analogy behavior is displayed by a driven, damped harmonic oscillator.

11.3 The Laplace transform

The last transform we will discuss in this class is the Laplace transform. It is, as we will see, well suited for solving initial value problems. We will introduce the Laplace transform using the Fourier transform.

The Laplace transform is concerned with functions $f(t)$ which are zero for $t<0$. The function $f(t)$ can be quite wild for $t>0$, but we will set a bound on how fast it can grow when $t\to\infty$.

We will therefore require $f(t)$ to be of exponential order as $t\to\infty$, This means that there are constants $c>0$, $k$ and $M$ such that

[TABLE]

Let now $\alpha$ be a constant such that $\alpha>k$. Then the Fourier transform of $e^{-\alpha t}f(t)$ will exist and is given by

[TABLE]

where the fact that $f(t)=0$ for $t\leq 0$ has been taken into account. For example if $f(t)=1$ we get

[TABLE]

In the Fourier transform we know that $\hat{\lambda}$ is a real number related to the spectrum of a certain differential operator. Observe, however, that the right hand side of (707) can be evaluated for complex $\hat{\lambda}$ also. For example

[TABLE]

Also observe that (707) can not be evaluated for all complex numbers, in fact for $\hat{\lambda}=-i\alpha$ the denominator in (707) is zero. However, this is the only problematic point and $\hat{F}(\hat{\lambda})$ can be extended to the whole complex plane otherwise

[TABLE]

In a similar way we will now prove that $\hat{F}(\lambda)$ can be extended to complex values for any $f(t)$ of exponential order.

Evaluating (706) at $\hat{\lambda}=\hat{\lambda}_{r}+i\hat{\lambda}_{i}\in\mathbb{C}$ gives us

[TABLE]

For this to make sense the integral has to converge. But

[TABLE]

and the last integral is finite only if

[TABLE]

So $\hat{F}(\lambda)$ can not be defined in the whole complex plane but only in a domain defined by (712).

Below the horizontal line $\hat{\lambda}=-i(\alpha-k)$ in the complex plane, we expect the function $\hat{F}(\hat{\lambda})$ to have singularities. For the special choice

[TABLE]

this was exactly what occured.

Change variables in the complex plane to $\lambda$ through

[TABLE]

and define a function $F(\lambda)$ by

[TABLE]

From (710) we then get

[TABLE]

For a function of exponential order $k$, $F(\lambda)$ is defined for complex arguments $\lambda=\lambda_{r}+i\lambda_{i}$, with

[TABLE]

$F(\lambda)$ is the Laplace transform of $f(t)$. Using the Fourier inversion formula we have

[TABLE]

Introducing $\hat{\lambda}=i(\lambda-\alpha)$ gives us $\mathrm{d}\hat{\lambda}=i\,\mathrm{d}\lambda$ so that (718) implies

[TABLE]

where $L$ is a parameterized curve

[TABLE]

The curve $L$ is indicated in Figure 60.

The formula (719) is the inversion formula for the Laplace transform. The inversion formula is a very effective and flexible way to calculate the inverse Laplace transform, but using it requires residue calculus from complex analysis, so we will not make further use of formula (719) in this class. Inverse Laplace transforms can in many common situations be found using a combination of formal properties of the inverse Laplace transform and tables of inverse transforms. Lists of such formal properties can be found in mathematical tables or on the web. Here we will only display the most important property we need when we solve initial value problems for differential equations using the Laplace transform

[TABLE]

This formula is easy to prove using integration by parts. From (721) it is clear that the $n$ initial conditions for a $n^{\text{th}}$ order differential equation can be included in a natural way. The Laplace transform is therefore very useful for solving initial value problem for both ODEs and PDEs. Here we will apply it to initial boundary value problems for PDEs.

We will now apply the Laplace transform to two initial value problems for PDEs.

11.3.1 The wave equation on a half interval

[TABLE]

Taking the Laplace transform of the equation we get

[TABLE]

Using the initial conditions we get

[TABLE]

The general solution to the family of ordinary differential equations (724) is

[TABLE]

Recall that we evaluate (725) in the right halfspace only, so $\lambda_{r}>0$. In order to satisfy the boundary condition

[TABLE]

We must clearly assume that

[TABLE]

Finally applying the boundary condition at $x=0$ we get

[TABLE]

We now must take the inverse Laplace transform of (728).

Observe that if $H(\xi)$ is the Heaviside function

[TABLE]

Then

[TABLE]

This is equal to (728) with $b=\frac{x}{c}$. Thus the inverse Laplace transform of (728) must be

[TABLE]

and we have solved the initial value problem.

11.3.2 The heat equation on a finite interval

Consider the problem

[TABLE]

Taking the Laplace transform of (732) we get

[TABLE]

This is now a boundary value problem. Equation (734) can be solved using variation of parameters from the theory of ODEs. Fitting the boundary conditions to the variation of parameter solution gives

[TABLE]

Finding the inverse transform of (735) is not easy, but we have a secret weapon. In the inversion formula for the Laplace transform (719) we can place the integration curve $L$ anywhere as long as $\lambda_{r}>k$. Thus we may assume in (735) that $|\lambda|$ is as large as we want.

For such $\lambda$ we have

[TABLE]

The last equality sign in (736) is valid because for $|\lambda|$ large enough \big{|}e^{-2l\sqrt{\lambda}/c}\big{|}<1 and therefore

[TABLE]

The hyperbolic sine under the integral signs can also be expressed in terms of exponentials. Thus we end up with an infinite sum where each term is an exponential.

From a table of Laplace transforms we have

[TABLE]

Inserting (738) into the infinite sum, formally interchanging summation and integration gives us finally

[TABLE]

where $G(y,t)$ is the fundamental solution to the heat equation. We have previously found the solution to (732) using a Fourier sine series

[TABLE]

Problem (732) has a unique solution so even if (739) and (740) appear to be very different they must in fact be identical.

The formulas can be proven to be identical using the Poisson summation formula[1].

The formulas are mathematically equivalent, but are useful for evaluating $u(x,t)$ for different ranges of $t$. Because of the exponential term in (740), evaluating $u(x,t)$ for large $t$ using this formula requires only a few terms.

However for small $t$, $j$ must be large before the exponential term kicks in. Thus we require many terms in (740) in order to evaluate $u(x,t)$ for small $t$.

On the other hand, for small $t$, the fundamental solution $G(y,t)$ is sharply peaked around $y=0$. Thus for small $t$, only a few terms in (739) contribute to $u(x,t)$. Thus for small $t$ (739) is fast to evaluate.

In general it is not easy to invert integral transforms. There are however powerful analytic approximation methods that can be brought to bear on this problem.

Observe that the formula for the Laplace transform

[TABLE]

clearly shows that $F(\lambda)$ for large $|\lambda|$ depends mainly on $f(t)$ for small $t$. This can be made more precise. From the body of results called Abelian asymptotic theory we have the following result:

If

[TABLE]

then

[TABLE]

The converse in also true, (743) implies (742). Thus if we are only interested in $f(t)$ for small $t$ we can find the expression for $f(t)$ by expanding $F(\lambda)$ in inverse powers of $\lambda$. This expansion gives us the coefficients $\alpha_{n}$ which we then insert into (742).

There are also results relating the behavior of $f(t)$ for large $t$ to the behavior of $F(\lambda)$ for small $\lambda$. This body of results is called Tauherian asymptotic theory. Describing the simplest results from this theory requires however complex analysis so we will not do it here.

11.4 Group velocity 1

Let us consider the Klein-Gordon equation

[TABLE]

Initial data is of the form

[TABLE]

Recall that the Klein-Gordon equation is of dispersive type. We have found this previously using normal mode analysis. Taking the Fourier transform of the Klein-Gordon equation we get

[TABLE]

The general solution to this equation is

[TABLE]

where $\omega=\omega(\lambda)$ is the dispersion relation for the Klein-Gordon equation

[TABLE]

Using the inverse Fourier transform we get the following general solution to (744)

[TABLE]

The functions $F_{+}(\lambda)$ and $F_{-}(\lambda)$ are determined from initial data.

From (745) the Fourier transform of the initial data is

[TABLE]

Using (747) we get

[TABLE]

Solving (751) gives us

[TABLE]

(752) together with (749) solves the initial value problem (745) for the Klein-Gordon equation (744).

The solution (749) is clearly a superposition of normal modes

[TABLE]

The phase-speed for the normal modes are

[TABLE]

Recall that the phase speed of a normal mode is the speed we have to move at in order for the phase

[TABLE]

to stay fixed. Thus we have

[TABLE]

The Klein-Gordon equation has the same principle part as the wave equation. It thus have the same characteristic speed $\pm\gamma$.

As we have seen, discontinuities move at characteristic speed $\pm\gamma$ and $\pm\gamma$ represents the upper speed limit for the Klein-Gordon equation.

But from (754) we see that the phase speed is bigger than $\gamma$ for all modes. The physical relevance of the phase speed is therefore in question. Furthermore, noting that the solution to the Klein-Gordon equation is a superposition of normal modes which all move at different phase velocities, it is not at all clear that the phase speed is a relevant quantity for this equation.

In order to investigate this problem we will consider special initial conditions of the form

[TABLE]

where $k_{0}$ is fixed. The function $e^{ik_{0}x}$ is clearly periodic of period

[TABLE]

In (758) we assume that $0<\epsilon<<1$ is a number so small that

[TABLE]

Then the amplitudes $p(\epsilon x)$ and $q(\epsilon x)$ vary very little over a period $L_{0}$. The initial conditions (758) are what we call wave packets. These are functions that locally look like plane waves but that has a (complex) amplitude that varies over scales much larger than $L_{0}$.

Using simple properties of the Fourier transform we find easily from (749) that

[TABLE]

where $P$ and $Q$ are the Fourier transform of $p$ and $q$.

From Figure 62 it is clear that in Fourier space, $F(\lambda)$ and $G(\lambda)$ are narrow "spikes" of width $\epsilon$, height $\frac{1}{\epsilon}$ and centered at $\lambda=k_{0}$.

The functions $F_{+}(\lambda)$, $F_{-}(\lambda)$ will also be of this form

[TABLE]

for some functions $A_{\pm}$.

The solution of the Klein-Gordon equation consists of two terms of a very similar structure. The first term is

[TABLE]

Consider the following picture

From figure 63 it is clear that if the wavepacket is narrow enough, the dispersion relation $\omega(\lambda)$ is well approximated by its first order Taylor polynomial within the support of $F_{+}(\lambda)$. Thus we have

[TABLE]

Changing integration variable in the integral

[TABLE]

gives us

[TABLE]

Defining

[TABLE]

we get the final expression

[TABLE]

We observe that the envelope of the wave packet, $a_{+}$, translates without changing shape at a speed

[TABLE]

$v_{g}$ is the group velocity for the wave packet. The group velocity, and not the phase velocity, is the relevant speed for a wave packet. This conclusion holds true for all equations of dispersive type.

For the Klein-Gordon equation

[TABLE]

The group velocity is always smaller than the characteristic velocity $\gamma$ for the Klein-Gordon equation

Note that the statement that the wave packet move at the group velocity is only approximately true. A closer investigation shows that eventually the wave packet spreads out and disperse.

We will shortly see that the group velocity plays a pivotal role in the description of general solutions to dispersive equations. However, before we do this we must describe a tool of outmost importance. This is the method of stationary phase.

11.5 Method of stationary phase

The method of stationary phase is part of a very large and important domain of applied mathematics called asymptotic approximations, or perturbation methods.

The method of stationary phase seeks to approximate integrals of complex integrands, where the phase of the integrand has a component that varies very quickly.

Let us consider the integral

[TABLE]

where $\varphi(t)$ is a real function and $f(t)$ can be a complex function. Both functions are assumed to be smooth.

We want to investigate the integral for large $k$. The formulas we derive will hold in the limit $k\to\infty$. In any actual application of the formulas, $k$ is a fixed number, that is large, but not infinite. The errors we make in applying the asymptotic formulas using a given finite $k$, are not so easy to estimate. That is why we call it method of stationary phase rather than theory of stationary phase: Proofs are few and far between in this domain of applied mathematics. The method of stationary phase and asymptotic methods in general tend however to be surprisingly accurate even for values of $k$ that are nowhere "close" to infinity.

Observe that when $k$ is very large the phase term

[TABLE]

performs very fast oscillations. Since $f(t)$ is assumed to be smooth, the phase (772) will oscillate on a scale much faster than the scale at which the amplitude, $f(t)$, varies.

In this situation, as the figure illustrates, the positive and negative contributions in the integral tends to cancel each other so that we get

[TABLE]

This argument can be made into precise mathematics using the Riemann-Lebesques lemma from complex analysis, but we will not pursue this any further here.

Let $t=t_{0}$. Close to $t_{0}$ we have

[TABLE]

The largest contribution to $I(k)$ comes close to points where

[TABLE]

because around such points the phase varies more slowly than around points where $\varphi^{\prime}(t_{0})\neq 0$ for $k$ large.

A point $t_{0}$ such that (775) hold is called a stationary point. If $t=t_{0}$ is a stationary point, the contribution to $I(k)$, which we call $I(k;t_{0})$, is approximately given by

[TABLE]

where we have used the fact that for a small enough interval $(t_{0}-\epsilon,t_{0}+\epsilon)$, $f(t)$ varies very little around the value $f(t_{0})$. This is where we use the smoothness of $f(t)$. If $f(t)$ has a singularity at $t=t_{0}$ the approximation (776) is not valid.

If $k$ is large enough the integrand will oscillate very fast beyond $(t_{0}-\epsilon,t_{0}+\epsilon)$, so by the same argument as the one leading up to (773), we can conclude that the contribution from the integral

[TABLE]

is negligible with respect to the contribution from the integral

[TABLE]

Thus we can without loss of generality write (776) in the form

[TABLE]

The integral in (779) can be evaluated exactly using complex analysis

[TABLE]

If you have several stationary points, each point gives a contribution of type (780) so that finally

[TABLE]

where $\varphi^{\prime}(t_{j})=0,\,j=1,2,\dots,n.$

The method of stationary phase can be extended to integrals over finite domains and to multidimensional integrals.

The importance of the method of stationary phase can hardly be overstated. For example, it is by using this method that we understand how our macroscopic everyday world, ruled by classical physics, arise from a substratum ruled by the laws of quantum physics.

11.6 Group velocity 2

Let us now return to the Klein-Gordon equation (744), whose solution satisfying the general initial data (745) is

[TABLE]

where $F_{+}$ and $F_{-}$ are defined in terms of the initial data in (752). Both terms in (782) have a similar structure, so let us focus our attention on the first of them

[TABLE]

We want to investigate this expression for both $x$ and $t$ large. For such large values of $x$ and $t$ the phase varies very quickly as a function of $\lambda$ and we will, according to the method of stationary phase, get the largest contributions to $I(x,t)$ close to values of $\lambda$ where the phase is stationary. Thus we seek $\lambda$ such that

[TABLE]

Thus for given $x$ and $t$ large we find the stationary points by solving the equation

[TABLE]

For the Klein-Gordon equation we have from (748)

[TABLE]

Thus the equation we must solve is

[TABLE]

Clearly this equation only has solutions if

[TABLE]

But the light cone for the Klein-Gordon equation is $|x|=\gamma|t|$,

so stationary points only exists for $(x,t)$ inside the light cone. Outside the light cone there are no stationary points and according to the Riemann Lebesques lemma $I(x,t)\approx 0$.

Assuming \big{|}\frac{x}{t}\big{|}<\gamma we solve (789) and find

[TABLE]

We do not get $\pm$ in this formula because according to (787) $\lambda$ and $\frac{x}{t}$ have the same sign.

Observe that

[TABLE]

so that

[TABLE]

and therefore

[TABLE]

Thus the stationary phase formula (780) gives us

[TABLE]

where $H(\lambda)=F_{+}(\lambda)(\omega^{\prime\prime}(\lambda))^{-\frac{1}{2}}$ and where the phase in (783) has been written as

[TABLE]

so $k$ in the stationary phase formula is identified with $t\gg 1$.

So what have we found? Observe that

[TABLE]

We have the picture

According to (795) and (779) for large $x$ and $t$ along the line $x=\omega^{\prime}(\lambda)t$, the wave field take the form of a plane wave with frequency $\omega(\lambda)$, wavelength $\frac{2\pi}{\lambda}$ and amplitude

[TABLE]

The amplitude of the plane wave will decay as $\frac{1}{\sqrt{t}}$. Or to put it another way; if we shift to a moving frame moving at speed $\omega^{\prime}(\lambda)$ with respect to the lab frame we would see a plane wave whose amplitude decays as $\frac{1}{\sqrt{t}}$ and whose frequency $\omega(\lambda)$ and wavelength $\frac{2\pi}{\lambda}$ will depend on our speed. Observe that $\omega^{\prime}(\lambda)$ is just the group velocity that we introduced in (769).

One can prove that for the Klein-Gordon equation, and also for many other dispersive equations, the energy in the wave field is transported at the group velocity.

The group velocity is the relevant velocity for dispersive wave equations. It is also a concept whose importance for the description of wave phenomena in general hardly can be overstated.

12 Projects

12.1 Project1

In this project the aim is to develop and implement on a computer a finite difference algoritm for the solution of the 1D wave equation. Deriving the finite difference algoritm and implementing it in a computer code is not really the challenging and time consuming step here. The real challenge is to make a convincing case for the correctness of the computer implementation. As I have stressed in my lectures we should never, ever, trust anything that comes out of a computer without extensive testing of the implemented code.

So how do we proceed to build the needed trust in our implementation? The best way to do this is to run the code in a situation where an analytic solution is known. Most of the time an analytic solution also needs to be implemented numerically, it typically include infinite sums of functions and integrals that can not be solved in terms of elementary functions. So what we are really doing is to compare the output of two different computer implementations for the same solution. If this is to make sense and build confidence, it is of paramount importance that the two codes implementing the finite difference algoritm and the analytical solution should be as independent as humanly possible. Under such circumstances closely similar output from the two codes will build confidence that our finite difference solution is correct.

So the challenge then is to find analytical solutions to the problem at hand, in our particular case this is the 1D wave equation subject to some boundary conditions that we will specify shortly. It is typically the case that we will not be able to find analytic solutions to the actual problem we are trying to solve but to a related problem or a special case of the problem where the domain, boundary conditions and/or equation is changed somewhat or fixed by some specific choises. After all, if we could find an analytic solution to the actual problem it would not be much point in spending a lot of effort in trying to solve it using finite differences. It is however important that changing from the actual problem to the related problem does not involve major modifications to the core of the finite difference code. The reason for this is that such major modifications are their own source of error and as a consequence our test will not be very convincing as a test of the actual problem.

So what is the actual problem then? Here it is:

[TABLE]

a)   Implement a finite difference method for this problem. The code should be flexible enough to handle arbitrary choises for the functions $h$,$k\,$,$\varphi$ and $\psi$.

b)   Show that

[TABLE]

is a solution to the problem

[TABLE]

In your finite difference code choose $l$ and $\varphi$ in such a way that the interval is much larger than the support of $\varphi$.( By the support of $\varphi$ I mean the set of points where the function is nonzero). A Gaussian function of the form

[TABLE]

is a good choise if the parameters $a$, $b$ and $x_{0}$ are choosen in the right way. Also choose $h=k=0$. As long as the support of the solution $u$ is well away from the boundary of the interval $[0,l]$ our exact solution and the finite difference solution should match. This is our first test.

c)    Our second test involve a rather general idea called an "artificial source". We consider the modified problem

[TABLE]

This equation would appear as the evolution equation for a string clamped at both ends and where there is a precribed force acting along the string. This kind of prescribed influence is called a "source". In the artificial source test we turn the problem around and rather assume a solution of a particular type and then compute what the corresponding source must be.

We thus pick a function $u(x,t)$ and calculate the source $\rho(x,t)$, boundary functions $h(t),k(t)$ and initial data $\varphi(x)$ and $\psi(x)$ so that $u(x,t)$ is in fact a solution to our modified problem. Thus

[TABLE]

We now have a analytic solution to our modified problem and can compare this solution to the solution produced by our finite difference code. You are free to choose the function $u(x,t)$ for yourself but here is a possible family of choises

[TABLE]

This is a family of Gaussian functions whose center move back and forth in a periodic manner determined by the parameters $x_{0}$ and $\omega$. Compare the numerical solution to the exact solution for several choises of parameters. Make sure that at least one of the test involve the boundary conditions in a significant way. By this I mean that $h(t)$ and $k(t)$ should be nonzero for some values of $t$.

d)   Now use you code to solve the problem

[TABLE]

This model what happend when you periodically move the left end of the string,whose right end is clamped, up and down. What do you expect to happend? Do your code verify your intuition?

12.2 Project2

In this project the aim is to develop a stable numerical scheme for the initial value problem

[TABLE]

where

[TABLE]

Remark: Since the equation is first order in time, one can actually not pose boundary conditions at both endpoints. Thus we should only use the condition at one end point, for example $u(-L,t)=0$. However since we are using center difference we really also need a boundary condition at $x=L$. In this problem I ask you to use the condition $u(-L,t)=0$. Using these two conditions you can construct a finite difference scheme using center difference in space. However, and this is of outmost importance, you must make sure that the solution does not become nonzero in the region close to the boundaries, if it does the scheme will go unstable no matter what we do. But this restriction is really ok since we really are trying to simulate the dynamics on an infinite domain, and thus the computations boundaries at $x=-L,L$ should be "invisible" at all time.

a)

Derive a numerical scheme for equation (799) using center differences for the time and space derivatives. Find the Von-Neuman condition for numerical stability of your scheme for the case when $a=a_{0},b=b_{0}$ are constants and $\rho=0$.

b)

Verify the numerical stability condition by running the your numerical scheme for some initial condition of your choosing. Run with grid parameters both in the Von Neumann stable and unstable regime Since our numerical scheme is meant to be a good approximation to solutions of equation (799) on an infinite interval, you must make sure that the initial condition is very small or zero close to the end points of the interval $[-L,L]$. A well localized Gaussian function will do the job.

c)

Solve the initial value problem

[TABLE]

using the method of characteristics. Plot some of the base characteristic curves in the $x-t$ plane and also plot the solution $u(x,t)$ as a function of $x$ for several $t>0$.

d)

Solve the initial value problem from c) using your finite difference scheme. Compare with the solution exact solution from c) by plotting them on top of each other for several times.

e)

Let the following quasi linear initial value problem be given

[TABLE]

Find exact formulas for the shock time for solutions to the initial value problem from e) corresponding to the initial data

1. i)

[TABLE] 2. ii)

[TABLE] 3. iii)

[TABLE]

Test the validity of your finite difference code, and exact formulas for the shock time, by running your code with the given initial conditions and verifying that the graph of the numerical solutions become approximately vertical for some $x$ as the shock time is approached.

12.3 Project3

In this project we are going to solve a heat equation numerically using the finite difference method and the finite Fourier transform.

[TABLE]

a)

Implement a finite difference method for this problem. The code should be flexible enough to handle arbitrary choises for the functions $f$, $g\,$, $\rho$ and $\varphi$.

b)

The function

[TABLE]

satisfy the boundary conditions. Calculate the corresponding artificial source $\rho$ and initial condition $\varphi$. Now run your finite difference code with the calculated $\rho$ and $\varphi$ and compare the numerical solution to the exact solution $u_{e}(x,t)$. Do this with a couple of choises for the functions $f(t),g(t)$ and $h(x,t)$ and for several values of the parameters in the problem. For each choise plot the numerical and exact solution in the same plot and thereby verify that they coincide.

c)

Solve (800) using the finite Fourier transform. Before you apply the finite Fourier transform you should convert (800) into a problem with homogenous boundary conditions like in equation 4.6.45 in our textbook. Use the same choises for $f(t)$ and $g(t)$ that you used in problem b) and choose your own source. Compare the finite Fourier transform solution with the finite difference solution for some choises of initial conditions. Make sure that the initial conditions is consistent with the boundary conditions

[TABLE]

d)

Let the source in problem (800) be the nonlinear source from equation 4.7.1 in our textbook.

[TABLE]

Use the finite difference code to verify the linear and nonlinear stability results found in section 4.7 of our textbook. In particular you should shown that

1. i)

$u(x,t)=0$ is stable for $\widehat{\lambda}<1$ 2. ii)

$u(x,t)$ is unstable for $1<$ $\widehat{\lambda}<4$ and approach the solution 4.7.24 as time increase to infinity.

Run the finite difference code for values of $\widehat{\lambda}$ higher than 4 and plot the resulting solutions. Any conclusions?