Enseee Maths IBVPs 1 PDF
Enseee Maths IBVPs 1 PDF
Enseee Maths IBVPs 1 PDF
Solving IBVPs
with the Laplace transform
∂2u ∂2u
(E) elliptic + 2 =0 u(x,y) potential Fourier
∂x2 ∂y
∂ 2 u ∂u
(P) parabolic − =0 u(x,t) heat diffusion Laplace, Fourier
∂x2 ∂t
∂2u ∂2u
(H) hyperbolic − 2 =0 u(x,t) wave propagation Laplace, Fourier
∂x2 ∂t
We will come back later to a more systematic treatment and classification of PDEs.
The one dimensional examples exposed below intend to display some basic features and
differences between parabolic and hyperbolic partial differential equations. The constitutive
assumptions that lead to the partial differential field equations whose mathematical nature is
of prime concern are briefly introduced.
Problems of mathematical physics, continuum mechanics, fluid mechanics, strength of ma-
terials, hydrology, thermal diffusion · · ·, can be cast as IBVPs:
IBVPs: Initial and Boundary Value Problems
We will take care to systematically define problems in that broad framework. To introduce
the ideas, we delineate three types of relations, as follows.
1
Posted November 22, 2008; Updated, March 26, 2009
5
6 Solving IBVPs with the Laplace transform
chapter will consider prototypes of the two later types, and emphasize their physical meaning,
interpretation and fundamental differences.
We can not stress enough that
(P) for a parabolic equation, the information diffuses at infinite speed, and progressively,
while
(H) for a hyperbolic equation, the information propagates at finite speed and discontinuously.
These field equations should be satisfied within the body, say Ω. They are not required to
hold on the boundary ∂Ω. Ω is viewed as an open set of points in space, in the topological
sense.
In this chapter, the solutions u(x, t) are obtained through the Laplace transform in time,
u(x, t) → U (x, p) = L{u(x, t)}(p), (I.1.1)
and we shamelessly admit that the Laplace transform and partial derivative in space operators
commute,
∂ ∂
L{u(x, t)}(p) = L{ u(x, t)}(p) . (I.1.2)
∂x ∂x
In contrast, we will see that the use of the Fourier transform adopts the dual rule, in as far
as the Fourier variable is the space variable x,
u(x, t) → U (α, t) = F{u(x, t)}(α) . (I.1.3)
A qualitative motivation consists in establishing a correspondence between the time variable,
a positive quantity, and the definition of the (one-sided) Laplace transform which involves an
integration over a positive variable.
In contrast, the Fourier transform involves the integration over the whole real set, which is
interpreted as a spatial axis.
Attention should be paid to the interpretation of the physical phenomena considered. For
example, hyperbolic problems in continuum mechanics involves the speed of propagation of
elastic waves. This wave speed should not be confused with the velocity of particles. To
understand the difference, consider a wave moving over a fluid surface. When we follow the
phenomenon, we attach our eyes to the top of the wave. At two subsequent times, the particles
at the top of the waves are not the same, and therefore the wave speed and particle velocity are
two distinct functions of space and time. In fact, in the examples to be considered, the wave
speeds are constant in space and time.
v 0 ( t ) x=0
x
¶u ¶u
( x , t*) snapshot at t=t* ( x*, t ) observer at x=x*
¶t ¶t
v0 v0
0 0
0 x*=c t* x 0 t*=x*/c t
Figure I.1 A semi-infinite elastic bar is subject to a velocity discontinuity. Spatial and time profiles
p
of the particle velocity. The discontinuity propagates along the bar at the wave speed c = E/ρ
where E > 0 is Young’s modulus and ρ > 0 the mass density.
where E > 0 and ρ > 0 are respectively the Young’s modulus and mass density of the material.
The governing equations of dynamic linear elasticity are,
∂2u 1 ∂2u
(FE) field equation − = 0, t > 0, x > 0;
∂x2 c2 ∂t2
∂u
(IC) initial conditions u(x, 0) = 0; (x, 0) = 0 ; (I.2.2)
∂t
∂u
(BC) boundary condition (0, t) = v0 (t); radiation condition (RC) .
∂t
The field equation is obtained by combining
∂σ ∂2u
momentum balance − ρ 2 = 0;
∂x ∂x
(I.2.3)
∂u
elasticity σ=E ,
∂x
where σ is the axial stress.
The first initial condition (IC) means that the displacement is measured from time t = 0,
or said otherwise, that the configuration (geometry) at time t = 0 is used as a reference. The
second (IC) simply means that the bar is initially at rest.
The radiation condition (RC) is intended to imply that the mechanical information prop-
agates in a single direction, and that the bar is either of infinite length, or, at least, that the
signal has not the time to reach the right boundary of the bar in the time window of interest.
Indeed, any function of the form f1 (x − c t) + f2 (x + c t) satisfies the field equation. A function
of x − c t represents a signal that propagates toward increasing x. To see this, let us keep the
eyes on some given value of f1 , corresponding to x − c t equal to some constant. Then clearly,
sine the wave speed c is a positive quantity, the point we follow moves in time toward increasing
x.
The solution is obtained through the Laplace transform in time,
and we admit that the Laplace transform and partial derivative in space operators commute,
∂ ∂
L{u(x, t)}(p) = L{ u(x, t)}(p) . (I.2.5)
∂x ∂x
Therefore,
∂ 2 U (x, p) 1 2 ∂u
(FE) − p U (x, p) − p u(x, 0) − (x, 0) =0
∂x2 c2 | {z } |∂t {z }
=0, (IC)
=0, (IC) (I.2.6)
p p
⇒ U (x, p) = a(p) exp(− x) + b(p) exp(+ x) .
c | {z c }
b(p)=0, (RC)
Benjamin LORET 9
The term exp(−p x/c) gives rise to a wave which propagates toward increasing value of x, and
the term exp(p x/c) gives rise to a wave which propagates toward decreasing value of x. Where
do these assertions come from? There is no direct answer, simply they can be checked on the
result to be obtained. Thus the radiation condition implies to set b(p) equal to 0.
In turn, taking the Laplace transform of the (BC),
x x ∂u x
u(x, t) = V0 (t − ) H(t − ), (x, t) = V0 H(t − ) . (I.2.9)
c c ∂t c
The analysis of the velocity of the particles reveals two main characteristics of a partial differ-
ential equation of the hyperbolic type in non dissipative materials, Fig. I.1:
(H1). the mechanical information propagates at finite speed, namely c, which is
therefore termed elastic wave speed;
(H2). the wave front carrying the mechanical information propagates undistorted
and with a constant amplitude.
Besides, this example highlights the respective meanings of elastic wave speed c and velocity
of the particles ∂u/∂t(x, t).
1.2 1.2
L=10-2 m D=10-9 m2/s tF=105 s
x/L=10-2
(T(x,t)-Tµ)/(T0-Tµ)
(T(x,t)-Tµ)/(T0-Tµ)
x/L=1
0.4 10-1 0.4
10-2 L=10-2 m
10-3 D=10-9 m2/s
t/tF=10-4 tF=105 s
0 0
0 0.002 0.004 0.006 0.008 0.01 0 20000 40000 60000 80000 100000
distance from left boundary x (m) time t(s)
Figure I.2 Heat diffusion on a semi-infinite bar subject to a heat shock at its left boundary x = 0.
The characteristic time tF that can be used to describe the information that reaches the position
x = L is equal to L2 /D.
10 Solving IBVPs with the Laplace transform
I.3 A parabolic PDE: diffusion of a heat shock
As a second example of PDE, let us consider the diffusion of a heat shock in a semi-infinite
body.
A semi-infinite plane, or half-space, or bar, is subject to a given temperature T = T 0 (t) at
its left boundary x = 0. The thermal diffusion is one dimensional. The initial temperature
along the bar is uniform, T (x, t) = T∞ . It is more convenient to work with the field θ(x, t) =
T (x, t) − T∞ than with the temperature itself.
For a rigid material, in absence of heat source, the energy equation links the divergence of
the heat flux Q [unit : kg/s3 ], and the time rate of the temperature field θ(x, t), namely
∂θ
energy equation div Q + C = 0, (I.3.1)
∂t
or in cartesian axes and using the convention of summation over repeated indices ( the index i
varies from 1 to n in a space of dimension n),
∂Qi ∂θ
+C = 0. (I.3.2)
∂xi ∂t
Here C > [unit : kg/m/s2 /◦ K] is the heat capacity per unit volume. The heat flux is related to
the temperature gradient by Fourier law
or componentwise,
∂θ
Qi = −kT , (I.3.4)
∂xi
where kT > 0 [unit : kg×m/s3 /◦ K] is the thermal conductivity.
Inserting Fourier’s law in the energy equation shows that a single material coefficient D
[unit : m2 /s],
kT
D= > 0, (I.3.5)
C
that we shall term diffusivity, appears in the field equation.
The initial and boundary value problem (IBVP) is thus governed by the following set of
equations:
∂2θ ∂θ
field equation (FE) D 2
− = 0, t > 0, x > 0;
∂x ∂t
initial condition (IC) θ(x, 0) = 0 ; (I.3.6)
The radiation condition is intended to imply that the thermal information diffuses in a single
direction, namely toward increasing x, and that the bar is either of infinite length.
The solution is obtained through the Laplace transform
then
x2
Z t
x 1
θ(x, t) = √ θ0 (t − u) exp − du . (I.3.11)
2 Dπ 0 u3/2 4 Du
With the further change of variable
x2
u→v with v2 = , (I.3.12)
4D u
the inverse Laplace transform takes the form,
x2
Z ∞
2
θ(x, t) = √ √ θ0 (t − ) exp(−v 2 ) dv , (I.3.13)
π x/2 D t 4 D v2
Note however, that a modification of the Fourier’s law, or of the energy equation, which
goes by the name of Cataneo, provides finite propagation speeds. The point is not considered
further here.
In contrast to wave propagation, where time and space are involved in linear expressions
x ± c t, here the space variable is associated with the square root of the time variable. Therefore
diffusion over a distance 2 L requires a time interval four times larger than the time interval of
diffusion over a length L.
and by the mass balance, which, in terms of concentration and absolute flux J = c v s , writes,
∂c
balance of mass + div J = 0 . (I.4.4)
∂t
Benjamin LORET 13
(c(x,t)-cµ)/(c0-cµ)
0.8 t/tF=1
0.6 0.5
0.4
0.2 Pe=10
0.2 0.1 tH/tF=0.1
10-2
0
0 2 4 6 8 10
distance from left boundary x/L
1
(c(x,t)-cµ)/(c0-cµ)
0.8
0.6
0.4 t/tF=1
0.5 Pe=2
0.2 0.2 tH/tF=1/2
0.1
10-2
0
0 2 4 6 8 10
distance from left boundary x/L
1.2
(c(x,t)-cµ)/(c0-cµ)
0.8
0.4 t/tF=1
Pe=0
0.5 tH/tF=µ
0.2
0.1
0
0 0.02 0.04 0.06 0.08 0.1
distance from left boundary x/L
Figure I.3 Advection-diffusion along a semi-infinite bar of a species whose concentration is subject at
time t = 0 to a sudden increase at the left boundary x = 0. The fluid is animated with a velocity
v such that the Péclet number Pe is equal to 10, 2 and 0 (pure diffusion) respectively. Focus
is on the events that take place at point x = L. The time which characterizes the propagation
phenomenon is tH = L/v while the characteristic time of diffusion is tF = L2 /D, with D the
diffusion coefficient. Therefore Pe = L v/D = tH /tF .
diffusive advective
terms term
z }| { z }| {
∂2c ∂c ∂c
(FE) field equation D 2− = v , t > 0, x > 0;
∂x ∂t ∂x (I.4.5)
(IC) initial condition c(x, 0) = ci ;
-v 0 (t)
t<0
t=0
particle particle
velocity=0 velocity - v 0 ( t )
t=t*>0
x*=c t*
wave front
x=0 x
Figure I.4 An elastic dart moving at speed −v0 hits a rigid target at time t = 0. The shock then
propagates along the dart at the speed of elastic longitudinal waves. The part of the dart behind
the wave front is set to rest and it undergoes axial compression: although the analysis here is
one-dimensional, we have visualized this aspect by a (mechanically realistic) lateral expansion.
The dart is assumed to be made in a linear elastic material, in which the longitudinal
waves propagate at speed c. The equations of dynamic linear elasticity governing the axial
displacement u(x, t) of points of the dart are,
∂2u 1 ∂2u
(FE) field equation − = 0, t > 0, x > 0;
∂x2 c2 ∂t2
∂u
(IC) initial conditions u(x, 0) = 0; (x, 0) = −v0 ; (1)
∂t
∂u
(BC) boundary condition u(0, t) = 0; radiation condition (RC) (∞, t) = 0 .
∂x
Find the displacement u(x, t) and give a vivid interpretation of the event.
Solution:
The solution is obtained through the Laplace transform in time,
and we admit that the Laplace transform and partial derivative in space operators commute,
∂ ∂
L{u(x, t)}(p) = L{ u(x, t)}(p) . (3)
∂x ∂x
16 Solving IBVPs with the Laplace transform
Therefore,
∂ 2 U (x, p) 1 2 ∂u
(FE) − p U (x, p) − p u(x, 0) − (x, 0) =0
∂x2 c2 | {z } |∂t {z }
=0, (IC) (4)
=−v0 , (IC)
d2 v0 p2 v0
U (x, p) + = U (x, p) + 2 .
dx2 p2 c2 p
The change of notation, from partial to total derivative wrt space, is intended to convey the
idea that the second relation is seen as an ordinary differential equation in space, where the
Laplace variable plays the role of a parameter. Then
v0 p p
U (x, p) + = a(p) exp( x) + b(p) exp(− x) . (5)
p2 c c
The first term in the solution above would give rise to a wave moving to left, which is prevented
by the radiation condition:
(RC) a(p) = 0
v0
(BC) 0+ = b(p) exp(0) (6)
p2
v0 p
U (x, p) = 2 exp(− x) − 1 .
p c
0
0 1 2 3 4 t
-1
f2 (t)
1
0
0 1 2 3 4 t
f3 (t)
1
0
0 a 2a 3a 4a t
Figure I.5 Periodic functions; for f3 (t), the real a > 0 is strictly positive.
Proof :
1. These functions are periodic, for t > 0, with period T = 2:
Z 2 Z 1 Z 2
1 −t p 1 −t p −t p
L{f1 (t)}(p) = e f1 (t) dt = e dt + e (−1) dt . (2)
1 − e−2p 0 1 − e−2p 0 1
x
x=0 x=L
Figure I.6 An elastic bar, fixed at its left boundary, is hit by a sudden load at its right boundary at
time t=0.
The bar is made of a linear elastic material, with a Young’s modulus E and a section S,
and the longitudinal waves propagate at speed c.
The equations of dynamic linear elasticity governing the axial displacement u(x, t) of the
points inside the bar are,
∂2u 1 ∂2u
(FE) field equation − = 0, t > 0, x ∈]0, L[;
∂x2 c2 ∂t2
∂u
(IC) initial conditions u(x, 0) = 0; (x, 0) = 0 ; (1)
∂t
∂u F0
(BC) boundary conditions u(0, t) = 0; (L, t) = H(t) .
∂x ES
Find the displacement u(L, t) of the right boundary and give a vivid interpretation of the
phenomenon.
Solution:
The solution is obtained through the Laplace transform in time,
and we admit that the Laplace transform and partial derivative in space operators commute,
∂ ∂
L{u(x, t)}(p) = L{ u(x, t)}(p) . (3)
∂x ∂x
Therefore,
∂ 2 U (x, p) 1 2 ∂u
(FE) − p U (x, p) − p u(x, 0) − (x, 0) =0
∂x2 c2 | {z } |∂t {z }
=0, (IC)
=0, (IC) (4)
d2 p2
U (x, p) = 2 U (x, p) .
dx2 c
The change of notation, from partial to total derivative wrt space, is intended to convey the
idea that the second relation is seen as an ordinary differential equation in space, where the
Laplace variable plays the role of a parameter. Then
p p
U (x, p) = a(p) cosh ( x) + b(p) sinh ( x) . (5)
c c
Benjamin LORET 19
The two unknown functions of p are defined by the two boundary conditions,
(BC)1 U (x, p) = 0 ⇒ a(p) = 0
∂U F0 1 p p
(BC)2 (L, p) = = b(p) cosh ( L)
∂x ES p c c (6)
p
F0 1 sinh ( c x)
U (x, p) = c .
ES p2 cosh ( p L)
c
Let T = L/c be the time for the longitudinal wave to travel the length of the bar. Then, the
displacement of the right boundary x = L is
F0 1 F0
u(L, t) = c L−1 { 2 tanh(T p)}(t) = 2 L f3 (t) , (7)
ES p ES
where f3 (t) is a periodic function shown on Fig. I.7. Use has been made of the previous exercise
with a = 2 T . Therefore the velocity of the boundary x = L reads,
∂u F0 df (t) F0
(L, t) = 2 L = ±c . (8)
∂t ES | dt
{z } ES
±1/2T
u(L,t) displacement ¶u
(L, t ) velocity
¶t
F
2L 0 F
c 0
ES ES
0 0
0 2T 4T 6T 8T t 0 2T 4T 6T 8T t
F
-c 0
ES
Figure I.7 Following the shock, the ensuing longitudinal wave propagates back and forth, giving rise
to a periodic motion with a period 4 T equal to four times the time required for the longitudinal
elastic wave to travel the bar.
The interpretation of these longitudinal vibrations is a bit tricky and it goes as follows:
- at time t = 0, the right boundary of the bar is hit, and the corresponding signal propagates
to the left at speed c;
- it needs a time T to reach the fixed boundary. When the wave comes back, it carries the
information that this boundary is fixed;
- this information is known by the right boundary after a time 2 T . Since it was not known
before, this point was moving to the right, with the positive velocity indicated by (8).
Immediately as the information is known, the right boundary stops moving right, and in
fact, changes the sign of its velocity, again as indicated by (8);
- the mechanical information “the right boundary is subject to a fixed traction” then travels
back to the left, and hits the fixed boundary at time 3 T . When the wave comes back,
it carries the information that this boundary is fixed. This information reaches the right
boundary at time 4 T , where in fact its displacement just vanishes;
20 Solving IBVPs with the Laplace transform
elongation
u(L,t)
wave front
0 0 T 2T 3T 4T t
wave front
0
2T elongation
3T
4T
Figure I.8 Two other equivalent illustrations of the back and forth propagation of the wave front in
the finite bar with left end fixed and right end subject to a given traction.
- the elongation increases linearly in time from t = 0 to reach its maximum at t = 2 T , and
then decreases up to t = 4 T where it vanishes;
- this succession of events, with periodicity 4 T , repeats indefinitely.
Figs. I.7 and I.8 display the position of the wave front and elongation at various times within
a period.
Benjamin LORET 21
x
x w(x)
x=0 x=L
¶w deformed beam
impact ( x , t = 0)
¶t
Figure I.9 The supports at the boundaries are bilateral, that is, they prevent up and down vertical
motions of the beam.
The equations of dynamic linear elasticity governing the transverse displacement w(x, t) of
the beam are,
∂2w 4
2 ∂ w
(FE) field equation + b = 0, t > 0, x ∈]0, L[;
∂t2 ∂x4
∂w x
(IC) initial conditions w(x, 0) = 0; (x, 0) = V0 sin π ; (1)
∂t L
∂2w ∂2w
(BC) boundary conditions w(0, t) = w(L, t) = 0; (0, t) = (L, t) = 0 .
∂x2 ∂x2
p
The coefficient b involved in the field equation is equal to EI/(ρ S), where ρ is the mass
density of the material and S the section of the beam, all quantities that we will consider as
constants.
The beam remains uncharged. The transverse shock generates transverse vibrations w(x, t).
Describe these so-called free vibrations.
Solution:
The solution is obtained through the Laplace transform in time,
w(x, t) → W (x, p) = L{w(x, t)}(p), (2)
and we admit that the operators transform and partial derivative in space Laplace commute,
∂ ∂
L{w(x, t)}(p) = L{ w(x, t)}(p) . (3)
∂x ∂x
Therefore,
∂ 4 W (x, p) ∂w
(FE) b2 + p2 W (x, p) − p w(x, 0) − (x, 0) =0
∂x4 | {z } |∂t {z }
=0, (IC)
=V0 sin (π x/L), (IC) (4)
d4 W (x, p) x
b2 + p2 W (x, p) = V0 sin (π ).
dx4 L
22 Solving IBVPs with the Laplace transform
The change of notation, from partial to total derivative wrt space, is intended to convey the
idea that the second relation is seen as an ordinary differential equation in space, where the
Laplace variable plays the role of a parameter.
The solution of this nonhomogeneous linear differential equation is equal to the sum of the
solution of the homogeneous equation (with zero rhs), and of a particular solution.
The solution of the homogeneous equation is sought in the format
p W (w, p) = C(p) exp(α x),
yielding four complex solutions α = ± (1 ± i) β, with β(p) = p/(2 b), and with an priori
complex factor C(p). Summing these four solutions, the real part may be rewritten in the
format,
W hom (x, p) = c1 (p) cos (β x) + c2 (p) sin (β x) exp(β x)
(5)
+ c3 (p) cos (β x) + c4 (p) sin (β x) exp(−β x) ,
x V0
W par (x, p) = c5 (p) sin (π ), c5 (p) = . (6)
L 2 π4 2
p +b 4
L
We will need the first two derivatives of the solution,
W (x, p) = c1 (p) cos (β x) + c2 (p) sin (β x) exp(β x)
+ c3 (p) cos (β x) + c4 (p) sin (β x) exp(−β x)
x
+ c5 (p) sin (π )
L
d
W (x, p) = β (c1 (p) + c2 (p)) cos (β x) + (−c1 (p) + c2 (p)) sin (β x) exp(β x)
dx
+ β (−c3 (p) + c4 (p)) cos (β x) − (c3 (p) + c4 (p)) sin (β x) exp(−β x) (7)
π x
+ c5 (p) cos (π )
L L
d2 2
W (x, p) = 2 β c 2 (p) cos (β x) − c 1 (p) sin (β x) exp(β x)
dx2
+ 2 β 2 − c4 (p) cos (β x) + c3 (p) sin (β x) exp(−β x)
π2 x
− c5 (p) 2
sin (π ) ,
L L
The four boundary conditions are used to obtain the four unknowns,
c1 (p) + c3 (p) = 0
c2 (p) − c4 (p) = 0
c1 (p) cos (β L) + c2 (p) sin (β L) exp(β L) + c3 (p) cos (β L) + c4 (p) sin (β L) exp(−β L) = 0
(c2 (p) cos (β L) − c1 (p) sin (β L) exp(β L) + (−c4 (p) cos (β L) + c3 (p) sin (β L) exp(−β L) = 0 ,
(8)
Benjamin LORET 23
yielding c1 = −c3 , c2 = c4 , and a 2 × 2 linear system for c1 and c2 ,
The determinant of this system, 4 (sinh 2 (β L) + sin 2 (β L)), does not vanish, and therefore,
ci = 0, i = 1 − 4, and finally,
V0 x
W (x, p) = W par (x, p) = sin (π ). (10)
2 2 π4 L
p +b
L4
Since
a
L{sin (a t)}(p) = , (11)
p2 + a2
the transverse displacement, and transverse velocity,
L2 b π2 x ∂w b π2 x
w(x, t) = V0 2
sin 2
t sin (π ), (x, t) = V0 cos 2
t sin (π ) , (12)
bπ L L ∂t L L
are periodic with a frequency,
π b
, (13)
2 L2
that is inversely proportional to the square of the length of the beam: doubling the length of the
beam divides by four its frequency of vibration. On the other hand, the higher the transverse
stiffness, the higher the frequency.
Since the solution displays a separation of the space and time variables, all points of the
beam vibrate in phase, and the beam keeps its spatial shape for ever.
24 Solving IBVPs with the Laplace transform
x=0 x=0
x=1 x=1
given
temperature
q(0,t)=0
Figure I.10 Heat diffusion in a finite bar subject to given temperature at its ends.
The temperature at the two ends of a finite bar x ∈ [0, 1] is maintained at a given value,
say T (0, t) = T (1, t) = T0 . The initial temperature along the bar T (x, 0)=T (x, 0) is a function
of space, say T (x, 0). It is more convenient to work with the field θ(x, t) = T (x, t) − T 0 than
with the temperature itself. Moreover, to simplify the notation, scaling of time and space has
been used so as to make the thermal conductivity equal to one.
The initial and boundary value problem (IBVP) is thus governed by the following set of
equations:
∂2θ ∂θ
field equation (FE) − = 0, t > 0, x ∈]0, 1[;
∂x2 ∂t
initial condition (IC) θ(x, 0) = θ0 sin (2 π x) ; (1)
Solution:
The solution is obtained through the Laplace transform
Therefore,
∂ 2 Θ(x, p)
(FE) − p Θ(x, p) − θ(x, 0) = 0. (3)
∂x2 | {z }
=θ0 sin (2 π x), (IC)
The solution to this nonhomogeneous linear equation is the sum of the solution to the ho-
mogeneous equation, and of a particular solution. The latter is sought in the form of the
inhomogeneity. Therefore,
θ0
Θpar (x, p) = c(p) sin (2 π x), c(p) = . (4)
p + 4 π2
Therefore,
√
px
√ θ0
Θ(x, p) = c1 (p) e + c2 (p) e− px
+ sin (2 π x) . (5)
p + 4 π2
Benjamin LORET 25
The boundary conditions imply the unknowns c 1 (p) and c2 (p) to vanish. Therefore
θ0
Θ(x, p) = sin (2 π x) , (6)
p + 4 π2
and
2
θ(x, t) = θ0 sin (2 π x) e−4 π t
H(t) . (7)
The spatial temperature profile remains identical in time, but its variation with respect to the
temperature imposed at the boundaries decreases and ultimately vanishes. In other words, the
information imposed at the boundaries penetrates progressively the body.
Since the temperature is imposed at the ends, the heat fluxes ∇θ(x, t) at these ends x = 0, 1
can be seen as the response of the structure to a constraint.
26 Solving IBVPs with the Laplace transform
Proof:
√
1.1 The function p should be made uniform by defining appropriate cuts. One can then
√
calculate the derivatives of F (p) = exp(− p),
√ √ √
√
− p d e− p d2 e− p e− p
F (p) = e , F (p) = − √ , 2
F (p) = + √ . (3)
dp 2 p dp 4p 4p p
Therefore,
d2 d
4p 2
F (p) + 2 F (p) − F (p) = 0 . (4)
dp dp
If f (t) has Laplace transform F (p), let us recall the rules,
d
L{t f (t)}(p) = − F (p),
dp
d2
L{t2 f (t)}(p) = (−1)2 2 F (p),
dp (5)
d 2 d2
L{ t f (t) }(p) = p L{t2 f (t)}(p) − t2 f (t) (t = 0) = p 2 F (p) .
dt | {z } dp
=0
That the second term on the rhs of the last line above is really zero need to be checked, once the
solution has been obtained. Collecting these relations, the differential equation in the Laplace
domain (4) can be transformed into a differential equation in time,
d2
4 t f (t) − 2 t f (t) − f (t) = 0 , (6)
dt
which, upon expansion of the first term, becomes,
df 3 1 1 1
+ − 2 dt = 0 , (7)
f 2t 4t
and thus can be integrated to,
c
f (t) = e−1/(4t) . (8)
t3/2
1.2 The constant c can be obtained as follows. We will need a generalized Abel theorem, that
we can state as follows:
Assume that, for large t, the two functions f (t) and g(t) are sufficiently close, f (t) ' g(t),
for t 1. Then their Laplace transforms t → p are also close for small p, namely F (p) ' G(p),
for p ∼ 0.
Benjamin LORET 27
Consider now,
√
c −1/(4t) d e− p
t f (t) = e , L{t f (t)}(p) = − F (p) = √ . (9)
t1/2 dp 2 p
Then
c Γ(1/2)
for t 1, t f (t) ' ⇒ L{t f (t)}(p) ' c ,
t1/2 p1/2
√ (10)
e− p 1
for p ∼ 0, √ ' √ .
2 p 2 p
Requiring the transforms in these two lines to be equal as indicated by (9), and by the gener-
√
alized Abel theorem, yields c Γ(1/2) = 1/2, and therefore c = 1/(2 π). Therefore
√ 1 1
L−1 (exp(− p)(t) = √ exp(− ) . (11)
2 π t3 4t
a2
1 Z t
√ a
L−1 exp(− p a) (t) = H(t − u) √ exp(− ) du
p 0 2 π u3 4u
Z ∞ (14)
2 −v 2 2 2
= √ √ e dv, with v = a /(4 u) .
π a/(2 t)
2
28 Solving IBVPs with the Laplace transform
where a > 0, b ≥ 0 are constants, x ≥ 0 is the space coordinate and p the Laplace variable
associated with time t. Show that the inverse reads,
s s
exp(b t) b x √ b x √
q(x, t) = exp (x ) erfc( √ + b t) + exp ( − x ) erfc( √ − b t) . (2)
2 a 2 at a 2 at
Proof:
The function is first made uniform by introducing a branch cut along the negative axis
< p ≤ 0, and the definitions,
√ q √
p = |p| exp(i θ), θ ∈] − π, π], p = |p| exp(i θ/2) (⇒ < p ≥ 0) . (3)
b x2
Z ∞
2
q(x, t) = exp(b t) √ √ exp(−v 2 − ) dv , (5)
π x/2 a t 4 a v2
since the Laplace transform of a convolution product is equal to the product of the Laplace
transforms.
This integral seems out of reach, except if b = 0, where
x
q(x, t) = exp(b t) erfc( √ ) . (6)
2 at
Still it is provided in explicit form in Abramowitz and Stegun [1964], p. 304, for any b.
Since we want to obtain the result on our own, even for b 6= 0, we take a step backward and
consider the inversion in the complex plane through an appropriate closed contour C = C(R, ),
that respects the branch cut and includes the pole p = b. A preliminary finite contour is shown
on Fig. I.11. The residue theorem yields,
Z Z Z Z Z
1 c+i R
+ + + + exp(tp) Q(x, p) dp
2iπ c−i R CR C C+ C−
I q (7)
1
= exp(t p) Q(x, p) dp = exp(b t − x b/a) ,
2iπ C
where c is an arbitrary real, which is required to be greater than b for a correct definition of the
inverse Laplace transform. For vanishingly small radius , the contour C does not contribute
as long as b 6= 0. If b = 0, then the inverse Laplace transfrm is read directly from (5). Moreover,
Benjamin LORET 29
Im p
c + iR
CR R
q
p = p exp(ip/2)
e
C+ b c
C- Re p
Ce
p = p exp(-i p/2)
c - iR
since Q(x, p) tends to 0 for large p in view of (3), the contour C R does not contribute either for
large R, in view of Jordan’s lemma. Consequently,
Z c+i R
1
q(x, t) = lim exp(t p) Q(x, p) dp
R→∞ 2iπ c−i R
q Z Z (8)
1
= exp(b t − x b/a) − lim + exp(t p) Q(x, p) dp .
R→∞, →0 2 i π C+ C−
√ √ √
On the branch cut, p = −r < 0, but p is equal to i r on the upper part and to −i r on the
lower part. Therefore, the integrals along the branch cut become,
Z p Z p
1 0 exp(−i x r/a) 1 ∞ exp(i x r/a)
I = exp(−t r) (−dr) + exp(−t r) (−dr)
2iπ ∞ −r − b 2iπ 0 −r − b
Z ∞
1 exp(−t r) q q
= ( exp(i x r/a) − exp(−i x r/a)) dr
2iπ 0 r+b
exp(−t ρ2 )
Z ∞
1 √ √
= ( exp(i x ρ/ a) − exp(−i x ρ/ a)) ρ dρ
iπ 0 ρ2 + b
Z ∞
1 √ ρ
= exp ( − t ρ2 + i x ρ/ a) 2 dρ .
i π −∞ ρ +b
(9)
In an attempt to see the error function emerging, we make a change of variable that transforms,
to within a constant, the argument of the exponential into a square, namely
√ x
ρ −→ v = tρ − i v0 , v0 ≡ √ . (10)
2 at
Then
x 1
Z ∞−i v0 exp(−v 2 ) exp(−v 2 ) √
I = exp(−( √ )2 ) × + dv, v ± ≡ v0 ± bt. (11)
2 at 2iπ −∞−i v0 v + i v+ v + i v−
30 Solving IBVPs with the Laplace transform
Re v Re v -iv- Re v
-iv-
-iv0 -iv0 -iv0
-iv+
Figure I.12 Contours of integration associated with the integrals (12). Note that (a) v + is always
larger than v0 > 0, while v − is (b) either between 0 and v 0 or (c) negative.
First observe that, by an appropriate choice of integration path, Fig. I.12, and application
of the residue theorem,
1
Z ∞−i v0 exp(−v 2 ) exp(−v 2 )
1 0
Z ∞ for v ± = v + or v − < 0,
±
dv = ±
dv +
2iπ −∞−i v0 v + iv 2iπ
−∞ v + i v exp((v − )2 ) for v ± = v − > 0,
(12)
we are left with integrals on the real line. The basic idea is to insert the following identity,
Z ∞
1
exp(−u X 2 ) du = 2 , X 6= 0 , (13)
0 X
in the integrals to be estimated and to exchange the order of integration. With this idea in
mind, the following straightforward transformations are performed,
exp(−v 2 )
Z ∞
1
dv
2iπ −∞ v + i v±
v± exp(−v 2 )
Z ∞
= − dv (use (13) )
π 0 v 2 + (v ± )2
v±
Z ∞ √
= − exp ( − v 2 − u (v 2 + (v ± )2 )) du dv, (du dv → dv du, v→v 1 + u) (14)
π 0
v± exp(−u (v ± )2 )
Z ∞
= − √ √ du (u → (sinh v)2 , |v ± | cosh v → w)
2 π 0 1+u
sgn(v ± )
Z ∞
± 2 2
= − exp ((v ) ) √ exp(−w2 ) dw ,
2 π |v± |
and therefore
exp(−v 2 ) sgn(v ± )
Z ∞
1
dv = − exp ((v ± )2 ) erfc(|v ± |) . (15)
2iπ −∞ v + i v± 2
Collecting the results (12) and (15), and using the property (I.3.17) 3 ,
1
Z ∞−i v0 exp(−v 2 ) − 12 exp ((v + )2 ) erfc(v + ) for v ± = v + ,
dv =
2iπ −∞−i v0 v + i v± − 1 exp ((v − )2 ) erfc(v − ) + exp ((v − )2 ) for v ± = v − ,
2
(16)
Finally, the inversion formula deduces from (8) and (16).