Download Partial Differential Equations: Chapter 14 and more Lecture notes Differential Equations in PDF only on Docsity!
Chapter 14
Partial Differential Equations
Our intuition for ordinary differential equations generally stems from the time evolution of phys- ical systems. Equations like Newton’s second law determining the motion of physical objects over time dominate the literature on such initial value problems; additional examples come from chemical concentrations reacting over time, populations of predators and prey interacting from season to season, and so on. In each case, the initial configuration—e.g. the positions and veloc- ities of particles in a system at time zero—are known, and the task is to predict behavior as time progresses. In this chapter, however, we entertain the possibility of coupling relationships between different derivatives of a function. It is not difficult to find examples where this coupling is necessary. For instance, when simulating smoke or gases quantities like “pressure gradients,” the derivative of the pressure of a gas in space, figure into how the gas moves over time; this structure is reasonable since gas naturally diffuses from high-pressure regions to low-pressure regions. In image process- ing, derivatives couple even more naturally, since measurements about images tend to occur in the x and y directions simultaneously. Equations coupling together derivatives of functions are known as partial differential equations. They are the subject of a rich but strongly nuanced theory worthy of larger-scale treatment, so our goal here will be to summarize key ideas and provide sufficient material to solve problems commonly appearing in practice.
14.1 Motivation
Partial differential equations (PDEs) arise when the unknown is some function f : R n^ → R m. We are given one or more relationship between the partial derivatives of f , and the goal is to find an f that satisfies the criteria. PDEs appear in nearly any branch of applied mathematics, and we list just a few below. As an aside, before introducing specific PDEs we should introduce some notation. In particu- lar, there are a few combinations of partial derivatives that appear often in the world of PDEs. If f : R^3 → R and ~v : R^3 → R^3 , then the following operators are worth remembering:
Gradient: ∇ f ≡
∂ f ∂ x 1
∂ f ∂ x 2
∂ f ∂ x 3
Divergence: ∇ · ~v ≡ ∂ v 1 ∂ x 1
∂ v 2 ∂ x 2
∂ v 3 ∂ x 3 Curl: ∇ × ~v ≡
∂ v 3 ∂ x 2
∂ v 2 ∂ x 3
∂ v 1 ∂ x 3
∂ v 3 ∂ x 1
∂ v 2 ∂ x 1
∂ v 1 ∂ x 2
Laplacian: ∇^2 f ≡
∂^2 f ∂ x^21
∂^2 f ∂ x^22
∂^2 f ∂ x^23
Example 14.1 (Fluid simulation). The flow of fluids like water and smoke is governed by the Navier- Stokes equations, a system of PDEs in many variables. In particular, suppose a fluid is moving in some region Ω ⊆ R^3. We define the following variables, illustrated in Figure NUMBER:
- t ∈ [0, ∞): Time
- ~v(t) : Ω → R^3 : The velocity of the fluid
- ρ (t) : Ω → R : The density of the fluid
- p(t) : Ω → R : The pressure of the fluid
- ~f (t) : Ω → R^3 : External forces like gravity on the fluid
If the fluid has viscosity μ , then if we assume it is incompressible the Navier-Stokes equations state:
ρ
∂ ~v ∂ t
= −∇p + μ ∇^2 ~v + ~f
Here, ∇^2 ~v = ∂^2 v 1 / ∂ x^21 + ∂^2 v 2 / ∂ x^22 + ∂^2 v 3 / ∂ x^23 ; we think of the gradient ∇ as a gradient in space rather than
time, i.e. ∇ f = ( (^) ∂∂ x^ f 1 , (^) ∂∂ x^ f 2 , (^) ∂∂ x^ f 3 ). This system of equations determines the time dynamics of fluid motion and actually can be constructed by applying Newton’s second law to tracking “particles” of fluid. Its statement, however, involves not only derivatives in time (^) ∂∂ t but also derivatives in space ∇, making it a PDE.
Example 14.2 (Maxwell’s equations). Maxwell’s equations determine the interaction of electric fields ~E and magnetic fields ~B over time. As with the Navier-Stokes equations, we think of the gradient, divergence, and curl as taking partial derivatives in space (and not time t). Then, Maxwell’s system (in “strong” form) can be written:
Gauss’s law for electric fields: ∇ · ~E = ρ ε 0 Gauss’s law for magnetism: ∇ · ~B = 0
Faraday’s law: ∇ × ~E = −
∂ ~B
∂ t
Amp`ere’s law: ∇ × ~B = μ 0
~J + ε 0^ ∂
~E
∂ t
Here, ε 0 and μ 0 are physical constants and ~J encodes the density of electrical current. Just like the Navier- Stokes equations, Maxwell’s equations related derivatives of physical quantities in time t to their derivatives in space (given by curl and divergence terms).
can be solved using the following PDE:
∇^2 f (~x) = 0 f (~x) = g(~x) ∀~x ∈ ∂ Ω
This equation is known as Laplace’s equation, and it can be solved using sparse positive definite linear methods like what we covered in Chapter 10. As we have seen, it can be applied to interpolation problems for irregular domains Ω; furthermore, E[ f ] can be augmented to measure other properties of f , e.g. how well f approximates some noisy function f 0 , to derive related PDEs by paralleling the argument above.
Example 14.4 (Eikonal equation). Suppose Ω ⊆ R n^ is some closed region of space. Then, we could take d(~x) to be a function measuring the distance from some point ~x 0 to ~x completely within Ω. When Ω is convex, we can write d in closed form: d(~x) = ‖~x − ~x 0 ‖ 2.
As illustrated in Figure NUMBER, however, if Ω is non-convex or a more complicated domain like a surface, distances become more complicated to compute. In this case, distance functions d satisfy the localized condition known as the eikonal equation:
‖∇d‖ 2 = 1.
If we can compute it, d can be used for tasks like planning paths of robots by minimizing the distance they have to travel with the constraint that they only can move in Ω. Specialized algorithms known as fast marching methods are used to find estimates of d given ~x 0 and Ω by integrating the eikonal equation. This equation is nonlinear in the derivative ∇d, so integration methods for this equation are somewhat specialized, and proof of their effectiveness is complex. Interestingly but unsurprisingly, many algorithms for solving the eikonal equation have structure similar to Dijkstra’s algorithm for computing shortest paths along graphs.
Example 14.5 (Harmonic analysis). Different objects respond differently to vibrations, and in large part these responses are functions of the geometry of the objects. For example, cellos and pianos can play the same note, but even an inexperienced musician easily can distinguish between the sounds they make. From a mathematical standpoint, we can take Ω ⊆ R^3 to be a shape represented either as a surface or a volume. If we clamp the edges of the shape, then its frequency spectrum is given by solutions of the following differential eigenvalue problem:
∇^2 f = λ f f (x) = 0 ∀x ∈ ∂ Ω,
where ∇^2 is the Laplacian of Ω and ∂ Ω is the boundary of Ω. Figure NUMBER shows examples of these functions on different domains Ω. It is easy to check that sin kx solves this problem when Ω is the interval [0, 2 π ], for k ∈ Z. In particular, the Laplacian in one dimension is ∂^2 / ∂ x^2 , and thus we can check:
∂^2 ∂ x^2
sin kx =
∂ x
k cos kx
= −k^2 sin kx sin k · 0 = 0 sin k · 2 π = 0
Thus, the eigenfunctions are sin kx with eigenvalues −k^2.
14.2 Basic definitions
Using the notation of CITE, we will assume that our unknown is some function f : R n^ → R. For equations of up to three variables, we will use subscript notation to denote partial derivatives:
fx ≡
∂ f ∂ x
fy ≡ ∂ f ∂ y
fxy ≡
∂^2 f ∂ x ∂ y
and so on. Partial derivatives usually are stated as relationships between two or more derivatives of f , as in the following:
- Linear, homogeneous: fxx + fxy − fy = 0
- Linear: fxx − y fyy + f = xy^2
- Nonlinear: f (^) xx^2 = fxy
Generally, we really wish to find f : Ω → R for some Ω ⊆ R n. Just as ODEs were stated as initial value problems, we will state most PDEs as boundary value problems. That is, our job will be to fill in f in the interior of Ω given values on its boundary. In fact, we can think of the ODE initial value problem this way: the domain is Ω = [0, ∞), with boundary ∂ Ω = { 0 }, which is where we provide input data. Figure NUMBER illustrates more complex examples. Boundary conditions for these problems take many forms:
- Dirichlet conditions simply specify the value of f (~x) on ∂ Ω
- Neumann conditions specify the derivatives of f (~x) on ∂ Ω
- Mixed or Robin conditions combine these two
14.3 Model Equations
Recall from the previous chapter that we were able to understand many properties of ODEs by ex- amining a model equation y′^ = ay. We can attempt to pursue a similar approach for PDEs, although we will find that the story is more nuanced when derivatives are linked together. As with the model equation for ODEs, we will study the single-variable linear case. We also will restrict ourselves to second-order systems, that is, systems containing at most the second derivative of u; the model ODE was first-order but here we need at least two orders to study how derivatives interact in a nontrivial way. A linear second-order PDE has the following general form:
∑ ij
aij ∂ f ∂ xi ∂ xj
bi ∂ f ∂ xi
- Dirichlet conditions for this equation simply specify f (a) and f (b); there is obviously a unique line that goes through (a, f (a)) and (b, f (b)), which provides the solution to the equation.
- Neumann conditions would specify f ′(a) and f ′(b). But, f ′(a) = f ′(b) = α for f (x) = α x + β. In this way, boundary values for Neumann problems can be subject to compatibility conditions needed to admit admit a solution. Furthermore, the choice of β does not affect the boundary conditions, so when they are satisfied the solution is not unique.
14.3.2 Parabolic PDEs
Continuing to parallel the structure of linear algebra, positive semidefinite systems of equations are only slightly more difficult to deal with than positive definite ones. In particular, positive semidefinite matrices admit a null space which must be dealt with carefully, but in the remaining directions the matrices behave the same as the definite case. The heat equation is the model parabolic PDE. Suppose f (0; x, y) is a distribution of heat in some region Ω ⊆ R^2 at time t = 0. Then, the heat equation determines how the heat diffuses over time t as a function f (t; x, y): ∂ f ∂ t = α ∇^2 f ,
where α > 0 and we once again think of ∇^2 as the Laplacian in the space variables x and y, that is, ∇^2 = ∂^2 / ∂ x^2 + ∂^2 / ∂ y^2. This equation must be parabolic, since there is the same coefficient α in front of fxx and fyy, but ftt does not figure into the equation. Figure NUMBER illustrates a phenomenological interpretation of the heat equation. We can think of ∇^2 f as measuring the convexity of f , as in Figure NUMBER(a). Thus, the heat equa- tion increases u with time when its value is “cupped” upward, and decreases f otherwise. This negative feedback is stable and leads to equilibrium as t → ∞. There are two boundary conditions needed for the heat equation, both of which come with straightforward physical interpretations:
- The distribution of heat f (0; x, y) at time t = 0 at all points (x, y) ∈ Ω
- Behavior of f when t > 0 at points (x, y) ∈ ∂ Ω. These boundary conditions describe behav- ior at the boundary of the domain. Dirichlet conditions here provide f (t; x, y) for all t ≥ 0 and (x, y) ∈ ∂ Ω, corresponding to the situation in which an outside agent fixes the temper- atures at the boundary of the domain. These conditions might occur if Ω is a piece of foil sitting next to a heat source whose temperature is not significantly affected by the foil like a large refrigerator or oven. Contrastingly, Neumann conditions specify the derivative of f in the direction normal to the boundary ∂ Ω, as in Figure NUMBER; they correspond to fixing the flux of heat out of Ω caused by different types of insulation.
14.3.3 Hyperbolic PDEs
The final model equation is the wave equation, corresponding to the indefinite matrix case:
∂^2 f ∂ t^2 − c^2 ∇^2 f = 0
The wave equation is hyperbolic because the second derivative in time has opposite sign from the two spatial derivatives. This equation determines the motion of waves across an elastic medium like a rubber sheet; for example, it can be derived by applying Newton’s second law to points on a piece of elastic, where x and y are positions on the sheet and f (t; x, y) is the height of the piece of elastic at time t. Figure NUMBER illustrates a one-dimensional solution of the wave equation. Wave behavior contrasts considerably with heat diffusion in that as t → ∞ energy may not diffuse. In particular, waves can bounce back and forth across a domain indefinitely. For this reason, we will see that implicit integration strategies may not be appropriate for integrating hyperbolic PDEs because they tend to damp out motion. Boundary conditions for the wave equation are similar to those of the heat equation, but now we must specify both f (0; x, y) and ft(0; x, y) at time zero:
- The conditions at t = 0 specify the position and velocity of the wave at the initial time.
- Boundary conditions on Ω determine what happens at the ends of the material. Dirichlet conditions correspond to fixing the sides of the wave, e.g. plucking a cello string, which is held flat at its two ends on the instrument. Neumann conditions correspond to leaving the ends of the wave untouched like the end of a whip.
14.4 Derivatives as Operators
In PDEs and elsewhere, we can think of derivatives as operators acting on functions the same way that matrices act on vectors. Our choice of notation often reflects this parallel: The derivative d f/dx looks like the product of an operator d/dx and a function f. In fact, differentiation is a linear operator just like matrix multiplication, since for all f , g : R → R and a, b ∈ R
d dx (a f (x) + bg(x)) = a
d dx f (x) + b
d dx g(x).
In fact, when we discretize PDEs for numerical solution, we can carry this analogy out com- pletely. For example, consider a function f on [0, 1] discretized using n + 1 evenly-spaced samples, as in Figure NUMBER. Recall that the space between two samples is h = 1 /n. In Chapter 12, we developed an approximation for the second derivative f ′′(x) :
f ′′(x) = f (x + h) − 2 f (x) + f (x − h) h^2
Suppose our n samples of f (x) on [0, 1] are y 0 ≡ f ( 0 ), y 1 ≡ f (h), y 2 ≡ f ( 2 h),... , yn = f (nh). Then, applying our formula above gives a strategy for approximating f ′′^ at each grid point:
y′′ k ≡
yk+ 1 − 2 yk + yk− 1 h^2
That is, the second derivative of a function on a grid of points can be computed using the 1— − 2—1 stencil illustrated in Figure NUMBER(a). One subtlety we did not address is what happens at y′′ 0 and y′′ n , since the formula above would require y− 1 and yn+ 1. In fact, this decision encodes the boundary conditions introduced in §14.2. Keeping in mind that y 0 = f ( 0 ) and yn = f ( 1 ), examples of possible boundary conditions for f ′ include:
algorithm in §13.4.2, one way to avoid these issues is to think of the derivatives as living on half gridpoints, illustrated in Figure NUMBER. In the one-dimensional case, this change corresponds to labeling the difference (^1) h (yk+ 1 − yk) as yk+ (^1) / 2. This technique of placing different derivatives on vertices, edges, and centers of grid cells is particularly common in fluid simulation, which maintains pressures, fluid velocities, and so on at locations that simplify calculations.
These subtleties aside, our main conclusion from this discussion is that if we discretize a func- tion f (~x) by keeping track of samples (xi, yi) then most reasonable approximations of derivatives of f will be computable as a product L~x for some matrix L. This observation completes the anal- ogy: “Derivatives act on functions as matrices act on vectors.” Or in standardized exam notation: Derivatives : Functions :: Matrices : Vectors
14.5 Solving PDEs Numerically
Much remains to be said about the theory of PDEs. Questions of existence and uniqueness as well as the possibility of characterizing solutions to assorted PDEs leads to nuanced discussions using advanced aspects of real analysis. While a complete understanding of these properties is needed to prove effectiveness of PDE discretizations rigorously, we already have enough to suggest a few techniques that are used in practice.
14.5.1 Solving Elliptic Equations
We already have done most of the work for solving elliptic PDEs in §14.4. In particular, suppose we wish to solve a linear elliptic PDE of the form L f = g. Here L is a differential operator; for example, to solve the Laplace’s equation we would take L ≡ ∇^2 , the Laplacian. Then, in §14. we showed that if we discretize f by taking a set of samples in a vector ~y with yi = f (xi), then a corresponding approximation of L f can be written L~y for some matrix L. If we also discretize g using samples in a vector~b, then solving the elliptic PDE L f = g is approximated by solving the linear system L~y = ~b.
Example 14.7 (Elliptic PDE discretization). Suppose we wish to approximate solutions to f ′′(x) = g(x) on [0, 1] with boundary conditions f ( 0 ) = f ( 1 ) = 0. We will approximate f (x) with a vector ~y ∈ R n sampling f as follows:
y 1 y 2 .. . yn
f (h) f ( 2 h) .. . f (nh)
where h = 1 /n+ 1. We do not add samples at x = 0 or x = 1 since the boundary conditions determine values there. We will use~b to hold an analogous set of values for g(x).
Given our boundary conditions, we discretize f ′′(x) as (^) h^12 L~y, where L is given by:
L ≡
Thus, our approximate solution to the PDE is given by ~y = h^2 L−^1 ~b.
Just as elliptic PDEs are the most straightforward PDEs to solve, their discretizations using matrices as in the above example are the most straightforward to solve. In fact, generally the discretization of an elliptic operator L is a positive definite and sparse matrix perfectly suited for the solution techniques derived in Chapter 10.
Example 14.8 (Elliptic operators as matrices). Consider the matrix L from Example 14.7. We can show L is negative definite (and hence the positive definite system −L~y = −h^2 ~b can be solved using conjugate gradients) by noticing that −L = D>^ D for the matrix D ∈ R (n+^1 )×n^ given by:
D =
This matrix is nothing more than the finite-differenced first derivative, so this observation parallels the fact that d^2 f/dx^2 = d/dx(d f/dx). Thus, ~x>^ L~x = −~x>^ D>^ D~x = −‖D~x‖^22 ≤ 0, showing L is negative semidefinite. It is easy to see D~x = 0 exactly when ~x = 0 , completing the proof that L is in fact negative definite.
14.5.2 Solving Parabolic and Hyperbolic Equations
Parabolic and hyperbolic equations generally introduce a time variable into the formulation, which also is differentiated but potentially to lower order. Since solutions to parabolic equations admit many stability properties, numerical techniques for dealing with this time variable often are stable and well-conditioned; contrastingly, more care must be taken to treat hyperbolic behavior and prevent dampening of motion over time.
Semidiscrete Methods Probably the most straightforward approach to solving simple time- dependent equations is to use a semidiscrete method. Here, we discretize the domain but not the time variable, leading to an ODE that can be solved using the methods of Chapter 13.
Example 14.9 (Semidiscrete heat equation). Consider the heat equation in one variable, given by ft = fxx, where f (t; x) represents the heat at position x and time t. As boundary data, the user provides a
Example 14.10 (Eigenfunctions of the Laplacian). Figure NUMBER shows eigenvectors of the matrix L from Example 14.7. Eigenvectors with low eigenvalues correspond to low-frequency functions on [0, 1] with values fixed on the endpoints and can be good approximations of f (x) when it is relatively smooth.
Fully Discrete Methods Alternatively, we might treat the space and time variables more demo- cratically and discretize them both simultaneously. This strategy yields a system of equations to solve more like §14.5.1. This method is easy to formulate by paralleling the elliptic case, but the resulting linear systems of equations can be large if dependence between time steps has a global reach.
Example 14.11 (Fully-discrete heat diffusion). Explicit, implicit, Crank-Nicolson. Not covered in CS 205A.
It is important to note that in the end, even semidiscrete methods can be considered fully discrete in that the time-stepping ODE method still discretizes the t variable; the difference is mostly for classification of how methods were derived. One advantage of semidiscrete techniques, however, is that they can adjust the time step for t depending on the current iterate, e.g. if objects are moving quickly in a physical simulation it might make sense to take more time steps and resolve this motion. Some methods even adjust the discretization of the domain of x values in case more resolution is needed near local discontinuities or other artifacts.
14.6 Method of Finite Elements
Not covered in 205A.
14.7 Examples in Practice
In lieu of a rigorous treatment of all commonly-used PDE techniques, in this section we provide examples of where they appear in practice in computer science.
14.7.1 Gradient Domain Image Processing
14.7.2 Edge-Preserving Filtering
14.7.3 Grid-Based Fluids
14.8 To Do
- More on existence/uniqueness
- CFL conditions
- Lax equivalence theorem
- Consistency, stability, and friends
14.9 Problems
- Show 1d Laplacian can be factored as D>^ D for first derivative matrix D
- Solve first-order PDE
