Partial differential equations
This chapter is an introduction to PDE with physical examples that allow straightforward numerical solution with Mathemat-
ica. Methods of solution of PDEs that require more analytical work may be will be considered in subsequent chapters.
General facts about PDE
Partial differential equations (PDE) are equations for functions of several variables that contain partial derivatives. Typical
PDEs are Laplace equation
∆φ@x, y, …D 0
(D is the Laplace operator), Poisson equation (Laplace equation with a source)
∆φ@x, y, …D f@x, y, …D,
wave equation
∂2t φ@t, x, y, …D − c2 ∆φ@t, x, y, …D 0,
heat conduction / diffusion equation
∂t φ@t, x, y, …D − κ∆φ@t, x, y, …D 0,
Schrödinger equation
∂t φ@t, x, y, …D + Ha∆ + bf@x, y, …DL φ@t, x, y, …D 0,
etc. There are systems of PDEs and nonlinear PDEs.
Solution of PDEs have more freedom than those of ODEs because integration "constants" are in fact functions. For instance,
the general solution of the second-order PDE
∂x,y f@x, yD 0
is
f@x, yD = F@xD + G@yD,
where F@xD and G@yD are arbitrary functions. The solution of the first-order PDE
∂t f@t, xD − v ∂x f@t, xD 0
is
f@t, xD = g@x − vtD
that describes a front of arbitrary shape moving in the positive direction if v > 0.
General analytical solutions of PDEs are available only in the simplest cases, and because of this freedom, they do not yet
solve the problem. The actual form of the solution is defined by the symmetry of the problem (if it exists) and boundary
conditions. If one of the variables is time, one usually speaks of initial conditions set at the initial time and of the boundary
conditions for spatial variables.
If there are initial conditions but no final conditions, the problem is evolutionary and it can be solved numerically starting
from the initial conditions and increasing time. The most efficient method to do it is the so-called "method of lines" used by
Mathematica. First, the problem is discretized in spacial variables and spatial derivatives are approximated by differences.
This reduces the PDE to a system of ODEs in time. Then the resulting system of ODEs is solved by one of high-performance
ODE solvers. In Mathematica, PDEs, as well as ODEs, are solved by NDSolve.
,2 Mathematical_physics-13-Partial_differential_equations.nb
However, currently Mathematica can only solve problems with a rectangular spatial region. For a typical PDE, in these cases
one also can find an analytical solution.
Mathematically, there is no difference between the time and other variables. If for a particular spacial variable a boundary
condition is set only at one end of the interval, one can treat this variable as time and the problem is evolutionary. Mathemat-
ica recognizes this situation and finds the solution. If boundary conditions are set at all ends of the interval (or infinity)
NDSolve does not find the solution and other methods have to be used. Most time-independent problems are like that.
Partial differential equations in physics
In physics, PDEs describe continua such as fluids, elastic solids, temperature and concentration distributions, electromag-
netic fields, and quantum-mechanical probabilities. Below we review these equations.
ü Continuity equation
Any conserved quantity described by the density r and the corresponding current j satisfies the continuity equation
∂ρ
+ divj 0. (1)
∂t
This is the continuity equation in the differential form. If „V is an infinitesimal volume, ∑t r„V is the rate of change of the
amount (corresponding to r) within this volume, whereas divj „V is the flux out of this volume. If the flux is positive, the
amount decreases. The integral form of the continuity equation can be obtained by integrating over a finite volume V
∂
‡ ρ V − ‡ divj V.
∂t
At the right, the integral over the volume can be expressed via a surface integral, the flux of j, according to the Gauss
theorem
∂
‡ ρ V − ‡ j ⋅ S.
∂t
As there are many types of densities in physics (mass density, energy density, charge density, etc.) there are also many types
of continuity equations.
ü Fluid dynamics
Motion of a fluid satisfies the continuity equation (1), there r is the mass density and the mass current j density is given by
j = ρv.
Liquids, unlike gases, are practically uncompressible, r = const, thus the continuity equation simplifies to
divv 0.
Dynamics of a non-viscous fluid is described by the Euler equation
∂v ∇P
+ Hv ∇L v − + g,
∂t ρ
where P is pressure and g is the gravity acceleration. Minus in front of the pressure gradient shows that the liquid accelerates
in the direction of a smaller pressure. The kinematic (streaming) term Hv “L v accounts for the change of the velocity at a
given point as a result of the media motion and makes the Euler equation nonlinear.
Flow of a viscous fluid obeys the Navier-Stokes equation. For an incompressible fluid it has the form
, Mathematical_physics-13-Partial_differential_equations.nb 3
∂v ∇P η
+ Hv ∇L v − +g+ ∆v,
∂t ρ ρ
where h is the viscosity coefficient.
Euler equation has a limited applicability because its solutions (for instance, for an object in a flow) have discontinuities of
the velocity in the stream. In these regions the viscosity becomes important and at high enough speeds the flow becomes
turbulent. Still, Euler equation generates Bernoilli equation that explains the lifting force on airplanes as well as the mecha-
nism of sound creation by voice bands.
Navier-Stokes equation is more realistic and it can be applied, for instance, to the problem of a drag force acting on a object
moving in a liquid. For small boby speeds u the quadratic terms can be dropped, equation becomes linear and can be solved
analytically for simple shapes such as sphere. The drag forces on a solid sphere and air bubble of radius R have the form
6 πηRu, solid sphere
F=µ
4 πηRu, air bubble.
The difference in the coefficient is due to different bounary conditions: The liquid sticks to the solid sphere but glides along
the surface of the air bubble that reduces the resistance.
ü Waves
Waves exist in compressible media such as fluids and solids. The main prerequisite for the waves' existence is a restoring
force. In fluids, there is only a restoring force with respect to change of the volume, thus the waves are waves of compression
/ expansion. In these waves the media displacement is collinear with the direction of the wave propagation, so that these
waves are called longitudinal waves. In solids, there is additionally a restoring force with respect to change of the shape of a
media element. The corresponding additional waves are transverse waves. Below we will consider small-amplitude longitudi-
nal waves in fluids.
In small-amplitude waves pressure and density only slightly deviate from their equilibrium values, so that one can define
P = P0 + δP, ρ = ρ0 + δρ
and linearize equations in small dP and dr, neglecting quadratic velocity terms because the velocity is small. The continuity
equation becomes
∂
δρ + ρ0 divv 0, (2)
∂t
whereas in the absence of viscosity from the Euler equation one obtains
∂v ∇ δP
− . (3)
∂t ρ0
The small changes of pressure and density are related by the adiabatic compressibility,
∂P
δP = δρ.
∂ρ S
One can obtain a single equation describing waves by differentiating the continuity equation over time and then using the
Euler equation. This yields the wave equation
∂P
∂2t δρ − c2 ∆δρ, c≡ ,
∂ρ S
where c is the speed of sound. Note that the speed of sound (that can be large) has no relation to the velocity of the media
(that is small). Since dP ~ dr, it satisfies the same equation,
, 4 Mathematical_physics-13-Partial_differential_equations.nb
∂2t δP − c2 ∆δP. (4)
Let us obtain the equation for the velocity. Taking rotor of the linearized Euler equation, one obtains
∂
∇ v 0,
∂t
thus velocity is a potential field, — ä v = 0, and it can be searched in the form
v = − ∇ φ.
Inserting it into the linearized Euler equation (3) and integrating it by removing the gradient one obtains
∑f dP c2
ã = dr. (5)
∑t r0 r0
Differentiating this over time and using the continuity equation, one obtains the same wave equation for the velocity potential,
∂2t φ − c2 ∆φ.
Let us now seek for the solution of the wave equation in the form of a plane wave
φ = ACos@ωt − k ⋅ r + φ0 D,
where A is the amplitude, w is the frequency and k is the wave vector that defines the direction of the wave propagation.
After substitution, differentiation, and cancellation one obtains the dispersion relation
ω2 − c2 k2 0.
The period T and wave length l are related to w and k by
2π 2π
ω= , k= .
T λ
One can see that w can be different while k adjusts to w as ck = w.
Let us now find the media velocity in a plane wave:
v = − ∇ φ = − A k Sin@ωt − k ⋅ r + φ0 D.
One can see that v is collinear to k, this is why such waves are called longitudinal waves.
Since the wave equation is linear, a superposition of plane waves moving in different directions with different frequencies
and wave lengths is also its solution. This illustrates the fact that there is much freedom in the solution of PDEs. The selec-
tion of the actual solution is done by the initial and boundary conditions.
In the important case of monochromatic waves whose time dependence is describe by a single frequency,
φ@r, tD = φ@rD Sin@ωtD,
or similar with Sin@wtD or Exp@ÂwtD, the wave equation simplifies to the time-independent stationary wave equation
ω
∆φ + k2 φ 0, k= .
c
In one dimension it becomes an ODE.
ü Thermal conductivity and diffusion
Temperature in an isolated body tends to equilibrate with time because the energy flows from warmer regions to colder
regions. Since the energy is conserved, its flow obeys the continuity equation
This chapter is an introduction to PDE with physical examples that allow straightforward numerical solution with Mathemat-
ica. Methods of solution of PDEs that require more analytical work may be will be considered in subsequent chapters.
General facts about PDE
Partial differential equations (PDE) are equations for functions of several variables that contain partial derivatives. Typical
PDEs are Laplace equation
∆φ@x, y, …D 0
(D is the Laplace operator), Poisson equation (Laplace equation with a source)
∆φ@x, y, …D f@x, y, …D,
wave equation
∂2t φ@t, x, y, …D − c2 ∆φ@t, x, y, …D 0,
heat conduction / diffusion equation
∂t φ@t, x, y, …D − κ∆φ@t, x, y, …D 0,
Schrödinger equation
∂t φ@t, x, y, …D + Ha∆ + bf@x, y, …DL φ@t, x, y, …D 0,
etc. There are systems of PDEs and nonlinear PDEs.
Solution of PDEs have more freedom than those of ODEs because integration "constants" are in fact functions. For instance,
the general solution of the second-order PDE
∂x,y f@x, yD 0
is
f@x, yD = F@xD + G@yD,
where F@xD and G@yD are arbitrary functions. The solution of the first-order PDE
∂t f@t, xD − v ∂x f@t, xD 0
is
f@t, xD = g@x − vtD
that describes a front of arbitrary shape moving in the positive direction if v > 0.
General analytical solutions of PDEs are available only in the simplest cases, and because of this freedom, they do not yet
solve the problem. The actual form of the solution is defined by the symmetry of the problem (if it exists) and boundary
conditions. If one of the variables is time, one usually speaks of initial conditions set at the initial time and of the boundary
conditions for spatial variables.
If there are initial conditions but no final conditions, the problem is evolutionary and it can be solved numerically starting
from the initial conditions and increasing time. The most efficient method to do it is the so-called "method of lines" used by
Mathematica. First, the problem is discretized in spacial variables and spatial derivatives are approximated by differences.
This reduces the PDE to a system of ODEs in time. Then the resulting system of ODEs is solved by one of high-performance
ODE solvers. In Mathematica, PDEs, as well as ODEs, are solved by NDSolve.
,2 Mathematical_physics-13-Partial_differential_equations.nb
However, currently Mathematica can only solve problems with a rectangular spatial region. For a typical PDE, in these cases
one also can find an analytical solution.
Mathematically, there is no difference between the time and other variables. If for a particular spacial variable a boundary
condition is set only at one end of the interval, one can treat this variable as time and the problem is evolutionary. Mathemat-
ica recognizes this situation and finds the solution. If boundary conditions are set at all ends of the interval (or infinity)
NDSolve does not find the solution and other methods have to be used. Most time-independent problems are like that.
Partial differential equations in physics
In physics, PDEs describe continua such as fluids, elastic solids, temperature and concentration distributions, electromag-
netic fields, and quantum-mechanical probabilities. Below we review these equations.
ü Continuity equation
Any conserved quantity described by the density r and the corresponding current j satisfies the continuity equation
∂ρ
+ divj 0. (1)
∂t
This is the continuity equation in the differential form. If „V is an infinitesimal volume, ∑t r„V is the rate of change of the
amount (corresponding to r) within this volume, whereas divj „V is the flux out of this volume. If the flux is positive, the
amount decreases. The integral form of the continuity equation can be obtained by integrating over a finite volume V
∂
‡ ρ V − ‡ divj V.
∂t
At the right, the integral over the volume can be expressed via a surface integral, the flux of j, according to the Gauss
theorem
∂
‡ ρ V − ‡ j ⋅ S.
∂t
As there are many types of densities in physics (mass density, energy density, charge density, etc.) there are also many types
of continuity equations.
ü Fluid dynamics
Motion of a fluid satisfies the continuity equation (1), there r is the mass density and the mass current j density is given by
j = ρv.
Liquids, unlike gases, are practically uncompressible, r = const, thus the continuity equation simplifies to
divv 0.
Dynamics of a non-viscous fluid is described by the Euler equation
∂v ∇P
+ Hv ∇L v − + g,
∂t ρ
where P is pressure and g is the gravity acceleration. Minus in front of the pressure gradient shows that the liquid accelerates
in the direction of a smaller pressure. The kinematic (streaming) term Hv “L v accounts for the change of the velocity at a
given point as a result of the media motion and makes the Euler equation nonlinear.
Flow of a viscous fluid obeys the Navier-Stokes equation. For an incompressible fluid it has the form
, Mathematical_physics-13-Partial_differential_equations.nb 3
∂v ∇P η
+ Hv ∇L v − +g+ ∆v,
∂t ρ ρ
where h is the viscosity coefficient.
Euler equation has a limited applicability because its solutions (for instance, for an object in a flow) have discontinuities of
the velocity in the stream. In these regions the viscosity becomes important and at high enough speeds the flow becomes
turbulent. Still, Euler equation generates Bernoilli equation that explains the lifting force on airplanes as well as the mecha-
nism of sound creation by voice bands.
Navier-Stokes equation is more realistic and it can be applied, for instance, to the problem of a drag force acting on a object
moving in a liquid. For small boby speeds u the quadratic terms can be dropped, equation becomes linear and can be solved
analytically for simple shapes such as sphere. The drag forces on a solid sphere and air bubble of radius R have the form
6 πηRu, solid sphere
F=µ
4 πηRu, air bubble.
The difference in the coefficient is due to different bounary conditions: The liquid sticks to the solid sphere but glides along
the surface of the air bubble that reduces the resistance.
ü Waves
Waves exist in compressible media such as fluids and solids. The main prerequisite for the waves' existence is a restoring
force. In fluids, there is only a restoring force with respect to change of the volume, thus the waves are waves of compression
/ expansion. In these waves the media displacement is collinear with the direction of the wave propagation, so that these
waves are called longitudinal waves. In solids, there is additionally a restoring force with respect to change of the shape of a
media element. The corresponding additional waves are transverse waves. Below we will consider small-amplitude longitudi-
nal waves in fluids.
In small-amplitude waves pressure and density only slightly deviate from their equilibrium values, so that one can define
P = P0 + δP, ρ = ρ0 + δρ
and linearize equations in small dP and dr, neglecting quadratic velocity terms because the velocity is small. The continuity
equation becomes
∂
δρ + ρ0 divv 0, (2)
∂t
whereas in the absence of viscosity from the Euler equation one obtains
∂v ∇ δP
− . (3)
∂t ρ0
The small changes of pressure and density are related by the adiabatic compressibility,
∂P
δP = δρ.
∂ρ S
One can obtain a single equation describing waves by differentiating the continuity equation over time and then using the
Euler equation. This yields the wave equation
∂P
∂2t δρ − c2 ∆δρ, c≡ ,
∂ρ S
where c is the speed of sound. Note that the speed of sound (that can be large) has no relation to the velocity of the media
(that is small). Since dP ~ dr, it satisfies the same equation,
, 4 Mathematical_physics-13-Partial_differential_equations.nb
∂2t δP − c2 ∆δP. (4)
Let us obtain the equation for the velocity. Taking rotor of the linearized Euler equation, one obtains
∂
∇ v 0,
∂t
thus velocity is a potential field, — ä v = 0, and it can be searched in the form
v = − ∇ φ.
Inserting it into the linearized Euler equation (3) and integrating it by removing the gradient one obtains
∑f dP c2
ã = dr. (5)
∑t r0 r0
Differentiating this over time and using the continuity equation, one obtains the same wave equation for the velocity potential,
∂2t φ − c2 ∆φ.
Let us now seek for the solution of the wave equation in the form of a plane wave
φ = ACos@ωt − k ⋅ r + φ0 D,
where A is the amplitude, w is the frequency and k is the wave vector that defines the direction of the wave propagation.
After substitution, differentiation, and cancellation one obtains the dispersion relation
ω2 − c2 k2 0.
The period T and wave length l are related to w and k by
2π 2π
ω= , k= .
T λ
One can see that w can be different while k adjusts to w as ck = w.
Let us now find the media velocity in a plane wave:
v = − ∇ φ = − A k Sin@ωt − k ⋅ r + φ0 D.
One can see that v is collinear to k, this is why such waves are called longitudinal waves.
Since the wave equation is linear, a superposition of plane waves moving in different directions with different frequencies
and wave lengths is also its solution. This illustrates the fact that there is much freedom in the solution of PDEs. The selec-
tion of the actual solution is done by the initial and boundary conditions.
In the important case of monochromatic waves whose time dependence is describe by a single frequency,
φ@r, tD = φ@rD Sin@ωtD,
or similar with Sin@wtD or Exp@ÂwtD, the wave equation simplifies to the time-independent stationary wave equation
ω
∆φ + k2 φ 0, k= .
c
In one dimension it becomes an ODE.
ü Thermal conductivity and diffusion
Temperature in an isolated body tends to equilibrate with time because the energy flows from warmer regions to colder
regions. Since the energy is conserved, its flow obeys the continuity equation
