C ha pter 1

Intr od uction

Partial diﬀerential equations, or PDEs for short, are an important type of diﬀerential equation that

arise as natural mathematical models in many physical problems. They allow us to describe the

behavior of a system in terms of fu nctions that depend on multiple variables, such as time and space.

For example, consider the classical heat equation, which describes the distribution of heat in a

conducting material over time. This equation can be derived from Fourier's law of heat conduction,

which states that the rate of heat transfer through a material is proportional to the temperature

gradient. By applying this law to an inﬁnitesimal volume element in the material, one obtains a

PDE for the temperature distribution.

Similarly, the wave equation can be derived from Newton's second law of motion, which relates

the force acting on a body to its acceleration. By applying this law to a small element of a string

or other vibrating object, one obtains a PDE for the displacement of the element as a function of

time and position.

In bo th of these examples, the PDEs are derived from fundamental physical laws and provide a

mathematical description of the underlying physical phenomenon. By solving these equations, we

can make predictions ab out how the system will behave under diﬀerent conditions, such as changes

in tempe rature or initial conditions.

Our main focus is to introduce four imp ortant types of PDEs at an entry-level technicality,

without delving into their theoretical aspects. These include:

1. The transport equation, which describes how the density of a ﬂuid ﬂow moving in the plane or

space changes according to a given velocity. This PDE is of ﬁrst-order. The following ﬁgure

illustrates the variation of mass density of a ﬂuid ﬂow in the xy-plane with the velocity ﬁeld

V = (¡y; x).

2. The heat problem, which involves determining how the temperature of a conductive medium

changes in time, given an initial heat distribution along a conductive r od or continuum. This

problem leads to a second-order PDE. The ﬁgure below shows the evolution of temperature

over time in a unit disk.

3. The wave equation, which describes how an elastic string or a 2D membrane responds to an

initial disturbance, leading to oscillations and the propagation o f the disturbance. This PDE

is also of second-order. The ﬁgure below illustrates how an initial disturbance propagates as

a bilateral wave, splitting into two branches moving with a given velocity to the left and right.

2 Introduction

-10 -5 0 5 10

0.2

0.4

0.6

0.8

1.2

-10 -5 0 5 10

0.2

0.4

0.6

0.8

1.2

4. The Poisson equation, which arises in the context of electrostatics and involves ﬁnding the

potential function generated in space by an electric charge distributed in a domain Ω in the

plane or space. This problem is reduced to a second-order PDE.

1.1 From ODEs to PDEs

In the following discussion, we will explore the natural extension of well-known ordinary diﬀerential

equations to partial diﬀerential equations. We will focus on the derivation of some simple PDEs

and highlight how they arise from the ODE versions of the same problems. This will illustrate the

process of extending mathematical models of physical phenomena from one dimension to multiple

dimensions.

1.1.1 From mass-spring system to vibrating string

Remember that the dynamics of a single or ﬁnite set of quantities that depends on a single indepen-

dent variable leads to a single or set of diﬀerential equations called ordinary diﬀerential equations

(ODEs). An i mportant example of this type of equation is the harmonic oscillator: consider a set

of n point masses m₁; :::; mₙ that are connected in series to n springs with stiﬀnesses k₁; :::; kₙ, as

shown in the ﬁgure below.

The equation that describes the dynamics of m

is derived by Newton's second l aw f = ma as

= ¡k

+ k

¡ x

);

where ¡k₁x₁ is the Hoo k's force exerted by the spring k₁, and k₂(x₂¡x₁) is the force exerted by k₂

on m₁. Similarly, the equation that describes the dynamics of m₂ is:

= ¡k

¡ x

) + k

¡ x

1.1 From ODEs to PDEs 3

The equation for m

; :::; m

are derived in a similar manner. In this way, we obtain a system of n

second-order ODEs

x¨

= ¡k

+ k

¡ x

)

x¨

= ¡k

¡ x

) + k

¡ x

)

x¨

= ¡k

¡ x

n¡1

)

Now, consider a vibrating string along the x-axis, which can be modeled a s an inﬁnite mass-spring

system connected in series. For a control length δx with density ρ, the control mass is ρδx. Let

u(x; t) denote the position of this control length. This piece is under the forces of masses ρ(x ¡δx)δx,

which is posited at u(x + δx; t), and ρ(x + δx)δx at the position u(x + δx; t).

u(x + δx; t)u(x; t)u(x ¡ δx; t)

The force on the mass at x is

F = ¡

δx

[u(x; t) ¡ u (x ¡ δx; t)] +

δx

[u(x + δx; t) ¡ u(x; t)];

where k is the stiﬀness of the string, and according to the control length δx, the ratio

δx

reﬂects

the impact of this length on the motion of mass in δx. Therefore, we can write

ρδx u

(x; t) = ¡k

u(x; t) ¡ u(x ¡ δx; t)

δx

+ k

u(x + δx; t) ¡ u(x; t)

δx

;

where u

(x; t) :=

is the acceleration of the control length δx. By using the relations

u(x; t) ¡ u(x ¡ δx; t)

δx

≈ u

(x ¡ δx; t);

u(x + δx; t) ¡ u(x; t)

δx

≈ u

(x; t);

we obtain

ρ δx u

(x; t) ≈ ¡ku

(x ¡ δx; t) + ku

(x; t);

and by dividing by δx, we obtain

ρ u

(x; t) ≈ k

(x; t) ¡ u

(x ¡ δx; t)

δx

;

where the right-hand side approaches ku

when δx!0. Finally, we derive the following diﬀerential

equation for the oscillation of the vibrating string:

(x; t) =

x x

(x; t):

Exercise 1.1. For a mass-spring system consists of n masses m

; :::; m

connected to n springs k

; :::;k

in series ,

the total energy is deﬁned as

E(t) =

j=1

jx_

(t)j

j=2

(t) ¡x

j ¡1

(t)j

Veri fy that the system conserves the ener gy, that is,

= 0, and conclude E(t) = E(0).

Exercise 1.2. The conser vation of energy also applies to the wave equation of a string. Let us consid er a string

of length L, ﬁxed at two endpoints x = 0 and x = L. We denote the position of t he point x at time t by u(x; t).

The st ring forms a dyna mic system i f either its initi al ki netic or potential energy is nonzero. The i nitial kinetic

energy of the string can be de ﬁned as follows:

K(t): =

ρ ju

(x; t)j

dx:

4 Introduction

The Hook's stretch potential energy is

U(t) =

(x; t)j

dx;

Veri fy that the tota l energy E(t)=K(t)+U(t) o f the str ing is constant, i.e .,

=0, and conclude t hat E(t)=E(0).

Thus, in order to solve a wave equation, we requir e the initial disturbance u(x; 0) and in itial velocity u

(x; 0).

1.1.2 General remarks

As we saw from the above example, ODEs and PDEs diﬀer in some signiﬁcant aspects. We can

summarize the diﬀerences a s follows.

• The solution to an ODE is a function or a set of functions of only one independent variable,

whereas for a PDE, the solution depends on two or more independent variables. For the wave

equation, the solution u depends on two variables, x and t, and represents the position o f

point x at time t. For this reason, the derivatives in a partial diﬀerential equation are partial

derivatives instead of ordinary ones.

• While ODEs are related to pointwise quantities, PDEs are mathematical models of distrib-

uted systems or continua. For this reason, PDEs are sometimes considered as an inﬁnite set

of ODEs.

• The general solution of an ODE depends on one or more constant parameters. For example,

the solution of the harmonic oscillator x¨ +

x = 0 can be expressed as a linear combination

of two fundamental solutions fc os(!

t); sin(!

t)g, where !

is the natural frequency

of the single mass-spring system. The solution is given by

x(t)=c

cos(!

t)+c

sin(!

t);

where c

and c

are arbitrary constant parameters. In contrast, the solution of PDEs usually

depends on arbitrary functions. For example, the solution of the wave equation u

x x

for constant k and ρ can be of the form u(x; t)= f(x¡ct)+ g(x+ct), where c=

, and f

and g are arbitrary smooth functions.

• The geometry of the solution to an ODE is a curve, while for a PDE it is a surface in space

or g enerally a hypersurface. For example, the graph of the solution u(x; t) of a wave equation

in the space (x; t; u) is a surface and not a curve.

• The solution of the wave equation can be interpreted from a physical perspective as a trav-

eling wave moving to the right as f(x ¡ c t) and left as g(x +c t) with velocity c =

In contrast, the ordinary diﬀerential equations for a c oupled mass-spring system do not

reﬂect such a traveling speed. For example, consider a system of three masses connected to

three springs in series, where the initial displacement of the ﬁrst mass is x

(0)=1 and the

other two are at rest without stretch. Assuming zero initial velocity for all three masses

(x_

(0)= x_

(0) = x_

(0)), the initial kinetic energy of the system is zero, and the system evolves

its dynamics based on its initial potential energy U(0) =

(0)j

The following ﬁgure shows the motion of all three masses. As we observe, the motion of m

aﬀects m

immediately, i.e., the pro pagation speed of the disturbance in the system is inﬁnite.

1.1 From ODEs to PDEs 5

0 2 4 6 8 10

-1.5

-1

-0.5

0.5

1.5

Exercise 1.3. For each of the following partial diﬀerential equation, verify that the given solution satisﬁes the

eq uation

a) y@

u ¡x@

u = 0, u = f (x

+ y

)

b) x@

u ¡ y@

u = 0, u = f(x; y)

c) @

u ¡@

u = ¡u, u = f (x + y) e

¡x

Exercise 1.4. For a multivariable function u = u(x; y), equations that involve only t he partial derivatives of o ne

variable are gene rally referred to as defective equations. One example is the equation u

+ u = 0, where u = u(x;

y). This equation can be solved by the methods outlined in ODEs, resulting in the solution u(x; y)= f(y)e

¡x

where f is an arbitrary function. The function f serves as a constant parame ter for integrating the ODE

+ u = 0:

The following ﬁgure show the geometry of th e above argument

u = f (c

) e

¡x

u = f(c

) e

¡x

u = f (c

) e

¡x

y = c

+ u = 0

y = c

+ u = 0

y = c

+ u = 0

Using the appro priate method, solve the following eq uations, where u is a two-variable function u = u(x; y):

a) u

+ c

u = 0, c > 0 a constant.

b) u

+ u = y

c) @

u + u = u

d) u

+ u

= 0

e) x

+ x u

+ u = y.

1.1.3 From Newton cooling law to heat equation

Remember from ODEs that the temperature of a quantity located in an ambient space of temper-

ature T

follows the Newton's cooling rule:

= α(T

¡ T );

6 Introduction

where α>0 is a proportionality factor that depends on the conductivity of the material, and T =T (t)

is the temperature function of the material as a function of time t. The function T

¡ T is the

temperature gradient of the material compared to the temperature of the ambient space, T

. The

solution of the above ODE is an exponential function:

T (t) = T

+ (T

¡ T

) e

¡αt

;

where T

= T (0) is the initial temperature of the material.

Now, consider a conductive ro d that is insulated from the surrounding environment, with an

initial temperature proﬁle given by the function f(x). The temperature at any point x of the rod

at time t is denoted by u(x; t). Consider the interval between points a and b on this rod. If q(x; t)

is the density of the thermal energy, then the total thermal energy in this interval is:

[a;b]

(t) =

cρq(x; t) dx;

where ρ is the mass density and c is the speciﬁc h eat capacity of the material. The change in energy

from t to t + δt is:

[a;b]

(t + δt) ¡ Q

[a;b]

(t) =

cρ [q(x; t + δ t) ¡ q(x; t)] dx;

and by thermodynamics:

q(x; t + δt)¡ q(x; t)≈u(x; t + δt)¡u(x; t):

In the limiting case as δt!0, we obtain:

[a;b]

(t) =

cρu

(x; t) dx:

On the other hand, if no thermal energy is produced or lost in the interval, the only factor that

contributes to the change of thermal energy in this interval is the energy escaping through the

endpoints x=a and x =b. The Fourier law states that the heat gradient is proportional to u

, that is:

[a;b]

(t) = αu

(b; t) ¡ αu

(a; t);

where α is the conductivity factor of the material. Assuming α is a constant, the above relation can

be written as:

αu

(b; t) ¡ αu

(a; t) =

αu

x x

(x; t) dx;

and ﬁnally:

[cρu

(x; t) + αu

x x

(x; t)] dx = 0:

This integral equation holds for arbitrary interval [a; b], and if the integrand is continuous, we

conclude the following partial for u(x; t):

(x; t) = ku

x x

(x; t);

for k =

ρc

. This i s the one-dimensional heat equation, which governs the evolution of temperature in

a homogeneous material. The equation states that the time rate of change of the thermal energy in

any interval is prop ortional to the rate of change of temperature and to the second spatial derivative

of temperature.

1.1 From ODEs to PDEs 7

The heat equation has important applications in various ﬁelds such as physics, engineering,

and ﬁnance. In physics, it is used to model the diﬀusion of heat in a medium, the propagation of

electromagnetic waves, and the behavior of quantum systems. In engineering, it is used to design

and analyze heat transfer systems, and to model the behavior of materials subjected to thermal

stress. In ﬁnance, it is used to model the pricing of ﬁ nancial derivatives, such as options and futures,

where the underlying asset price follows a random diﬀusion process.

Exercise 1.5. Consider the heat equation u

= ku

where k > 0 is a constant.

a) Verify that the following function called the fundamental solution satisﬁes t he equation for t > 0

Φ(x; t) =

4πkt

p e

4kt

b) Assume that f (x) is a smooth funct ion. Verify the following relation

lim

t! 0

¡1

Φ(x; t) f (x) dx = f (0):

as long as the integral exists. The above relation implies that Φ(x; t) behaves like a Dirac delta function

when t !0.

Exercise 1.6. By Newton's cooling law, the h eat ﬂows from hot spots to cold s pot.

a) Explain how this fact is embedded in the derivation of the heat equation. In particular, explain why k

must be a positive constant.

b) If k < 0, the equation is called a reve rse heat equation. Explain why in this case heat ﬂows from cold spots

to hot spots.

Exercise 1.7. Consider a conductive rod of length L which is insulated along its length and are insulated at

the end points x = 0 and x = L. The i nsulation can be expressed mathematically as u

(0; t) = 0 and u

(L; t) = 0,

where u

denotes the temperature gradient. If so, show that the thermal energy in the system is conserved, that

is,

d Q

= 0, where Q is

[0;L]

(t) =

u(x; t) dt:

1.1.4 From population d ynam ics to migration

The exponential equation, formulated in the eighteenth century, was the ﬁrst mathematical model

for population dynamics. It describes the growth of a living species in a certain region through the

equation

= αP , where P =P (t) denotes the population and α is the growth rate. However, this

model has proven to b e inaccurate, and more sophisticated models, such as the logistic model, have

been suggested.

In this cont ext, we consider a possible extension of population dynamics for a living species that

migrates along a long strip. Speciﬁcally, we examine a living zone in the shape of a strip populated

with animals that have a fertility rate of α. We assume that these animals migrate along the strip

with a constant velocity of c. Our o bjective is to determine the population in an arbitrary segment

of the strip at any given time.

To achieve this, we need to deﬁne the population density function ρ(x;t) such that the population

in a segment [a; b] at time t is equal to

[a;b]

(t) =

ρ(x; t) dx:

Here, we present the derivation through the following exercise.

8 Introduction

Exercise 1.8. Based on the above settings

a) Suppose th e oﬀspring rat e is zero, α = 0. The rate of cha nge of P

[a;b]

in time is

[a;b]

ρ(x; t) dx =

(x; t) dx:

By c onservation of population, this rate of change is equal to the rate of moving in this segment minus

the rate of animals leaving this segment. Show that

[a;b]

= cρ(a; t) ¡cρ(b; t);

and conclude

[ρ

(x; t) + cρ

(x; t)] dx = 0:

The above relation holds for any segment [a; b], and thus we conclude the f ollowing diﬀerential equation

(x; t) + cρ

(x; t) = 0:

b) Verify that the function ρ = f(x ¡ ct) solves th e above PDE for ar bitrary smooth function f . Conclude

that if ρ(x; 0) = ρ

¡x

, then ρ(x; t) = ρ

¡(x¡ct)

. If c = 1, ﬁnd the population in the segment [¡1; 1] at

time t = 3.

c) Now, let us assume that α > 0 and

dρ

= αρ. Verify that the PDE in this case is

+ cρ

= αρ:

Veri fy that t he functi on ρ = f (x ¡ct) e

αt

sol ve s the equation. For ρ(x; 0) = 100e

¡x

, c = 1 and α = 0.3,

draw the solu tion ρ(x; t) at t = 0; 1; 2; 3.

1.2 Diﬀerential operators in higher dimensions

In order to extend the study of partial diﬀerential equations to higher dimensional spaces, we require

the use of appropriate diﬀerential operators. Among the most crucial of these operators are the

gradient, divergence, and Laplacian. These operators enable us to express and manipulate partial

derivatives in these spaces, which is essential for solving many important problems in ﬁelds such as

physics, engineering, and mathematics.

Before we begin, let's establish some general notation used in this book. The set of real numbers

is denoted by R, and t he set of n-tuples by R

. In R

, a point is written a s (x; y) ; in R

, as (x; y; z);

and in R

, as x=(x

; :::; x

). The standard unit vectors in R

are denoted by

; :::; e^

;

where e^

is a vector with zeros in all coordinates except the j-th coordinate, which is 1. Thus, a

point x can be considered as a vector x = x

+ ··· + x

. The set of all vectors in R

forms a

vector space with scalar multiplication given by

λ (x

; :::; x

) = (λx

; :::; λx

); λ 2 R

and vector addition deﬁned as

; :::; x

) + (y

; :::; y

) = (x

+ y

; :::; x

+ y

The magnitude of a vector x = (x

; :::; x

) is denoted as jxj and deﬁned as

jxj =

j=1

1.2 Differential oper ators in higher dimensions 9

The standard dot product between two vectors x and y is denoted as x · y and deﬁned as:

x · y =

j=1

Alternatively, we may use the more general notation h;i for the dot product, i.e., hx; yi=x · y. The

distance between two points x and y in R

is deﬁned as

jx¡ yj = hx¡ y; x¡ y i

= (x₁¡ y₁)² +···+(xₙ ¡ yₙ)²

For a point a 2 Rⁿ and r>0, Bᵣ(a) denotes the r-ball centered at a. The unit ball with the center

at the origin is denoted by B. A set Ω⊂Rⁿ is open if for any a 2Ω, there exists an r > 0 such

that B

(a) ⊂ Ω. Equivalently, a set is open if that every point in the set has a neighborhood

that is contained entirely within the set. The boundary of a set Ω is denoted by bnd(Ω), and

its closure by Ω

. For example, the set Ω = f(x; y) ; r

< jx ¡ yj < r

g is open, with the closure

Ω

= f(x; y); r

≤ jx ¡ yj ≤ r

g. The set Ω

= f(x; y); r

< jx ¡ yj ≤ r

g is neither open nor closed.

We use the notation u

; u

for the partial derivatives of a two-variable function u=u(x; y):

: =

; u

On occasion, we also use the notation @

u and @

u for the partial derivatives. For a function u of

n variables, u=u(x₁; :::; xₙ), the partial derivatives are denoted by @ⱼu for

. The second-order

partial derivatives are similarly denoted by @

u =

1.2.1 Diﬀerential operat o rs in Cartesian coordinate

The diﬀerential operator r, called nabla, in R

is deﬁned as:

r := e^

+ ··· + e^

Let u = u(x

; :::; x

) be a smooth scalar function deﬁned on an open set Ω ⊂ R

. The gradient of u,

denoted by ru, is deﬁned as:

ru = @

u e^

+ ··· + @

u e^

For a ﬁxed point x

2 Ω, ru(x

) is a vector, and ru(x) for x 2 Ω deﬁnes a vector ﬁeld over Ω

which assigns a vector to each point in Ω. This vector ﬁeld is called the gradient vector ﬁeld. The

gradient of u, is a vector ﬁeld that points in the direction of maximum increase of u at each point in

Ω. The following ﬁgures illustrate the graph of the function u=xy exp(¡x

¡ y

) and the grad ient

trajectories which are the solutions of the following system:





= ru(x; y);

or equivalently

= (1 ¡ 2x

)ye

¡x

¡y

= (1 ¡ 2y

)xe

¡x

¡y

10 Introduction

The direction of the gradient ﬂow lines in heat ﬂow or ﬂuid dynamics is opposite to that of heat

or ﬂuid ﬂow. This is because heat and ﬂuid tend to move from regions with higher temperatures or

pressures to those with lower temperatures or pressures, whereas the gradient vector ﬁeld points in

the direction of maximum increase of the scalar f unction at each p oint.

In the context of a two-dimensional elastic membrane, the function u(x; y) represents the

vertical displacement of the membrane at point (x; y). The tension T (x; y) on the boundary of the

membrane, bnd(Ω), is deﬁned as the product of a constant τ and the gradient of u(x; y), ru(x; y).

This means that the tension at a point on the boundary is proportional to the rate of change of the

displacement at that point.

The vertical Hook's force exerted at a point (x; y) in the interior of the membrane is deﬁned as

the limit of the average tension along the boundary of a shrinking region Ω around the point, as the

size of Ω approaches zero. The force is given by the line integral of the tension over the boundary

bnd(Ω), with respect to the outward unit normal vector ν, divided by the area of Ω

F ( x; y) = lim

jΩj!0

A(Ω)

bnd( Ω)

τ ru · νdS:

In other words, the Hook's force measures the rate of change of the tension with respect to the area

of the membrane. This force is important in studying the behavior of elastic membranes, such as

in the design of structures that use membranes as a supporting element.

Another important fact that we will use extensively in the next chapter of this book is the time

or mass derivative of a scalar function along a smooth curve. For a smooth parametrized curve

γ(t)=(x(t); y(t)), the derivative of u(x; y) along γ is deﬁned as:

(γ(t)) = ru(γ(t)) · γ

(t);

where γ

(t) is the tangent vector to γ at time t. This expression represents the instantaneous rate of

change of u(x; y) at the point γ(t) along the curve γ(t), and it takes into account both the direction

and magnitude of the tangent vector γ

(t). In other words, it tells us how quickly u(x; y) is changing

as we move along the curve γ(t).

In the context of ﬂuid mechanics, if γ(t) denotes the trajectory of a ﬂuid particle in the xy-

plane, then γ

(t) represents its velocity vector, and

represents the instantaneous rate of change

of u(x; y) at each point on the trajectory of the particle as it moves through the ﬂuid.

1.1. It appears that the natural world operates more justly than our modern economic system. In the context of money,

wealth tends to ﬂow from poor individuals to those who already possess more, following the gradient of wealth inequality.

However, the laws of physics dictate the opposite behavior: heat, pressure, mass, and other quantities ﬂow from regions of high

density to those of lower density.

1.2 Differential oper ators in higher dimensions 11

Exercise 1.9. Con sider the temperat ure distribution in the xy-plane given by T (x; y)=10 e xp(x

¡y

). Suppose

a runner is moving along the curve γ(t)=(cos(t); sin(t)), and weari ng a hand clock marked f rom 0 to 60. What is

the angular velocity of the hand of the runner's watch as they move along the curve? Find the angular velocity

of the hand clock if the runner runs given by the following funct ion γ(t) = (cos(!

t); sin(!

t)) for some !

> 0.

The divergence of a vector ﬁeld F =(f ₁; :::; fₙ) on a domain Ω⊂Rⁿ is a scalar function denoted

by r · F , deﬁned by

r · F =

j=1

Geometrically, the divergence of a vector ﬁeld measures the extent to which the ﬁeld “ﬂows out”

from a given point in space. If r · F (x

)> 0, then the point x

b ehaves like a source for the ﬁeld; that

is, the ﬁeld ﬂows outward from x

. If r · F (x

) > 0, then x

b ehaves like a sink, and the ﬁeld ﬂows

inward toward x

. If r · F (x

) =0, then x

is a neutral point for the ﬁeld. The concept of divergence

is fundamental in ﬂuid dynamics, where it is used to study the behavior of ﬂuids in motion.

For example, the vector ﬁeld F = ru, for u = xy exp(¡x

¡ y

), is equal to

r · F (x; y) = 4xy(x

+ y

¡ 3) exp(¡x

¡ y

For a vector ﬁeld x 7! F (x) = (f

(x); :::; f

(x)) for x 2 Ω ⊂ R

, the divergence denoted by r · F is

deﬁned as the scalar function

r · F =

j=1

The divergence of a vector ﬁeld measures how much the vector ﬁeld “ﬂows out” from a given point

in space. If r · F (x

) is positive, t hen x

b ehaves like a source for the given ﬁeld. If r · F (x

) is

negative, it behaves like a sink and if r · F (x

) is zero, the point is a neutral point for the ﬁeld.

The following ﬁgure depicts the vector ﬁeld F = ru, where u = xy exp(¡x

¡ y

), and the

corresponding contours of r · F .

-2 -1 0 1 2

-2

-1.5

-1

-0.5

0.5

1.5

One of the important results related to the divergence of a vector ﬁeld is the Gauss or divergence

theorem. The theorem relates the integral of the divergence of a vector ﬁeld over a volume to the

ﬂux o f the vector ﬁeld through the surface e nclosing the volume. It is a fundamental theorem in

vector calculus and has many applications in physics and engineering.

Theorem 1.1. (Gauss) Assume that F is a smooth vector ﬁeld in an open bounded set Ω ⊂ R

The following equality holds:

Ω

r · FdV =

bnd( Ω)

F ·νdS;

12 Introduction

where dV is the volume element, dS is the surface element, ν is the unit normal vector to the surface,

and the dot product F · ν is the ﬂux of the vector ﬁeld through the surface.

Exercise 1.10. Let B

denote the ball o f radius R in R

and S

n¡1

the surface of the ball. Use divergence

theorem to show that

V (B

) =

A(S

n¡1

)

where V (B

) and A(S

n¡1

) are the volume and surface area of B

and S

n¡ 1

respectively. Verify your formula

for n = 2; 3. (H int: consider the function u = x

+ ···+ x

Remark 1.1. The divergence of a vector ﬁeld at a point x can be deﬁned using the divergence

theorem. Speciﬁcally, let F be a smooth vector ﬁeld in an open set containing x, and let Ω be a

small, op en, three-dimensional ball centered at x. Then, the divergence of F at x is given by:

r · f(x) = lim

jΩj!0

vol(Ω)

bnd( Ω)

F ·νdS;

where vol(Ω) is the volume of Ω, and the limit is taken as the size of Ω approaches zero. In other

words, the divergence of F at x is equal to the ﬂux of F across the boundary of a small ball centered

at x, divided by the volume of the ball.

Consider a homogeneous ﬂuid with a constant density distributed in the regionx

+ y

< 1, which

is shown in the ﬁgure below:

Suppose the ﬂuid moves a ccording to the velocity vector ﬁeld V

= (¡x; ¡y). As expected from

the shape of the vector ﬁeld, the ﬂuid will concentrate at the origin at later times, and the density

of the ﬂuid will increase around the origin. In fact, the ﬂow lines of this vector ﬁeld are the solutions

of the following system of ODEs

= ¡x

= ¡y

;

which is solved for γ(t) = (x

¡t

; y

¡t

). The top ﬁgure in the following image depicts this scenario.

If the velocity ﬁeld changes to V

= (x; y), the ﬂuid will spread out, according to the system

= x

= y

;

and the density will decrease around the origin, as shown in the bottom ﬁgure.

1.2 Differential oper ators in higher dimensions 13

Note that the behavior of the ﬂuid is determined by the divergence of the velocity ﬁeld. In the

ﬁrst scenario, the divergence of V

is negative, indicating that the ﬂuid is ﬂ owing inward toward

the origin, while in the second scenario, the divergence of V

is positive, indicating that the ﬂuid

is ﬂowing outward from the origin. This observation is consistent with the divergence theorem

presented earlier, which shows the relationship between the ﬂow of a vector ﬁeld and its divergence.

Exercise 1.11. Let the initial density of a ﬂuid ﬂow is given by the function

f(x; y) =

(

xy exp(¡x

¡ y

) x

+ y

< 1

0 otherwise

Find the ﬂow lines of the ﬂuid if the velocity vect or ﬁeld is given by V = (¡y; x). Round the following code in

Matlab and explain the change of density function in time.

[thet,r]=meshgrid(-pi:pi/1000:pi,0:0.005:1);

[x,y]=pol2cart(thet,r);

xt=@(t) x*cos(t)+y*sin(t);

yt=@(t) -x*sin(t)+y*cos(t);

zt=@(t) (xt(t).*yt(t).*exp(-x.^2-y.^2)).*(x.^2+y.^2<1)+0.*(x.^2+y.^2>=1);

for i=0:1:2

subplot(2,2,i+1)

surf(x,y,zt(i*pi/4)); shading interp; view(2); axis equal tight; grid off; colormap jet

clim([-0.25,0.25]); colorbar;

title(sprintf('density at $t=%.2f$',i*pi/4),'interpreter','latex','fontsize',10);

end

Exercise 1.12. A gradient vector ﬁeld of the form F = ¡ru for a scalar ﬁled is also called conservative ﬁeld s.

The reason is due to the following fact. Consider mass particle m moving along an arbitrary path γ(t) in this

ﬁeld. The total energy of the mass is

E =

m jv j

+ u(x):

Use the Newton's second law and show that the derivative of E along γ(t) is zero. Therefore, a conservative

force ﬁeld conserves the total energy of a mass particle.

14 Introduction

Exercise 1.13. The line integral of a vector ﬁeld F alo ng an arbitrary smooth path γ(t) for t 2(a; b) is deﬁned

as follows

F :=

F (γ(t)) · γ

(t) dt:

Show that if F = ¡ru, then the integral is independent of γ and is equal to u(γ(b)) ¡ u(γ(a)). This integral

shows the amount o f work or kinet ic energy spend to move a particle from point γ(a) to γ(b) and as it is seen,

this amo unt is equal to the diﬀerence between the po tential at these two po ints.

Exercise 1.14. Show that the line integral of a vect or ﬁeld is independent of the parametrizatio n, that is, if γ

(c; d) !C the n

F (γ(t)) · γ

(t) dt =

F (γ

(t)) · γ

(t) dt:

Exercise 1.15. The diver ge nce of a two dimensional smooth vector ﬁled F = (f ; g) is

r·F = f

+ g

Let R be the rectangle R := [¡a; a] ×[¡b; b]. We prove the relation

r·F (0 ; 0) = lim

a;b! 0

lim

a;b! 0

A(R)

bnd(R)

F ·ν dl

a) Show that the net ﬂux passing through the boundary of R is

¡b

ff (a; y) ¡ f (¡a; y)gdy +

¡a

fg(x; b) ¡ g(x; ¡b)gdx:

b) Use the r elations

f(a; y) ¡ f(¡a; y) =

¡a

(x; y)dx; g(x; b) ¡ g(x; ¡b) =

¡b

(x; y)dy;

and show the relation

¡b

ff (a; y) ¡ f (¡a; y)gdy +

¡a

fg(x; b) ¡ g(x; ¡b)gdx = 4abff

(ξ

; η

) + g

(ξ

; η

for some ¡a < ξ

; ξ

< a, ¡ b < η

; η

< b a nd conclu de

lim

a;b!0

A(R)

bnd(R)

F ·νdl = f

(0; 0) + g

(0; 0) = r·F (0; 0):

c) Use the same argument a nd show that for the cube C := [¡a; a] × [¡b; b] × [¡c; c] a nd the smooth ﬁeld

F = (f ; g; h), the followin g relation holds

lim

a;b;c!0

vol(C)



bnd(C)

F ·νdS = f

+ g

+ h



(0;0;0)

Exercise 1.16. Let B

⊂ R

be a ball of radius a in n-dimensional space. Consider the vector ﬁeld F (x) = x,

and use t he divergence theorem to conclude that the volum e of B

to its surface are a is equal to

Exercise 1.17. Let Ω ⊂ R

be a domain with sm ooth boundary su rface bnd(Ω). Use the divergence theorem

and show the fo llowing relations

a) Let F be a smooth vector ﬁeld u a smoot h scalar functi on, then

div (uF ) = u div (F ) + F ·ru;

where div sta nds f or r·.

b) Use the above result and show the following relation

ZZZ

Ω

u div (F ) dV =



bnd(Ω)

u F ·ν dS ¡

ZZZ

Ω

F ·ru dV :

c) If φ; are smooth functions de ﬁned on Ω, then

ZZZ

Ω

∆φdV =



bnd(Ω)

φdS ¡

ZZZ

Ω

rφ ·r dV ;

1.2 Differential oper ators in higher dimensions 15

where @

φ stands for rφ ·ν.

d) Use the above result and show the following relation called the Green's formu la:

ZZZ

Ω

(φ∆ ¡ ∆φ) dV =



bnd(Ω)

[φ @

¡ @

φ] dS:

The Laplacian operator ∆ := r · r is a diﬀerential operator that is commonly used in mathe-

matics and physics. As it is evident from its deﬁnition, the Laplacian of a scalar function u is equal

to the divergence of the gradient of u, and in the coordinate form is equal to

∆u =

j=1

for a smooth function u = u(x

; :::; x

The Laplacian ope rator has important applications in physics, such as in the study of the diﬀu-

sion e quation, wave equation, and S chrödinger equation. In these contexts, the Laplacian operator

describes the be havior of p hysical quantities such as temperature, pressur e, electric potential, and

wave functions. In the following discussion, we present an example of how the Laplacian operator

is applied in the ﬁeld of electrostatics.

Let's consider an electric charge Q located at the origin of a three-dimensional space. The

potential ﬁeld generated by this charge at any po int r = (x; y; z) =/ (0; 0; 0) is given by

Φ(r) =

4π"

jrj

;

where "

is the p ermittivity of the air and jrj is the distance of x to the origin

jrj = x

+ y

+ z

This potential ﬁeld generates an electric ﬁeld

E(r) := ¡rΦ(r) =

4π"

jrj

;

with a magnitude of jE(r)j =

4π"

jrj

and direction of r^ =

jrj

. Let B

b e the ball with radius

R centered at the origin. We can calculate the net amount of electric ﬁeld passing through the

boundary of this ball, denoted by S

, as follows:



E(r) · ν(r) dS ;

where ν = r^ for r =/ 0. Hence,



E(r) · ν(r) dS =



4π"

dS =

4π"



dS =

As observed, the net amount of electric ﬁeld passing through S

is independent of R and is equal

Now, let's consider the divergence of the electric ﬁeld E(r) for r =/ 0. By direct calculation, we

have

∆Φ = r · E(r) =



4π"

jrj





4π"

jrj





4π"

jrj



4π"



jrj

¡ 3x

jrj

¡ 3y

jrj

¡ 3z

jrj



= 0

16 Introduction

Therefore, ∆Φ = 0 everywhere except at r = 0.

Note that we can apply the divergence theorem to the integral of ∆Φ over the region between

two concentric balls B

and B

, which gives

ZZZ

∆Φ dV =

ZZZ

∆Φ d V +

ZZZ

¡B

∆Φ dV :

∆Φ = 0

By the result ∆Φ = 0 for r =/ 0, we conclude

ZZZ

∆ΦdV =

ZZZ

∆ΦdV ;

and by the Gauss theorem, we obtain



E · νdS =



E · νdS:

One can use the above results and solve the following problems.

Exercise 1.18. Let Ω be any open b ounded domain with the smooth boundary bnd(Ω), and let and Φ be the

potential ﬁeld of a pointwise Q-charge. Show that if Q is located insi de Ω, then

ZZZ

Ω

∆ΦdV =

;

and if Q is located outsi de Ω

, then

ZZZ

Ω

∆ΦdV = 0:

Exercise 1.19. Let two elec tric charges Q; ¡Q are loca ted inside an open bounded domain Ω. Show

ZZZ

Ω

∆ΦdV = 0:

Exercise 1.20. Suppos e there is an open bounded set Ω in R

, and it contains an electric charge density q. The

total electri c cha rge inside Ω can be expressed as

Q =

ZZZ

Ω

q(r) dV :

Show the following relation

ZZZ

Ω

∆ΦdV =

Hint : consider a diﬀerential elem ent qdV and le t E

(r) be the electric ﬁeld generated by qdV at point r.

1.2.2 Diﬀerential operat o rs in polar and spherical coordinates

Remember that the polar coordinate is deﬁned through the transformation

T (r; θ) = (x; y) = (r cos(θ); r sin(θ) );

1.2 Differential oper ators in higher dimensions 17

where r 2 [0;1) and θ 2[¡π; π]. In order to ﬁ nd the form of the operator nabla r in this coordinate,

wen need ﬁrst to determine the vectors r^; θ

, the unit vector along r and θ directions. For ﬁxed θ = θ

The transformation T (r; θ

) is a parameterized curve and then T

(r; θ

) is tangent to this curve.

The unit vector r

is then obtained as

r^ =

(r; θ

)

(r; θ)j

= (cos(θ

); sin(θ

)):

As it is observed, the unit vector r^ depends on the chosen angle θ

, and then

= r

(θ) = (cos(θ); sin(θ)):

Now, ﬁx r = r

, and determine θ

by the relation

; θ)

; θ)j

= (¡ sin(θ); cos(θ)):

As it is seen again, θ

= θ

(θ), is a function of θ and indepe ndent of r.

By the above calculations, it is evident that r^ · θ

= 0, and thus they are orthogonal. Transfor-

mations for which this property holds are said orthogonal transformations. Moreover, we have

@ r

@θ

= θ

;

@ θ

@θ

= ¡r^;

and

@ r^

@ θ

= 0. Note that

@ θ

@θ

= ¡r^ contributes to the centrifugal force experienced by a particle

moving in a circular path around the origin in the polar coordinate system. The vector

@ r^

@θ

= θ

, on

the other hand, points in the tangential direction to the curve traced out by a particle moving in a

circular path around the origin in the polar co ordinate system.

Now, let u(r; θ) be a given scalar function in polar coordinates. We can write ru in polar

coordinates as:

ru = α(r; θ) r^ + β(r; θ) θ

;

where we need to determine the appropriate functions α and β. To ﬁnd α, we use the orthogonality

condition, which yields:

α(r; θ) = ru · r^ = (u

; u

) · (cos(θ); sin(θ)) = u

cos(θ) + u

sin(θ):

Using the chain rule, we can also write:

= u

+ u

= u

cos(θ) + u

sin(θ);

and so we see that α(r; θ)=u

. Similarly, we can ﬁnd β as:

β(r; θ) = ru · θ

= ¡u

sin(θ) + u

cos(θ) =

Therefore, we obtain:

(r;θ)

u(r; θ) = u

;

and we can express the diﬀerential op erator r

(r;θ)

as:

(r;θ)

= r

;

where @

; @

denote partial derivatives with respect to r and θ, respectively.

18 Introduction

With the nabla operator in polar coordinates, we can determine the Laplacian ∆ := r · r, which

is crucial for our discussion of second-order partial diﬀerential equations. By u sing the Laplacian in

polar coordinates, we can solve a variety of problems involving circular or cylindrical domains. For

example, we can model the behavior of electric or magnetic ﬁelds inside a cylindrical conductor, or

the ﬂow of ﬂuids inside a circular pipe. The Laplacian also plays an important role in the study of

harmonic functions, which are solutions to the Laplace equation ∆u=0.

Based on the relation ∆:=r·r and the derived r

(r;θ)

, we can obtain the Laplacian in polar

coordinates as follows:

∆ = r · r =



r^ @





r^ @



Using the orthogonality relation r^ · θ

= 0, and the partial derivative relations:

r^ = @

= 0

r^ = θ

= ¡r^

;

we can simplify the Laplacian to the following form in polar coordinates:

∆

(r;θ)

= @

θθ

Exercise 1.21. Find th e form of nabla r in th e cylindr ical coordinate and then ﬁnd the Laplacian in this

coor dinat e syste m.

Let's now turn to the spherical co ordinate system. In this system, the transformation is deﬁned

through the following relations:

T (ρ; θ; φ) = (x; y; z) = (r sin(θ) cos(φ); r sin(θ) sin(φ); r cos(θ)):

Here, r 2[0; 1) is the distance from the origin, θ 2[0; π] is the angle with respect to the z-axis, a nd

φ2[¡π; π] is the angle with respect to the x-axis.

The unit vector r^ is deﬁned as

r^ =

(r; θ; φ)

(r; θ; φ)j

= (sin(θ) cos(φ); sin(θ) sin(φ); cos(θ)):

Similarly, we have the following expressions for θ

and φ

(r; θ; φ)

(r; θ; φ)j

= (cos(θ) cos(φ); cos(θ) sin(φ); ¡sin(θ));

(r; θ; φ)

(r; θ; φ)j

= (¡sin(φ); cos(φ); 0):

1.2 Differential oper ators in higher dimensions 19

It can be easily veriﬁed that r^; θ

, and φ

are mutually orthogonal, and thus the transformation T

deﬁnes an orthogonal transformation.

Given a scalar function u = u(r; θ; φ), we can determine its gradient using the following

expression:

ru = α(r; θ; φ) r^ + β(r; θ; φ) θ

+ γ(r; θ; φ) φ

;

where functions α; β, and γ need to be determined. To ﬁnd these functions, we can use the

orthogonality conditions:

α = ru · r^ = u

sin(θ) cos(φ) + u

sin(θ) sin(φ) + u

cos(θ) = u

;

β = ru · θ

= u

cos(θ) cos(φ) + u

cos(θ) sin(φ) ¡ u

sin(θ) =

;

γ = ru · φ

= ¡u

sin(φ) + u

cos(φ) =

r sin(θ)

Using these expressions, we can derive the form of r

(r;θ;φ)

= r

r sin(θ )

This can be used to determine the Laplacian in this coordinate system.

Exercise 1.22. To derive the Laplacian in the spherical coordinate:

a) Show the following relations f or the spherical tran sformation

r^ = 0; @

r^ = θ

; @

r^ = sin(θ) φ

= 0; @

= ¡r^; @

= cos (θ)φ

= 0; @

= ¡sin(θ) r^ ¡cos(θ) θ

b) Find the dive rgence of a vector ﬁled F = α(r; θ; φ)r^ + β(r; θ; φ) θ

+ γ(r; θ; φ)φ

in this coordinate syste m.

c) Use r

(r;θ ;φ)

and the result of the above problem, and derive ∆

(r;θ; φ)

which is

∆

(r;θ ;φ)

) +

sin (θ)

(sin(θ) @

) +

sin

(θ)

φφ

Exercise 1.23. In this exercise, we wi ll derive the form of the Laplacian for a general orthogonal curvilinear

coor dinat e syste m deﬁ ned by the transformation:

T (q

; q

) = (x( q

; q

); y(q

; q

); z(q

; q

));

where q

; q

, and q

ar e the coordinates of a point i n the curvili near system, and x; y, and z are their Cartesian

counterpart s.

a) Find the unit vectors q^

; q^

, and q^

of the coordin ate system.

b) Use the ortho gonality condition and show the following relation for the nabla operator r

)

;

where @

, and T

c) We use symbols γ

i; j

for the coordinates of the partial derivatives of

@ q^

as follows

@ q^

= γ

^ + γ

^ :

Show that γ

= 0.

d) Use the above notation and conclude

∆ := r·r=

i=1





i; j=1

jjT

i;j

e) For a given vector ﬁeld F = P q^

+ Q q^

+ R q^

, ﬁnd div (F ).

20 Introduction