Search in book...
Toggle Font Controls
Create new playlist

Name your new playlist

Playlist description (optional)
Sign In

Email address

Password

Forgot Password?

or

Continue with Facebook

Continue with Google
Sign Up

Full Name

Email address

Confirm Email Address

Password

or

Continue with Facebook

Continue with Google

Chapter 8

Differential Equations

8.1 Important Differential Equations

In this section we will look at some very important differential equations that we often come across in practical physics and fluid mechanics. We will not derive each equation here. For this, the interested reader is referred to textbooks that discuss the underlying physics to the differential equations. However, several of these equations will be treated throughout this book and will therefore be derived. Wherever possible the equations will be briefly introduced and then given in a standard form using f denoting the actual function to solve, c as the constants, and the independent variables x, y, and t. Occasionally, we will use the denotation c² instead of simply writing c. In this case we want to indicate that the constant has to be positive as c² will always be.

8.1.1 Transport Equation

The first differential equation we will introduce is the transport equation. The transport equation is the equation that describes the movement of a function along time and space. Imagine throwing a stone in the water. The stone will cause water waves that will travel away from the point where the stone sank. If you look very closely you will see that these waves are perturbations of the surface that simply travel along the surface.

As you sit outside of the water and watch the waves there are two potential ways of looking at this phenomenon. First, you could focus on a given point of the water surface that is at rest initially. As the wave propagates toward this point you will see that the point is lifted upwards first and then drops downwards again, pretty much like a sine wave. If there are a couple of waves this up and down movement will occur a couple of times until the waves have completely passed and the surface is yet again at rest.

The second way of looking at it would be to focus, e.g., the highest point of the first wave and move your eyes along with this point. You would then see that the waves (which look like a sine function) will actually not change over time. They will seem "fixed" and merely travel along the surface. If you look carefully you will see that the height of the highest point of the wave will slowly drop over time. This is due to the fact that the wave will experience damping and therefore decays over time. Without decay, the wave would travel on forever.

This is just a very simple example of the transport equation. Imagine any function, any physical shape or contour which has a "fixed" shape. Let us denote this shape or profile by the function f (x). This profile travels along a physical direction. If damping is considered, the function will slowly decay over time. This effectively makes the function dependent both on time t and space x as its independent variables. Therefore we are looking at a function f (x, t). Whenever you have such an effect, you are looking at a phenomenon that is described by the transport equation.

The transport equation is given by

∂f(x,t)∂t+c∂f(x,t)∂x+af(x,t)∂f∂t+c∂f∂x+af==00 $\begin{array}{l} \frac{\partial f (x, t)}{\partial t} + c \frac{\partial f (x, t)}{\partial x} + a f (x, t) & = & 0 \\ \frac{\partial f}{\partial t} + c \frac{\partial f}{\partial x} + a f & = & 0 \end{array}$

si36_e (Eq. 8.1)

where c is the velocity at which the profile travels and a is the damping factor. The damping factor will let the function decay over time. As we will see, this decay process is really only a function of time.

8.1.2 Heat Conduction

Heat conduction is among the oldest problems that yields partial differential equations. The resulting differential equation is referred to as the heat conduction equation. It is derived in section 6.9.2 and given by Eq. 6.82 as

∂T∂t∂f∂t==αΔTc2Δf $\begin{array}{l} \frac{\partial T}{\partial t} & = & α Δ T \\ \frac{\partial f}{\partial t} & = & c^{2} Δ f \end{array}$

si37_e (Eq. 8.2)

One-Dimensional Heat Conduction

Reduced to one dimension, the heat conduction equation becomes

∂T∂t∂f∂t==α∂2T∂x2c2∂2f∂x2 $\begin{array}{l} \frac{\partial T}{\partial t} & = & α \frac{\partial^{2} T}{\partial x^{2}} \\ \frac{\partial f}{\partial t} & = & c^{2} \frac{\partial^{2} f}{\partial x^{2}} \end{array}$

si38_e (Eq. 8.3)

8.1.3 Diffusion Equation

Very closely related to heat conduction is mass transport by means of diffusion. The underlying PDE is derived in section 6.9.3 and given by Eq. 6.90 as

∂ρ∂t∂f∂t==DΔρc2Δf $\begin{array}{l} \frac{\partial ρ}{\partial t} & = & D Δ ρ \\ \frac{\partial f}{\partial t} & = & c^{2} Δ f \end{array}$

si39_e (Eq. 8.4)

Comparing Eq. 8.2 and Eq. 8.4 the analogy becomes obvious: heat and mass transport are governed by the same underlying PDE which simply has different physical meanings. Therefore, finding solutions to one equation will inevitably also yield solutions to the other.

One-Dimensional Diffusion

Reduced to one dimension, the diffusion equation becomes

∂ρ∂t∂f∂t==D∂2ρ∂x2c2∂2f∂x2 $\begin{array}{l} \frac{\partial ρ}{\partial t} & = & D \frac{\partial^{2} ρ}{\partial x^{2}} \\ \frac{\partial f}{\partial t} & = & c^{2} \frac{\partial^{2} f}{\partial x^{2}} \end{array}$

si40_e (Eq. 8.5)

8.1.4 Wave Equation

8.1.4.1 Introduction

The wave equation is one of the most important equations in mechanics. It describes not only the movement of strings and wires, but also the movement of fluid surfaces, e.g., water waves. The wave equation is surprisingly simple to derive and not very complicated to solve although it is a second-order PDE.

We will derive the wave equation using the model of the suspended string (see Fig. 8.1). We will apply a few simplifications. First, the string is only assumed to move along the direction of the y-axis. Second, we will assume the problem to be planar, i.e., there is no movement in the z-direction. Third, we will assume the string to be perfectly bendable, i.e., the tangential forces F₁ and F₂ are effectively tangential to the shape of the string.

f08-01-9781455731411 — Fig. 8.1 Derivation of the wave equation using the suspended string of length L.

8.1.4.2 Force Balance

The force balance along the y-axis is given as

FyF2,yF1,y===F2,y−F1,yF2sinβF1sinα $\begin{array}{l} F_{y} & = & F_{2, y} - F_{1, y} \\ F_{2, y} & = & F_{2} sin β \\ F_{1, y} & = & F_{1} sin α \end{array}$

si41_e (Eq. 8.6)

As there is (per definition) no movement in the x-direction the forces along the x-axis balance to a constant value, F_2,x = F_1,x = F₀

F2,x−F1,xF2,x==0F1,x=F0 $\begin{array}{l} F_{2, x} - F_{1, x} & = & 0 \\ F_{2, x} & = & F_{1, x} = F_{0} \end{array}$

si42_e

Geometrically we can see that

F1,y=F1,xtanα=F0tanα $F_{1, y} = F_{1, x} tan α = F_{0} tan α$

(Eq. 8.7)

F2,y=F2,xtanβ=F0tanβ $F_{2, y} = F_{2, x} tan β = F_{0} tan β$

(Eq. 8.8)

Inserting Eq. 8.7 and Eq. 8.8 into Eq. 8.6 yields

Fy=F0(tanβ−tanα) $F_{y} = F_{0} (tan β - tan α)$

(Eq. 8.9)

8.1.4.3 Newton’s Second Law of Motion

Now we apply Newton’s¹ second law of motion (see Eq. 11.1), which states

dF==dmadm∂2u∂t2 $\begin{array}{l} d F & = & dm a \\ = & dm \frac{\partial^{2} u}{\partial t^{2}} \end{array}$

si46_e (Eq. 8.10)

which we expressed as a function of the string’s contour line u (x). The mass of the string in the interval dx between points 1 and 2 is given by

dm=ρmdx $dm = ρ_{m} d x$

(Eq. 8.11)

where ρ_m is the density per unit length x with the SI unit g cm^− 3 m^− 1.

Inserting Eq. 8.10 and Eq. 8.11 into Eq. 8.9 yields

d(F0tanβ−tanα)d(tanβ−tanα)dx==ρm∂2u∂t2dxρmF0∂2u∂t2 $\begin{array}{l} d (F_{0} tan β - tan α) & = & ρ_{m} \frac{\partial^{2} u}{\partial t^{2}} d x \\ \frac{d (tan β - tan α)}{d x} & = & \frac{ρ_{m}}{F_{0}} \frac{\partial^{2} u}{\partial t^{2}} \end{array}$

si48_e (Eq. 8.12)

for which we will need to express the change in the angles α and β as a function of x. The angle α is defined by the inclination of u at the location x = x₁ and is therefore given by

tanα=dux∣∣x=x1 ${tan α = \frac{d u}{x}|}_{x = x_{1}}$

whereas the angle β is given by the inclination of u at the location x = x₂, given by

tanβ=dux∣∣x=x2=u(x1)+dudx ${tan β = \frac{d u}{x}|}_{x = x_{2}} = \frac{u (x_{1}) + d u}{d x}$

si50_e

in which case we find

d(tanβ−tanα)dx=d(u(x1)+du−u(x1)dx)dx=ddudxdx=∂2u∂x2 $\frac{d (tan β - tan α)}{d x} = \frac{d (\frac{u (x_{1}) + d u - u (x_{1})}{d x})}{d x} = \frac{d \frac{d u}{d x}}{d x} = \frac{\partial^{2} u}{\partial x^{2}}$

si51_e (Eq. 8.13)

where we also account for the fact that u is a function both of x and the time t.

8.1.4.4 The Wave Equation

Finally inserting Eq. 8.13 into Eq. 8.12 yields

∂2u∂t2=F0ρm∂2u∂x2 $\frac{\partial^{2} u}{\partial t^{2}} = \frac{F_{0}}{ρ_{m}} \frac{\partial^{2} u}{\partial x^{2}}$

(Eq. 8.14)

∂2u∂t2=c2∂2u∂x2 $\frac{\partial^{2} u}{\partial t^{2}} = c^{2} \frac{\partial^{2} u}{\partial x^{2}}$

(Eq. 8.15)

which is referred to as the wave equation. Eq. 8.14 is the wave equation in physical notation, whereas Eq. 8.15 uses a standard notation denoting the constant expression F0ρm $\frac{F_{0}}{ρ_{m}}$ by the constant c². The constant is written such in order to denote that this value is always positive.

General Formulation

The most commonly used general notation for the wave equation is given by

∂2f∂t2=c2Δf $\frac{\partial^{2} f}{\partial t^{2}} = c^{2} Δ f$

where we have used the dependent variable f. Here we have used the Laplace operator, which we formally introduced in section 7.1.3.5.

One-Dimensional Wave Equation

In the case of the suspended string it would be sufficient to solve the wave equation in the one-dimensional form, which is given by

∂2f∂t2=c2∂2f∂x2 $\frac{\partial^{2} f}{\partial t^{2}} = c^{2} \frac{\partial^{2} f}{\partial x^{2}}$

(Eq. 8.16)

Three-Dimensional Wave Equation

The most general case of the wave equation is the three-dimensional wave equation, given by

∂2f∂t2=c2(∂2f∂x2+∂2f∂y2+∂2f∂z2) $\frac{\partial^{2} f}{\partial t^{2}} = c^{2} (\frac{\partial^{2} f}{\partial x^{2}} + \frac{\partial^{2} f}{\partial y^{2}} + \frac{\partial^{2} f}{\partial z^{2}})$

si57_e

8.1.5 Poisson Equation

The Poisson equation is an equation with which we will work a lot. The most important equation in fluid mechanics, the Navier-Stokes equation can, in many cases, be simplified to a Poisson equation. This equation is derived in section 11 and given by Eq. 11.40. If considering only stationary incompressible flow conditions, the Navier-Stokes equation can be simplified to

Δv→=1η(k→−∇→p) $Δ \vec{v} = \frac{1}{η} (\vec{k} - \vec{\nabla} p)$

si58_e

In the general formulation the Poisson equation is given by

Δf=g(x,y) $Δ f = g (x, y)$

(Eq. 8.17)

Two-Dimensional Poisson Equation

In two dimensions, the Poisson equation is given by

∂2f∂x2+∂2f∂y2=g(x,y) $\frac{\partial^{2} f}{\partial x^{2}} + \frac{\partial^{2} f}{\partial y^{2}} = g (x, y)$

si60_e

As we will see, finding solutions to Poisson equations is very important in fluid mechanics.

8.1.6 Laplace Equation

The Laplace equation is a special case of the Poisson equation. It is the homogeneous Poisson equation for which g (x, y) = 0 is assumed. It is given by

Δf=0 $Δ f = 0$

(Eq. 8.18)

Two-Dimensional Laplace Equation

In two dimensions it is given by

∂2f∂x2+∂2f∂y2=0 $\frac{\partial^{2} f}{\partial x^{2}} + \frac{\partial^{2} f}{\partial y^{2}} = 0$

si62_e

Gravitation

Newton’s law of gravity states that the attractive force between two masses m and M, which are separated by a distance r, is given as

F→=GmMr2e→r $\vec{F} = G \frac{m M}{r^{2}} {\vec{e}}_{r}$

(Eq. 8.19)

where G is the earth gravitational constant with the value 6.673 84 × 10^− 11 m³ kg^− 1 s^− 2. Please note that this force acts in radial direction, thus the unit vector e→r ${\vec{e}}_{r}$ . By dividing Eq. 8.19 by M we obtain

F→Ma→m==−Gmr2e→r−Gmr2e→r $\begin{array}{l} \frac{\vec{F}}{M} & = & - G \frac{m}{r^{2}} {\vec{e}}_{r} \\ {\vec{a}}_{m} & = & - G \frac{m}{r^{2}} {\vec{e}}_{r} \end{array}$

si65_e

which gives the gravitational acceleration a→m ${\vec{a}}_{m}$ created by the mass m. The sign has been introduced to correct the direction of the acceleration. We can understand a→m ${\vec{a}}_{m}$ as a field for which we can derive the total flux through the surface of a sphere with radius R given by

∮a→m⋅n→dA==∮a→mdA−∮GmR2dA $\begin{array}{l} \oint {\vec{a}}_{m} \cdot \vec{n} d A & = & \oint {\vec{a}}_{m} d A \\ = & - \oint G \frac{m}{R^{2}} d A \end{array}$

si68_e

where we can now insert the infinitesimal change in a sphere given by Eq. 3.107 to find

∮a→m⋅n→dA==−GmR28π∮rdr−4πGm $\begin{array}{l} \oint {\vec{a}}_{m} \cdot \vec{n} d A & = & - G \frac{m}{R^{2}} 8 π \oint r d r \\ = & - 4 π G m \end{array}$

si69_e (Eq. 8.20)

Eq. 8.20 is referred to as the first formulation of Newton’s law in Gauss’s notation. By using the Gauss’s theorem (see Eq. 7.10) we can now set

∮a→m⋅n→dA=∫∇→⋅a→mdV $\oint {\vec{a}}_{m} \cdot \vec{n} d A = \int \vec{\nabla} \cdot {\vec{a}}_{m} d V$

(Eq. 8.21)

By combining Eq. 8.20 and Eq. 8.21 we find

∫∇→⋅a→mdV=−4πGm $\int \vec{\nabla} \cdot {\vec{a}}_{m} d V = - 4 π G m$

(Eq. 8.22)

where we can exploit the fact that we can express the mass m by means of the density ρ as m=∫ρdV $m = \int ρ d V$ , which simplifies Eq. 8.22 to

∫∇→⋅a→mdV∇→⋅a→m==−4πG∫ρdV−4πGρ $\begin{array}{l} \int \vec{\nabla} \cdot {\vec{a}}_{m} d V & = & - 4 π G \int ρ d V \\ \vec{\nabla} \cdot {\vec{a}}_{m} & = & - 4 π G ρ \end{array}$

si73_e (Eq. 8.23)

which gives the second formulation of Newton’s law in Gauss’s notation. In the last step, we introduce the gravitational potential Ψ. We know from Eq. 7.4 that the divergence of the vector a→m ${\vec{a}}_{m}$ can be created by a vector potential if the latter is source- and sink-free. For gravitation this is the case. Therefore we can obtain a→m ${\vec{a}}_{m}$ from Ψ as

a→m=∇→⋅Ψ ${\vec{a}}_{m} = \vec{\nabla} \cdot Ψ$

(Eq. 8.24)

When combining Eq. 8.23 and Eq. 8.24 we obtain

ΔΨ==∇→⋅∇→⋅Ψ∇→⋅a→m=−4πGρ $\begin{array}{l} Δ Ψ & = & \vec{\nabla} \cdot \vec{\nabla} \cdot Ψ \\ = & \vec{\nabla} \cdot {\vec{a}}_{m} = - 4 π G ρ \end{array}$

si77_e (Eq. 8.25)

Eq. 8.25 tells us that the gravitational field is related to the mass distributions by a Laplace equation.

8.1.7 Helmholtz Equation

The Helmholtz equation is given by

Δf=−c2f $Δ f = - c^{2} f$

(Eq. 8.26)

For c = 0 it results in the Laplace equation. The Helmholtz equation often occurs during a substitution and separation of variables approach to a differential equation.

One-Dimensional Helmholtz Equation

In one dimension the Helmholtz equation is given by

∂2f∂x2=−c2f $\frac{\partial^{2} f}{\partial x^{2}} = - c^{2} f$

8.1.8 Classification

There is a common classification scheme used to distinguish PDEs. Referring to PDEs with the independent variables time t and space x the most general form of a PDE can be written as

A(t,x)∂2f∂x2+B(t,x)∂2f∂x∂t+C(t,x)∂2f∂t2+D(t,x)∂f∂x+E(t,x)∂f∂t+F(t,x)f=G(t,x) $A (t, x) \frac{\partial^{2} f}{\partial x^{2}} + B (t, x) \frac{\partial^{2} f}{\partial x \partial t} + C (t, x) \frac{\partial^{2} f}{\partial t^{2}} + D (t, x) \frac{\partial f}{\partial x} + E (t, x) \frac{\partial f}{\partial t} + F (t, x) f = G (t, x)$

where A (t, x), B (t, x), C (t, x), D (t, x), E (t, x), F (t, x), and G (t, x) are the coefficients and functions of the independent variables. The expression B² − 4 A C is referred to as the discriminant of the equation. Depending on the value of this discriminant PDEs are classified as being the following:

1. elliptic if B² − 4 A C < 0, e.g., the Laplace equation (A = C = 1)

2. parabolic if B² − 4 A C = 0, e.g., the heat equation (A > 0, B = C = 0)

3. hyperbolic if B² − 4 A C > 0, e.g., the wave equation (A = 1, C < 0, B = 0)

8.2 General Solutions to Selected Ordinary Differential Equations

8.2.1 Introduction

In this section we will look at solutions to some of the most commonly encountered ODEs. As it turns out, there are a few techniques that are used to solve these equations and we will briefly discuss them using suitable examples.

8.2.2 Separation of Variables and Partial Integration

Before venturing into more complicated differential equations, we will start with the most simple case for solving a differential equation, which we have already used in the example of the barrier lake. This method is generally referred to as separation of variable. Let us consider Eq. 3.6 that we used during the derivation of the solution to the volume equation of the barrier lake. This differential equation was given by

dVlakedt=c $\frac{d V_{lake}}{d t} = c$

which we will write in standard form

dfdt=c $\frac{d f}{d t} = c$

In the next step, we now separate the variables, i.e., the dependent from the independent variables, resulting in

df=cdt $d f = c d t$

(Eq. 8.27)

where we can now integrate both sides of the equation to result in

∫dff==c∫dtct+c1 $\begin{array}{l} \int d f & = & c \int d t \\ f & = & c t + c_{1} \end{array}$

si84_e

More often than expected, simple separation of the variable will allow solving a differential equation.

Higher-Order Integration

Sometimes, the notation of higher-order integrations during separation of variables causes a bit of trouble. Let us consider the following example:

d2fdx2=c $\frac{d^{2} f}{d x^{2}} = c$

(Eq. 8.28)

This notation sometimes causes difficulties in rewriting this equation the way we did for Eq. 8.27. However, if we rewrite Eq. 8.28 as

d2fdx2=ddx(dfdx)=c $\frac{d^{2} f}{d x^{2}} = \frac{d}{d x} (\frac{d f}{d x}) = c$

si86_e

it becomes much clearer to perform the separation, which would read

ddx(dfdx)∫d(dfdx)dfdx∫dff=====cc∫dxcx+c1∫(cx+c1)dxc2x2+c1x+c2 $\begin{array}{l} \frac{d}{d x} (\frac{d f}{d x}) & = & c \\ \int d (\frac{d f}{d x}) & = & c \int d x \\ \frac{d f}{d x} & = & c x + c_{1} \\ \int d f & = & \int (c x + c_{1}) d x \\ f & = & \frac{c}{2} x^{2} + c_{1} x + c_{2} \end{array}$

si87_e

8.2.3 Linear Differential Equations With Constant Coefficients

The second commonly encountered case of differential equations is differential equations that only contain linear terms. An ordinary inhomogeneous differential equation for the function f (x) with order n and constant coefficients is given by

dnfxn+==an−1dn−1fxn−1+⋯+a1dfx+a0xBm(x)bmxm+bm−1xm−1+⋯+b1x+b0 $\begin{array}{l} \frac{d^{n} f}{x^{n}} & + & a_{n - 1} \frac{d^{n - 1} f}{x^{n - 1}} + \dots + a_{1} \frac{d f}{x} + a_{0} x \\ = & B_{m} (x) \\ = & b_{m} x^{m} + b_{m - 1} x^{m - 1} + \dots + b_{1} x + b_{0} \end{array}$

si88_e (Eq. 8.29)

where B_n (x) is referred to as disturbance function. For B_m (x) = 0 the differential equation becomes a homogenous differential equation. aⁱ_… and bⁱ_… are constants. If a_n ≠ 1 the function is first divided by a_n.

8.2.3.1 Homogeneous and Particular Solution

The solution to an inhomogeneous differential equation is given by the sum of the homogeneous f_hom (x) and the inhomogeneous solution (referred to as the particular solution) f_part (x):

f(x)=fhom(x)+fpart(x) $f (x) = f_{hom} (x) + f_{part} (x)$

(Eq. 8.30)

This concept is very important as we will often proceed this way when solving differential equations. First we find the homogeneous solution; then in a second step, we derive the particular solution.

Homogeneous Solution

The approach for the solution of the homogeneous case for Eq. 8.29 is

fhom(x)=eλx $f_{hom} (x) = e^{λ x}$

(Eq. 8.31)

Inserting Eq. 8.31 into the homogeneous version of Eq. 8.29 results in

λneλx+an−1λn−1eλx+⋯+a1λneλx+a0eλxλn+an−1λn−1+⋯+a1λn+a0==00 $\begin{array}{l} λ^{n} e^{λ x} + a_{n - 1} λ^{n - 1} e^{λ x} + \dots + a_{1} λ^{n} e^{λ x} + a_{0} e^{λ x} & = & 0 \\ λ^{n} + a_{n - 1} λ^{n - 1} + \dots + a_{1} λ^{n} + a_{0} & = & 0 \end{array}$

si91_e (Eq. 8.32)

where Eq. 8.32 is referred to as the characteristic polynomial. It will have n solutions λn,λn−1,⋯,λ1,λ0 $λ_{n}, λ_{n - 1}, \dots, λ_{1}, λ_{0}$ which are referred to as roots of the characteristic polynomial. The sum of the individual solutions will form the solution to the homogeneous differential equation, and therefore

fhom(x)=cneλnx+cn−1eλn−1x+⋯+c1eλ1x+c0eλ0x $f_{hom} (x) = c_{n} e^{λ_{n} x} + c_{n - 1} e^{λ_{n - 1} x} + \dots + c_{1} e^{λ_{1} x} + c_{0} e^{λ_{0} x}$

where the (integration) constants cn,cn−1,⋯,c1,c0 $c_{n}, c_{n - 1}, \dots, c_{1}, c_{0}$ must be determined by the boundary conditions. Please note that λⁱ_… can be complex values, in which case the solutions may be written as sine and cosine functions using Euler’s formula, Eq. 3.38:

fhom(x)=cne(a+ib)x=cneaxeibx=cneax(cos(bx)+isin(bx)) $f_{hom} (x) = c_{n} e^{(a + i b) x} = c_{n} e^{ax} e^{i b x} = c_{n} e^{ax} (cos (b x) + i sin (b x))$

8.2.3.2 Inhomogeneous Solution

The approach for an inhomogeneous solution uses a polynomial of one order higher than the disturbance function. The approach in general is

fpart(x)=cm+1xm+1+⋯+c1x+c0 $f_{part} (x) = c_{m + 1} x^{m + 1} + \dots + c_{1} x + c_{0}$

The particular solution is then inserted into Eq. 8.29, and the coefficients are compared.

8.2.3.3 Characteristic Polynomials With Recurring Roots

In many cases, a characteristic polynomial may have recurring solutions. As an example, consider the differential equation

d2fx2=0 $\frac{d^{2} f}{x^{2}} = 0$

(Eq. 8.33)

which has the characteristic polynomial λ² = 0. In this case, there are two solutions: λ₁ = λ₂ = 0. Whenever a solution λ to the characteristic polynomial is recurring, the solution is a polynomial of type

fr(x)=(∑n=0r−1cnxn)eλx $f_{r} (x) = (\sum_{n = 0}^{r - 1} c_{n} x^{n}) e^{λ x}$

si98_e (Eq. 8.34)

where the parameter r is the number of recurrences of the solution. In the given case, the solution to the ODE Eq. 8.33 would be

f(x)=c0x0e0x+c1x1e0x=c0+c1x $f (x) = c_{0} x^{0} e^{0 x} + c_{1} x^{1} e^{0 x} = c_{0} + c_{1} x$

By setting this equation into Eq. 8.33 one can easily find this to be correct.

8.2.3.4 Characteristic Polynomials With Complex Roots

The characteristic polynomial may have complex solutions. As an example, consider the following equation:

d2fx2+3f=0 $\frac{d^{2} f}{x^{2}} + 3 f = 0$

(Eq. 8.35)

for which the characteristic polynomial is λ² = − 3, which yields the complex solutions λ1=3i−−√ $λ_{1} = \sqrt{3 i}$ and λ2=−3i−−√ $λ_{2} = - \sqrt{3 i}$ . In all cases of complex solutions we will always have a conjugated complex pair of solutions denoted by λ₁ = α + i β and λ₂ = α − i β. In this case our solution will be

f(x)==c0e(α+iβ)x+c1e(α−iβ)xc0eαxeiβx+c1eαxe−iβx $\begin{array}{r} f (x) & = & c_{0} e^{(α + i β) x} + c_{1} e^{(α - i β) x} \\ = & c_{0} e^{α x} e^{i β x} + c_{1} e^{α x} e^{- i β x} \end{array}$

si103_e

where we now apply Euler’s formula, Eq. 3.38 and Eq. 3.39, to convert them to sine and cosine functions, in which case we find

f(x)===c0eαx(cos(βx)+isin(βx))+c1eαx(cos(βx)−isin(βx))eαx(cos(βx)(c0+c1)+isin(βx)(c0−c1))eαx(c3cos(βx)+ic4sin(βx)) $\begin{array}{l} f (x) & = & c_{0} e^{α x} (cos (β x) + i sin (β x)) + c_{1} e^{α x} (cos (β x) - i sin (β x)) \\ = & e^{α x} (cos (β x) (c_{0} + c_{1}) + i sin (β x) (c_{0} - c_{1})) \\ = & e^{α x} (c_{3} cos (β x) + i c_{4} sin (β x)) \end{array}$

si104_e

Here it is important to remember that both the real and the imaginary parts of the function obtained are solutions to our differential equation. This is why we find in the case of a characteristic polynomial that has the complex roots ± β, the solution to be

f(x)=c0eαxcos(βx)+c1eαxsin(βx) $f (x) = c_{0} e^{α x} cos (β x) + c_{1} e^{α x} sin (β x)$

(Eq. 8.36)

Regarding our example ODE, Eq. 8.35, we find α = 0 and β=3–√ $β = \sqrt{3}$ , which yields

f(x)=c0cos(3–√x)+c1sin(3–√x) $f (x) = c_{0} cos (\sqrt{3} x) + c_{1} sin (\sqrt{3} x)$

as our solutions.

8.2.3.5 Characteristic Polynomials With Recurring Complex Roots

Combining both special cases we have just studied, whenever you come across a polynomial that has complex roots λ = α ± iβ, which recurs r times, Eq. 8.36 is simply multiplied by the terms used for recurring roots (see Eq. 8.34), yielding

f(x)=(∑n=0r−1cnxneαx)(c0eαxcos(βx)+c1eαxsin(βx)) $f (x) = (\sum_{n = 0}^{r - 1} c_{n} x^{n} e^{α x}) (c_{0} e^{α x} cos (β x) + c_{1} e^{α x} sin (β x))$

si108_e (Eq. 8.37)

8.2.3.6 Examples

In the following, we will look at some examples.

Second-Order Inhomogeneous Ordinary Differential Equation

The following differential equation is given:

d2fdx2−2dfdx=x−2 $\frac{d^{2} f}{d x^{2}} - 2 \frac{d f}{d x} = x - 2$

(Eq. 8.38)

boundary conditions

f(0)dfdx∣∣x=0==10 $\begin{array}{l} f (0) & = & 1 \\ {\frac{d f}{d x}|}_{x = 0} & = & 0 \end{array}$

si110_e

We first solve the homogeneous differential equation for which we find the characteristic polynomial to be

λ2−2λ=0 $λ^{2} - 2 λ = 0$

for which we find the two solutions λ₁ = 0 and λ₂ = 2. Therefore the homogeneous solution to Eq. 8.38 is

fhom(x)=c1+c2e2x $f_{hom} (x) = c_{1} + c_{2} e^{2 x}$

We now proceed to finding the partial solution. We need a polynomial of order 3 using the approach

fpart(x)=c3x3+c2x2+c1x+c0 $f_{part} (x) = c_{3} x^{3} + c_{2} x^{2} + c_{1} x + c_{0}$

(Eq. 8.39)

Inserting Eq. 8.39 into Eq. 8.38 yields

6c3+2c2−6c3x2−4c2x−2c1−6c3+(6c3−4c2)x+(2c2−2c1)==x−2x−2 $\begin{array}{l} 6 c_{3} + 2 c_{2} - 6 c_{3} x^{2} - 4 c_{2} x - 2 c_{1} & = & x - 2 \\ - 6 c_{3} + (6 c_{3} - 4 c_{2}) x + (2 c_{2} - 2 c_{1}) & = & x - 2 \end{array}$

si114_e

from which we find the following system of equations

−6c36c3−4c22c2−2c1===0→c3−1→c2−2→c1===0−1434 $\begin{array}{l} - 6 c_{3} & = & 0 \to c_{3} & = & 0 \\ 6 c_{3} - 4 c_{2} & = & - 1 \to c_{2} & = & - \frac{1}{4} \\ 2 c_{2} - 2 c_{1} & = & - 2 \to c_{1} & = & \frac{3}{4} \end{array}$

si115_e

Therefore the general solution to Eq. 8.38 is

f(x)==yhom(x)+ypart(x)c1+c2e2x−14x2+34x $\begin{array}{l} f (x) & = & y_{hom} (x) + y_{part} (x) \\ = & c_{1} + c_{2} e^{2 x} - \frac{1}{4} x^{2} + \frac{3}{4} x \end{array}$

si116_e

Using the boundary conditions, we can determine c2=−38 $c_{2} = - \frac{3}{8}$ and c1=2−c2=118 $c_{1} = 2 - c_{2} = \frac{11}{8}$ .

Solving the Differential Equation for Meniscus Formation at Liquid/Solid Interfaces

During the derivation of the shape of a meniscus forming due to surface tension, we have a very practical example of a linear ordinary differential equation with constant coefficients (see Eq. 21.9). The equation is given by

d2z(x)dx2−ρgγz(x)=0 $\frac{d^{2} z (x)}{d x^{2}} - \frac{ρ g}{γ} z (x) = 0$

si119_e

This is a homogenous second-order differential equation with constant coefficients. We proceed by finding the characteristic polynomial, which is given by

λ2−ρgγ=0λ=±ρgγ−−√ $\begin{array}{l} λ^{2} - \frac{ρ g}{γ} = 0 \\ λ = \pm \sqrt{\frac{ρ g}{γ}} \end{array}$

si120_e

Therefore we find the general solution to this differential equation to be

z(x)=c1eρgγ√x+c2e−ρgγ√x $z (x) = c_{1} e^{\sqrt{\frac{ρ g}{γ}} x} + c_{2} e^{- \sqrt{\frac{ρ g}{γ}} x}$

si121_e

8.2.4 Linear Differential Equations With Nonconstant Coefficients

After having studied methods for solving differential equations with constant coefficients we will now look at differential equations that take functions instead of constants. For linear differential equations of first order, these can be solved conveniently.

8.2.4.1 Introduction

We will look at differential equations of the type

a(x)df(x)dx+b(x)f(x)adfdx+bf==c(x)c $\begin{array}{l} a (x) \frac{d f (x)}{d x} + b (x) f (x) & = & c (x) \\ a \frac{d f}{d x} + b f & = & c \end{array}$

si122_e (Eq. 8.40)

with the coefficients a, b, and c being functions of x instead of constants: a (x), b (x), and c (x). These types of differential equations occur frequently, and they are very easy to solve. In the first step, we rewrite Eq. 8.40 so that the factor in front of the derivative becomes 1, resulting in

dfdx+baf=ca $\frac{d f}{d x} + \frac{b}{a} f = \frac{c}{a}$

(Eq. 8.41)

8.2.4.2 Integration Factor

In the next step we apply the product rule of differentiation (see Eq. 3.13). We assume that the right-hand side of the equation is the result of the differentiation of the product r (x) f (x) where we introduced a new coefficient r (x), which we still have to determine. This coefficient is referred to as an integration factor. We will see in a moment where this name comes from. If we differentiate r f with respect to x we are left with

d(rf)dx=rdfdx+drdxf $\frac{d (r f)}{d x} = r \frac{d f}{d x} + \frac{d r}{d x} f$

(Eq. 8.42)

If we multiply Eq. 8.41 with r we find

rdfdx+rbaf=rca $r \frac{d f}{d x} + r \frac{b}{a} f = r \frac{c}{a}$

(Eq. 8.43)

Now the left-hand side of Eq. 8.43 already looks very similar to the differentiated product in Eq. 8.42. In order to match exactly, we would need to ensure that

rba=drdx $r \frac{b}{a} = \frac{d r}{d x}$

by comparing the coefficients. From this equation, we can derive r (x) which is given by

drrlnrr(x)===badx∫badxe∫b(x)a(x)dx $\begin{array}{l} \frac{d r}{r} & = & \frac{b}{a} d x \\ ln r & = & \int \frac{b}{a} d x \\ r (x) & = & e^{\int \frac{b (x)}{a (x)} d x} \end{array}$

si127_e (Eq. 8.44)

We have now determined the integration factor r (x). Using this factor we can now rewrite Eq. 8.43 as

d(rf)dxd(rf)rff(x)====rcarcadx∫rcadx1r(x)∫r(x)c(x)a(x)dx $\begin{array}{l} \frac{d (r f)}{d x} & = & r \frac{c}{a} \\ d (r f) & = & r \frac{c}{a} d x \\ r f & = & \int r \frac{c}{a} d x \\ f (x) & = & \frac{1}{r (x)} \int r (x) \frac{c (x)}{a (x)} d x \end{array}$

si128_e (Eq. 8.45)

8.2.4.3 Solution

Eq. 8.44 and Eq. 8.45 give us the solution to the first-order differential equation with nonconstant coefficients, which is given by

f(x)=1e∫b(x)a(x)dx∫e∫b(x)a(x)dxc(x)a(x)dx $f (x) = \frac{1}{e^{\int \frac{b (x)}{a (x)} d x}} \int e^{\int \frac{b (x)}{a (x)} d x} \frac{c (x)}{a (x)} d x$

si129_e (Eq. 8.46)

We used the integration factor r (x) in order to combine the differential dfdx $\frac{d f}{d x}$ and the function f on the left-hand side of our original differential equation into one combined differential. This is the reason why these factors are referred to as integration factors. Integration factors allow expressing an (unknown) differential by an integral that depends on the coefficients.

8.2.4.4 Example

In the following we will solve the following differential equation using the method we just discussed:

2xdfdx−3f=9x3 $2 x \frac{d f}{d x} - 3 f = 9 x^{3}$

(Eq. 8.47)

where a (x) = 2 x, b (x) = − 3, and c (x) = 9 x³. According to Eq. 8.46 the solution to this ODE is

f(x)==1e∫−32dxx∫e∫−32dxx9x32xdx1e−32lnx+c0∫e−32lnx+c092x2dx $\begin{array}{l} f (x) & = & \frac{1}{e^{\int - \frac{3}{2} \frac{d x}{x}}} \int e^{\int - \frac{3}{2} \frac{d x}{x}} \frac{9 x^{3}}{2 x} d x \\ = & \frac{1}{e^{- \frac{3}{2} ln x + c_{0}}} \int e^{- \frac{3}{2} ln x + c_{0}} \frac{9}{2} x^{2} d x \end{array}$

si132_e

Here, we can chose any value for c₀. For convenience we will set it to 0, in which case we find

f(x)======1e−32lnx∫e−32lnx92x2dx1x−32∫x−3292x2dxx32∫x−3292x2dxx3292∫x12dxx32(3x32+c1)3x3+c1x32 $\begin{array}{l} f (x) & = & \frac{1}{e^{- \frac{3}{2} ln x}} \int e^{- \frac{3}{2} ln x} \frac{9}{2} x^{2} d x \\ = & \frac{1}{x^{- \frac{3}{2}}} \int x^{- \frac{3}{2}} \frac{9}{2} x^{2} d x \\ = & x^{\frac{3}{2}} \int x^{- \frac{3}{2}} \frac{9}{2} x^{2} d x \\ = & x^{\frac{3}{2}} \frac{9}{2} \int x^{\frac{1}{2}} d x \\ = & x^{\frac{3}{2}} (3 x^{\frac{3}{2}} + c_{1}) \\ = & 3 x^{3} + c_{1} x^{\frac{3}{2}} \end{array}$

si133_e (Eq. 8.48)

Let us quickly verify that this is in fact a solution to our differential equation. Using Eq. 8.48 in Eq. 8.47 yields

2x(9x2+32c1x12)−3(3x3+c1x32)18x3+3c1x32−9x3−3c1x32==9x39x3 $\begin{array}{l} 2 x (9 x^{2} + \frac{3}{2} c_{1} x^{\frac{1}{2}}) - 3 (3 x^{3} + c_{1} x^{\frac{3}{2}}) & = & 9 x^{3} \\ 18 x^{3} + 3 c_{1} x^{\frac{3}{2}} - 9 x^{3} - 3 c_{1} x^{\frac{3}{2}} & = & 9 x^{3} \end{array}$

si134_e

which is obviously correct.

8.2.5 Laplace Transform for Solving the Harmonic Oscillation ODE

8.2.5.1 Introduction

In section 5.2 we have introduced the Laplace transform, which is an extremely useful tool to solving complex ODEs. In this section we will use the Laplace transform for solving the ODE of the damped harmonic oscillation. The details of this problem are introduced in section 27.3.1. The resulting equation is a second-order homogeneous ODE given by Eq. 27.32 as

md2sdt2(t)+kv(t)+Ds(t)=0 $m \frac{d^{2} s}{d t^{2}} (t) + k v (t) + D s (t) = 0$

(Eq. 8.49)

The independent variable is the time t, and the dependent variable is the position or path s, which is a function of time. Obviously, the velocity v is the first derivative of s with respect to time, and the acceleration is the second derivative with respect to time. The constants are the oscillating mass m, the damping constant k, and the spring constant D.

As in section 27.3.1 we assume the following initial values: s (t = 0) = S and v (t = 0) = 0. At the beginning of the experiment, the mass is at rest but displaced to the maximum amplitude S from which it is left to oscillate.

We will now derive the solution to this equation using the Laplace transform. The analytical solutions are derived in section 27.3.1 using a different approach. We will compare the results found with the Laplace transform with these equations.

8.2.5.2 Laplace Transform of the ODE

We first need to look up the Laplace transforms of the individual terms from Tab. 5.2, from which we find Eq. 5.44 that gives the Laplace transforms of derivatives in time domain as

d2sdt2(t)dsdt(t)s(t)===f2S(f)−dsdt(0)−fs(0)=f2S(f)−fSfS(f)−s(0)=fS(f)−SS(f) $\begin{array}{l} \frac{d^{2} s}{d t^{2}} (t) & = & f^{2} S (f) - \frac{d s}{d t} (0) - f s (0) = f^{2} S (f) - f S \\ \frac{d s}{d t} (t) & = & f S (f) - s (0) = f S (f) - S \\ s (t) & = & S (f) \end{array}$

si136_e

If we insert these equations into Eq. 8.49 we find

m(f2S(f)−fS)+k(fS(f)−S)+DS(f)=0 $m (f^{2} S (f) - f S) + k (f S (f) - S) + D S (f) = 0$

(Eq. 8.50)

which we have to solve for S (f), finding

mf2S(f)+(fk+D)S(f)−S(fm+k)S(f)(mf2+kf+D)S(f)===0S(fm+k)S(fm+k)mf2+kf+D $\begin{array}{l} m f^{2} S (f) + (f k + D) S (f) - S (f m + k) & = & 0 \\ S (f) (m f^{2} + k f + D) & = & S (f m + k) \\ S (f) & = & \frac{S (f m + k)}{m f^{2} + k f + D} \end{array}$

si138_e (Eq. 8.51)

We have now found the solution to the ODE of Eq. 8.49. Unfortunately, this solution is still in the frequency domain and we first have to reconvert it to the time domain.

8.2.5.3 Method of Partial Fraction Expansion

As we have already discussed in section 5.2, the inversion of the Laplace transform in generally the problematic step. If we look at Tab. 5.2 we can see that we do not find an obvious transform pair that would match Eq. 8.51.

However, we can see that there are numerous transform pairs with have terms of the type (s − a) in the denominator. Here, a is a constant. It is a root of the denominator, i.e., a constant that makes the denominator equal to 0. These result from a process which is referred to as the method of partial fraction expansion.

If we look at Tab. 5.2 we can see that we have numerous transform pairs where the function in the frequency domain is of the form of simple fraction terms such as cf−a $\frac{c}{f - a}$ . We also know from section 5.2.3 that we can very conveniently retransform sums of transformed functions from the frequency domain to the time domain.

Finding Roots

The method of partial fractions allows us to transform a rather complicated looking polynomial fraction such as the one we derived in Eq. 8.51 to a product of several simple fractions. As an example, take

2f−8f2−5f+6=4f−2−2f−3 $\frac{2 f - 8}{f^{2} - 5 f + 6} = \frac{4}{f - 2} - \frac{2}{f - 3}$

(Eq. 8.52)

You can check this by finding a common denominator of the two fractions, which is done by multiplying the denominator of one fraction with the denominator of the other, respectively. Mathematically speaking, we have decomposed the polynomial f² + 5f + 6 into a product of its roots, which are λ₁ = 2 and λ₁ = 3. These can be calculated using the formula for finding the roots of the quadratic equation ax² + bx + c = 0, given by

x1,2=−b±b2−4ac√2a $x_{1, 2} = \frac{- b \pm \sqrt{b^{2} - 4 a c}}{2 a}$

si141_e (Eq. 8.53)

Approaches for Partial Fractions

Obviously, the roots do not necessarily have to be real numbers, they may also be complex. Depending on the order of the denominator and the type of roots found, several common approaches for partial fractions are used. The most common ones are listed in Tab. 8.1. As you can see, the approaches are lacking the constants A_i, B, which need to be determined. We will carry out such a partial fraction expansion for Eq. 8.51.

Tab. 8.1

Common approaches for partial fraction expansion that is required, e.g., for the inversion of the Laplace transform. The polynomial to decompose is P(f)Q(f) $\frac{P (f)}{Q (f)}$ , where the degree of the denominator is bigger than the degree of the numerator. Please note that the denominator polynomial must be normalized, i.e., the coefficient of the first term must be 1

Polynomial and expansion approach	Roots of denominator	Comment
P(f)Q(f)=P(f)f2+af+b=Af−λ1+Bf−λ2 $\frac{P (f)}{Q (f)} = \frac{P (f)}{f^{2} + a f + b} = \frac{A}{f - λ_{1}} + \frac{B}{f - λ_{2}}$	λ₁, λ₂	denominator with two different noncomplex roots
P(f)Q(f)=P(f)f2+af+b=A(f−λ)2 $\frac{P (f)}{Q (f)} = \frac{P (f)}{f^{2} + a f + b} = \frac{A}{{(f - λ)}^{2}}$	λ₁, λ₂	denominator with one recurring noncomplex roots
P(f)Q(f)=P(f)∑n−1i=0aifi=∑n−1i=0Aif−λi $\frac{P (f)}{Q (f)} = \frac{P (f)}{\sum_{i = 0}^{n - 1} a_{i} f^{i}} = \sum_{i = 0}^{n - 1} \frac{A_{i}}{f - λ_{i}}$	λ_i	denominator with n different real roots
P(f)Q(f)=P(f)f2+af+b=A(s−R{λ})+B(f−R{λ})2−I{β}2 $\frac{P (f)}{Q (f)} = \frac{P (f)}{f^{2} + a f + b} = \frac{A (s - ℜ \{λ\}) + B}{{(f - ℜ \{λ\})}^{2} - ℑ {\{β\}}^{2}}$	λ₁, λ₂, λ1¯¯¯¯ $\bar{λ_{1}}$ , λ2¯¯¯¯ $\bar{λ_{2}}$	denominator with two complex-conjugated roots

Performing Partial Fraction Expansion

In many cases, partial fraction expansion allows us to transform a fraction with polynomial denominator to a sum of very simple polynomials with constant denominators. As an example, take the right-hand side of Eq. 8.52. This is significantly easier to retransform to the time domain than the left-hand side. The process is very simple. If we consider Eq. 8.52 we find the following roots:

λ1λ2==5+25−24√2=35−25−24√2=2 $\begin{array}{l} λ_{1} & = & \frac{5 + \sqrt{25 - 24}}{2} = 3 \\ λ_{2} & = & \frac{5 - \sqrt{25 - 24}}{2} = 2 \end{array}$

si142_e

which we obtain by applying Eq. 8.53. We could then rewrite Eq. 8.52 as

2f−8f2−5f+6=Af−λ1+Bf−λ2=Af−2+Bf−3 $\frac{2 f - 8}{f^{2} - 5 f + 6} = \frac{A}{f - λ_{1}} + \frac{B}{f - λ_{2}} = \frac{A}{f - 2} + \frac{B}{f - 3}$

(Eq. 8.54)

where A and B are constants. We determine these by comparison with the left-hand side of Eq. 8.54. For this we reexpand the right-hand side of Eq. 8.54 on a common denominator, which we do by multiplying the first term with (f − 3) and the second term with (f − 2), yielding

2f−8f2−5f+6==2f−8(f−2)(f−3)=A(f−3)+B(f−2)(f−2)(f−3)f(A+B)−3A−2B(f−2)(f−3) $\begin{array}{l} \frac{2 f - 8}{f^{2} - 5 f + 6} & = & \frac{2 f - 8}{(f - 2) (f - 3)} = \frac{A (f - 3) + B (f - 2)}{(f - 2) (f - 3)} \\ = & \frac{f (A + B) - 3 A - 2 B}{(f - 2) (f - 3)} \end{array}$

si144_e (Eq. 8.55)

We now simply have to compare coefficients to find that A+B=!2→B=2−A $A + B \overset{!}{=} 2 \to B = 2 - A$ and −3A−2B=−3A−(4−2A)=!−8 $- 3 A - 2 B = - 3 A - (4 - 2 A) \overset{!}{=} - 8$ , which is satisfied for A = 4; in this case B = − 2, which is exactly the right-hand side of Eq. 8.52.

When to Use Partial Fraction Expansion

Partial fraction expansion allowed us to decompose the rather complicated polynomial to two simple polynomials that can be conveniently retransformed using Tab. 5.2. To be exact, we would use Eq. 5.29 with which we can write down the retransformed as

h(t)=4e2t−2e3t $h (t) = 4 e^{2 t} - 2 e^{3 t}$

In many cases, partial fraction expansion significantly simplifies the retransformation step. However, as we will see in a moment, if a comprehensive look-up table is available, it may not be required to perform the full partial fraction expansion since intermediary polynomials may already be available for look-up.

8.2.5.4 Performing the Retransform

As stated, we require the roots for the denominator of Eq. 8.51, i.e., the values f_i for which the denominator becomes zero:

0f1,2==mf2+kf+D−k±k2−4mD√2m $\begin{array}{l} 0 & = & m f^{2} + k f + D \\ f_{1, 2} & = & \frac{- k \pm \sqrt{k^{2} - 4 m D}}{2 m} \end{array}$

si148_e

where we applied Eq. 8.53 for finding the roots. If you compare this calculation with the derivation of the analytical solution in section 27.3.1 you will see that we derived the same characteristic polynomial as given by Eq. 27.33, which obviously has the same roots. We have to make a case distinction here because depending on the value in the square root we either get complex roots or not. In both cases we have to use different approaches (see Tab. 8.1)

Overdamped Oscillation: k² − 4 m D > 0

The first case we will consider is the overdamped oscillation. It occurs when the momentum of the moving mass is consumed by the damping in the system and thus the oscillation will decay. At this point, we will insert numbers for the variables in order to simplify the calculations a little. We set D = 1, k = 3, and m = 1, which are the values used in section 27.3.1 as well as in listing 27.7. We also set S = 1. In this case, we obtain the roots

f1f2==−k+k2−4mD√2m=−3+5√2−k−k2−4mD√2m=−3−5√2 $\begin{array}{l} f_{1} & = & \frac{- k + \sqrt{k^{2} - 4 m D}}{2 m} = \frac{- 3 + \sqrt{5}}{2} \\ f_{2} & = & \frac{- k - \sqrt{k^{2} - 4 m D}}{2 m} = \frac{- 3 - \sqrt{5}}{2} \end{array}$

si149_e

in which case we can decompose the denominator polynomial of Eq. 8.51 according to

S(fm+k)mf2+kf+D=1(f+3)(f−f1)(f−f2)=f(f−f1)(f−f2)+3f(f−f1)(f−f2) $\frac{S (f m + k)}{m f^{2} + k f + D} = \frac{1 (f + 3)}{(f - f_{1}) (f - f_{2})} = \frac{f}{(f - f_{1}) (f - f_{2})} + 3 \frac{f}{(f - f_{1}) (f - f_{2})}$

si150_e

Obviously, we could further decompose the term into full polynomials but this is not required at this point. If we look at Tab. 5.2 we can see that we have suitable transform pairs we can use here (see Eq. 5.33 and Eq. 5.34). Using these equations, we find the solution

s(t)f1−f2s(t)f1+3f2+3s(t)======f1ef1t−f2ef2tf1−f2+3ef1t−ef2tf1−f2−3+5√2+3+5√2=5–√15√((f1+3)ef1t−(f2+3)ef2t)−3+5√2+3=3+5√2−3−5√2+3=3−5√25√10((3+5–√)e−3+5√2t+(−3+5–√)e−3−5√2t) $\begin{array}{l} s (t) & = & \frac{f_{1} e^{f_{1} t} - f_{2} e^{f_{2} t}}{f_{1} - f_{2}} + 3 \frac{e^{f_{1} t} - e^{f_{2} t}}{f_{1} - f_{2}} \\ f_{1} - f_{2} & = & \frac{- 3 + \sqrt{5}}{2} + \frac{3 + \sqrt{5}}{2} = \sqrt{5} \\ s (t) & = & \frac{1}{\sqrt{5}} ((f_{1} + 3) e^{f_{1} t} - (f_{2} + 3) e^{f_{2} t}) \\ f_{1} + 3 & = & \frac{- 3 + \sqrt{5}}{2} + 3 = \frac{3 + \sqrt{5}}{2} \\ f_{2} + 3 & = & \frac{- 3 - \sqrt{5}}{2} + 3 = \frac{3 - \sqrt{5}}{2} \\ s (t) & = & \frac{\sqrt{5}}{10} ((3 + \sqrt{5}) e^{\frac{- 3 + \sqrt{5}}{2} t} + (- 3 + \sqrt{5}) e^{\frac{- 3 - \sqrt{5}}{2} t}) \end{array}$

si151_e (Eq. 8.56)

where we have found the retransformed solution in the time domain. If you compare Eq. 8.55 with the analytical solution (see Eq. 27.38) you will see that this is in fact the correct result. As you can see, with a good look-up table the conversion of the Laplace transformed solution into the frequency domain is simple.

Critically-Damped Oscillation: k²− 4 m D = 0

The second case is referred to as critically damped oscillation. It is the case in which the oscillation decays most quickly. We chose the following values for this case: D = 1, k = 2, and m = 1. We now have are curring root that is given by

f1=f2=−k2m=−22=−1 $f_{1} = f_{2} = \frac{- k}{2 m} = - \frac{2}{2} = - 1$

We can therefore decompose the denominator polynomial of Eq. 8.51 according to

S(fm+k)mf2+kf+D=1(f+2)(f−f1)2=f(f−f1)2+21(f−f1)2 $\frac{S (f m + k)}{m f^{2} + k f + D} = \frac{1 (f + 2)}{{(f - f_{1})}^{2}} = \frac{f}{{(f - f_{1})}^{2}} + 2 \frac{1}{{(f - f_{1})}^{2}}$

si153_e

Referring to Tab. 5.2 we can find Eq. 5.31 and Eq. 5.30, which gives the following retransformed solutions:

s(t)===(1+f1t)ef1t+2tef1t(1−t)e−t+2te−t(1+t)e−t $\begin{array}{l} s (t) & = & (1 + f_{1} t) e^{f_{1} t} + 2 t e^{f_{1} t} \\ = & (1 - t) e^{- t} + 2 t e^{- t} \\ = & (1 + t) e^{- t} \end{array}$

si154_e

which is the solution in the time domain. If you compare this equation with the analytical solution derived in Eq. 27.40 you will see that this is the correct result.

Underdamped Oscillation: k² − 4 m D < 0

The last case we will study is called underdamped oscillation, in which the damping does not consume the momentum of the moving mass fast enough to hinder the system from oscillating. We chose the following values for this case: D = 1, k = 1, and m = 5. We find Eq. 8.51 as

S(f)===S(fm+k)mf2+kf+D5f+15f2+f+1f+15f2+15f+15 $\begin{array}{l} S (f) & = & \frac{S (f m + k)}{m f^{2} + k f + D} \\ = & \frac{5 f + 1}{5 f^{2} + f + 1} \\ = & \frac{f + \frac{1}{5}}{f^{2} + \frac{1}{5} f + \frac{1}{5}} \end{array}$

si155_e

The reason for the last rearrangement is that, as indicated in Tab. 8.1, we have to scale the equation such that the factor of the f² term in the denominator polynomial is 1. We now find the following conjugated complex roots:

f1f2==−110+−19√10=−110+i−19√10−110−−19√10=f1¯¯¯=−110−i−19√10 $\begin{array}{l} f_{1} & = & - \frac{1}{10} + \frac{\sqrt{- 19}}{10} = - \frac{1}{10} + i \frac{\sqrt{- 19}}{10} \\ f_{2} & = & - \frac{1}{10} - \frac{\sqrt{- 19}}{10} = \bar{f_{1}} = - \frac{1}{10} - i \frac{\sqrt{- 19}}{10} \end{array}$

si156_e

Referring to Tab. 8.1 we will use the following approach for the decomposition:

f+15f2+15f+15=!A(f+110)+B(f+110)2+19100 $\frac{f + \frac{1}{5}}{f^{2} + \frac{1}{5} f + \frac{1}{5}} \overset{!}{=} \frac{A (f + \frac{1}{10}) + B}{{(f + \frac{1}{10})}^{2} + \frac{19}{100}}$

si157_e (Eq. 8.57)

It takes a moment to see that the denominators are actually equal. We will shortly perform the multiplication of the right-hand side of Eq. 8.57, for which we find

(f+110)2+19100=f2+15f+1100+19100=f2+15f+15 ${(f + \frac{1}{10})}^{2} + \frac{19}{100} = f^{2} + \frac{1}{5} f + \frac{1}{100} + \frac{19}{100} = f^{2} + \frac{1}{5} f + \frac{1}{5}$

We now simply have to determine the missing constants A and B by comparing coefficients in the numerator of Eq. 8.57 for which we find

f+15=!A(f+110)+B $f + \frac{1}{5} \overset{!}{=} A (f + \frac{1}{10}) + B$

from which we find A = 1 and B=110 $B = \frac{1}{10}$ , which completes the partial fraction decomposition; for the decomposition we find

f+15f2+15f+15==(f+110)+110(f+110)2+19100f+110(f+110)2+19100+110(f+110)2+19100 $\begin{array}{l} \frac{f + \frac{1}{5}}{f^{2} + \frac{1}{5} f + \frac{1}{5}} & = & \frac{(f + \frac{1}{10}) + \frac{1}{10}}{{(f + \frac{1}{10})}^{2} + \frac{19}{100}} \\ = & \frac{f + \frac{1}{10}}{{(f + \frac{1}{10})}^{2} + \frac{19}{100}} + \frac{\frac{1}{10}}{{(f + \frac{1}{10})}^{2} + \frac{19}{100}} \end{array}$

si161_e (Eq. 8.58)

From Tab. 5.2 we use Eq. 5.41 and Eq. 5.42 for retransforming Eq. 8.58 into the time domain for which we find

s(t)===e−t10cos(19√10t)+119√e−t10sin(19√10t)e−t10(cos(19√10t)+119√sin(19√10t))e−t10(19√19sin(19√10t)+cos(19√10t)) $\begin{array}{l} s (t) & = & e^{- \frac{t}{10}} cos (\frac{\sqrt{19}}{10} t) + \frac{1}{\sqrt{19}} e^{- \frac{t}{10}} sin (\frac{\sqrt{19}}{10} t) \\ = & e^{- \frac{t}{10}} (cos (\frac{\sqrt{19}}{10} t) + \frac{1}{\sqrt{19}} sin (\frac{\sqrt{19}}{10} t)) \\ = & e^{- \frac{t}{10}} (\frac{\sqrt{19}}{19} sin (\frac{\sqrt{19}}{10} t) + cos (\frac{\sqrt{19}}{10} t)) \end{array}$

si162_e

which is the solution in the time domain. If you compare this equation with the analytical solution derived in Eq. 27.42 you will see that this is the correct result.

Visualization

At this point, we will not bother visualizing our results because we will do this when studying numerical methods for which we will use this example as well. Please refer to section 27.3.1 for details on the solution and to Fig. 27.12 for a visualization of the results.

8.2.5.5 Summary

In this section we have briefly studied the Laplace transform as a method for solving ODEs in the time domain using damped oscillation as an example. As we have seen, the transformation is simple and we quickly come up with a solution of the ODE in the frequency domain. The critical aspect about the Laplace transform is the retransformation. We have used a rather complicated example which gave rise to three distinct cases each of which produced a different solution in the frequency domain. Retransforming these solutions into the time domain is possible using adequate partial fraction expansions, which allow us to derive polynomial fractions for which we can look up the corresponding equations in the time domain from reference tables.

The Laplace transform is a very important method to solve ODEs on semi-infinite domains, i.e., whenever the independent variable is bound from 0 to ∞. In this example, the time t was the independent variable and we were interested in solutions to the ODE for values of t ≥ 0. Even though we usually refer to the original domain as the time domain we could also refer to it as the original domain. In many cases, the independent variable may not be time but, e.g., a space coordinate such as x. Then the frequency domain would better be referred to as the image-domain.

Besides ODEs the Laplace transform is a very important tool to solve PDEs. We will study an example where the Laplace transform is used to solve a PDE in section 8.3.7 for solving.

8.3 General Solutions to Selected Partial Differential Equations

8.3.1 Introduction

Having introduced methods to solve ODEs we will turn to methods for solving PDEs. As it turns out, PDEs are significantly more challenging to solve and there is only a very limited number of general solutions that can be obtained for PDEs. Luckily, some of the most relevant PDEs that we will come across can be solved. In the following we will look at a selection of these equations and discuss the approaches used to solve them. We will see most of these equations again in later chapters and the methods for solving them form the basis of some of the concepts that we will employ later for solving the fundamental equations in fluid mechanics.

Rationale for Solving PDEs

With very few exceptions, the general method for solving a PDE is to first turn it into a system of several independent ODEs. For each independent variable, we will require a separate ODE. The main trick in solving PDEs lies in the method for splitting it into ODEs. Once these ODEs are obtained, the task is straightforward and we may apply any suitable method discussed in section 8.2 to solve the ODEs. After the individual solutions are found, we will reassemble these equations into one single equation which will (if we made no mistakes) be a solution to the PDE.

8.3.2 The Method of Characteristics: The Transport Equation

In this section we will solve the transport equation (see section 8.1.1).

8.3.2.1 Characteristic Lines

The transport equation is very easy to solve, which is why we will quickly do it. We begin by forming the total differential (see section 3.1.5.1) of the function f, which is given by

df=∂f∂tdt+∂f∂xdx $d f = \frac{\partial f}{\partial t} d t + \frac{\partial f}{\partial x} d x$

(Eq. 8.59)

As x is a function of the time, we can rewrite dx as

dx=∂x∂tdt $d x = \frac{\partial x}{\partial t} d t$

in which case we find from Eq. 8.59

dfdfdt==∂f∂tdt+∂f∂x∂x∂tdt∂f∂t+∂f∂x∂x∂t $\begin{array}{l} d f & = & \frac{\partial f}{\partial t} d t + \frac{\partial f}{\partial x} \frac{\partial x}{\partial t} d t \\ \frac{d f}{d t} & = & \frac{\partial f}{\partial t} + \frac{\partial f}{\partial x} \frac{\partial x}{\partial t} \end{array}$

si165_e (Eq. 8.60)

If we compare Eq. 8.59 with Eq. 8.1 we first find that the whole equation must be equal to 0. Second, we can simply derive the following equations by comparison of the coefficients:

c=∂x∂t=dxdt→x(t)=ct+c0 $c = \frac{\partial x}{\partial t} = \frac{d x}{d t} \to x (t) = c t + c_{0}$

(Eq. 8.61)

a=−1fdfdt→f(t)=e−at+c1=c2e−at $a = - \frac{1}{f} \frac{d f}{d t} \to f (t) = e^{- a t + c_{1}} = c_{2} e^{- at}$

(Eq. 8.62)

Eq. 8.61 gives a very important information. It states that x (t) − ct = c₀. If you recall the second version of looking at the traveling water wave, you will see that this is exactly what this equation tells you: As long as you move your eyes along with the velocity c you will not see a relative shift in the height of the wave. The point you chose to look at is given by c₀. There obviously is an infinite number of points you could chose to look at. Any point would do as long as you move your eyes at the correct velocity. The lines along which you move your eyes are referred to as characteristic lines, which is why this method of solving the differential equation is sometimes referred to as the method of characteristics.

8.3.2.2 Solution

We can now use Eq. 8.62 to formulate f (x, t) as

f(x,t)=c2e−at $f (x, t) = c_{2} e^{- a t}$

(Eq. 8.63)

for which we have to find c₂ by using a boundary condition. Let this boundary condition be at t = 0, for which we find

f(x(t),0)=f(ct+c0,0)=f(c⋅0+c0,0)=f(c0,0)=c2e−a⋅0=c2 $f (x (t), 0) = f (c t + c_{0}, 0) = f (c \cdot 0 + c_{0}, 0) = f (c_{0}, 0) = c_{2} e^{- a \cdot 0} = c_{2}$

where we have inserted Eq. 8.61 to convert x (t). If we solve Eq. 8.61 for c₀ we can write

f(c0,0)=f(x−ct,0):=F(x−ct) $f (c_{0}, 0) = f (x - c t, 0) : = F (x - c t)$

(Eq. 8.64)

where we have introduce, the function F (x − ct), which is any chosen function. If we now insert Eq. 8.64 in Eq. 8.62 and Eq. 8.63 we find

f(x,t)=F(x−ct)e−at $f (x, t) = F (x - c t) e^{- a t}$

(Eq. 8.65)

which is the solution to the transport equation. The transport equation can be solved by any function F (x − ct) that is multiplied with an exponential term that describes the decay over time. As discussed, the decay is a function only of time, and it is independent of x.

On one hand it seems strange that we may really use "any" function as solution. However, in principle any physical shape or contour can be transported, no matter which function. This may be a traveling rectangular pulse, a sine wave, or any other function we may come up with. The only condition they have to fulfill is that they are defined as functions of the independent variable term x − ct.

8.3.2.3 Visualization

It may be easier to understand the implications of this solution when looking at some practical examples. In Fig. 8.2a you can see a sine wave that is traveling at constant speed c = 1. This may be a water wave from our introductory example. It is defined as

f08-02-9781455731411 — Fig. 8.2 Sine wave function as a solution to the transport equation. The function travels at the speed c = 1 toward the right. The function is shown for the time points t = 0, t = 1, t = 3, t = 5, and t = 10. a) Undamped case. The function does not decay over time. b) Damped case assuming a damping factor a = 0.01. As you can see, the function decays over time.

f1(x,t)={sin(x−ct)00≤x−ct≤2πotherwise $f_{1} (x, t) = \{\begin{cases} sin (x - c t) & 0 \leq x - c t \leq 2 π \\ 0 & otherwise \end{cases}$

si172_e

As can be seen from the plot, the water wave travels to the right, e.g., it moves in space without actually changing its shape. It is "transported" along: A classical example of a transport phenomena. For Fig. 8.2a the damping factor was set to a = 0, which is why you do not see a decrease in the function over time. In Fig. 8.2b you can see the same function, but this time it is damped using a damping factor a = 0.01. Obviously, the function still moves to the right but it slowly decays over time. This is a somewhat more realistic scenario since we would expect to see damping, also in the introductory example of the water waves.

Just to give you an idea of how a different solution function to the transport equation would look, Fig. 8.3 shows the pulse function transported without damping (see Fig. 8.3a) and with damping (see Fig. 8.3b). This function is defined as

f08-03-9781455731411 — Fig. 8.3 Pulse function as a solution to the transport equation. The function travels at the speed c = 1 toward the right. The function is shown for the time points t = 0, t = 1, t = 3, t = 5, and t = 10. a) Undamped case. The function does not decay over time. b) Damped case assuming a damping factor a = 0.01. As you can see, the function decays over time.

f2(x,t)={100≤x−ct≤0.75otherwise $f_{2} (x, t) = \{\begin{cases} 1 & 0 \leq x - c t \leq 0.75 \\ 0 & otherwise \end{cases}$

si173_e

As you can see, this function is a solution to the transport equation.

Wave Collision

A phenomenon we often observe is two waves that run in opposite directions and eventually "collide." They temporarily fuse into a single wave, pass each other, and travel onwards. This effect can be visualized quite conveniently using the sine wave we have used in Fig. 8.2. For this we simply defined two waves traveling in opposite directions, i.e., for one wave we assume c = 1 while for the other we assume c = − 1. The effect of the waves "passing each other" is shown in Fig. 8.4.

f08-04-9781455731411 — Fig. 8.4 Two waves "colliding." The waves travel in the opposite direction with one wave at c = 1 while the other wave is at c = − 1. The second wave is offset to start at the right side.

8.3.3 Substitution and Separation of Variables: The One-Dimensional Heat Equation

The next method we will look into is often also referred to as "separation of variables" because it actually uses this method. However, it uses a very important trick before performing separation of variables, which is why we will refer to it as substitution and separation of variables. We will discuss this method using a practical example: The one-dimensional heat equation (see Eq. 8.3).

∂f∂t=c2∂2f∂x2 $\frac{\partial f}{\partial t} = c^{2} \frac{\partial^{2} f}{\partial x^{2}}$

(Eq. 8.66)

As we can see, f depends on both of our independent variables x and t.

8.3.3.1 Substitution

We will assume that f is actually a function that consists of two functions X and T. Each of these functions is assumed to only depend on one of the independent variables. We assume

f(x,t)=X(x)T(t) $f (x, t) = X (x) T (t)$

(Eq. 8.67)

which allows us to find the two partial derivatives in Eq. 8.66 as

∂f∂t=X∂T∂t+∂X∂tT=X∂T∂t $\frac{\partial f}{\partial t} = X \frac{\partial T}{\partial t} + \frac{\partial X}{\partial t} T = X \frac{\partial T}{\partial t}$

(Eq. 8.68)

∂2f∂x2=∂∂x(T∂X∂x+X∂T∂x)=∂∂x(T∂X∂x)=T∂2X∂x2+∂X∂x∂T∂x=T∂2X∂x2 $\frac{\partial^{2} f}{\partial x^{2}} = \frac{\partial}{\partial x} (T \frac{\partial X}{\partial x} + X \frac{\partial T}{\partial x}) = \frac{\partial}{\partial x} (T \frac{\partial X}{\partial x}) = T \frac{\partial^{2} X}{\partial x^{2}} + \frac{\partial X}{\partial x} \frac{\partial T}{\partial x} = T \frac{\partial^{2} X}{\partial x^{2}}$

(Eq. 8.69)

where we have applied the product rule (see Eq. 3.13). Please note that as X is only a function of x we know that ∂X∂t=0 $\frac{\partial X}{\partial t} = 0$ . Likewise, we know ∂T∂x=0 $\frac{\partial T}{\partial x} = 0$ as T is only a function of t.

8.3.3.2 Separation of Variables

We can now insert Eq. 8.68 and Eq. 8.69 into Eq. 8.69 and obtain

X∂T∂t∂T∂tc2T==c2T∂2X∂x2∂2X∂x2X $\begin{array}{l} X \frac{\partial T}{\partial t} & = & c^{2} T \frac{\partial^{2} X}{\partial x^{2}} \\ \frac{\frac{\partial T}{\partial t}}{c^{2} T} & = & \frac{\frac{\partial^{2} X}{\partial x^{2}}}{X} \end{array}$

si180_e (Eq. 8.70)

which we have rewritten such that all terms involving X are collected on the right-hand side whereas all terms involving T have been collected on the left-hand side of Eq. 8.70. This is the stated "separation" part.

8.3.3.3 Separation Into Two Differential Equations

After separation of variables we see that the left-hand side of Eq. 8.70 only depends on the independent variable t, whereas the right-hand side of the equation only depends on the independent variable x. Still for all values of t and x these two expressions must be equal. If this must be true for all potential values that the two independent variables could have, then there is only one possible solution: the expressions must both be a constant. We denote this constant − λ, which is, by convention, a negative term. This now turns our differential equation (8.66) into two differential equations:

∂T∂tc2T∂2X∂x2X==−λ−λ $\begin{array}{l} \frac{\frac{\partial T}{\partial t}}{c^{2} T} & = & - λ \\ \frac{\frac{\partial^{2} X}{\partial x^{2}}}{X} & = & - λ \end{array}$

si181_e

As in both cases, the respective functions are only functions of one independent variable. We can rewrite these equations as

dTdtc2T=−λ $\frac{\frac{d T}{d t}}{c^{2} T} = - λ$

si182_e (Eq. 8.71)

d2Xdx2X=−λ $\frac{\frac{d^{2 X}}{d x^{2}}}{X} = - λ$

si183_e (Eq. 8.72)

8.3.3.4 Solving T (t)

Eq. 8.71 is very easy to solve. We rewrite it as

dTdtdTdt+λc2T==−λc2T0 $\begin{array}{l} \frac{d T}{d t} & = & - λ c^{2} T \\ \frac{d T}{d t} + λ c^{2} T & = & 0 \end{array}$

si184_e (Eq. 8.73)

which yields a linear differential equation with constant coefficients (section 8.2.3). We find the characteristic polynomial (using a instead of λ in order to avoid confusion) to be

a−λc2a==0−λc2 $\begin{array}{l} a - λ c^{2} & = & 0 \\ a & = & - λ c^{2} \end{array}$

si185_e

Therefore we find

T(t)=e−λc2t $T (t) = e^{- λ c^{2} t}$

(Eq. 8.74)

8.3.3.5 Solving X (x)

Now that we have solved T (t) we shall proceed with solving X (x), which can be solved in exactly the same manner. From Eq. 8.72 we find

d2Xdx2=−λX $\frac{d^{2} X}{d x^{2}} = - λ X$

(Eq. 8.75)

Eq. 8.75 is a Helmholtz equation (see Eq. 8.26), which is a linear differential equation with constant coefficients. We find the characteristic polynomial (again using a instead of λ) to be

a2=−λ $a^{2} = - λ$

(Eq. 8.76)

There are three cases we need to distinguish here.

Case 1: λ = 0

For λ = 0 we find a₁ = a₂ = 0 to be a recurring root (see Eq. 8.34) in which case we find the solution to be

Xλ=0(x)=c1x+c2 $X_{λ = 0} (x) = c_{1} x + c_{2}$

Case 2: λ > 0

For λ > 0 Eq. 8.76 has two complex solutions: a1=iλ−−√ $a_{1} = i \sqrt{λ}$ and a2=−iλ−−√ $a_{2} = - i \sqrt{λ}$ . In this case we find the solution according to Eq. 8.36 to be

X1,λ>0(x)X2,λ>0(x)==c1cos(λ−−√x)c2sin(λ−−√x) $\begin{array}{l} X_{1, λ > 0} (x) & = & c_{1} cos (\sqrt{λ} x) \\ X_{2, λ > 0} (x) & = & c_{2} sin (\sqrt{λ} x) \end{array}$

si192_e

Case 3: λ < 0

For λ > < Eq. 8.76 has two noncomplex solutions, which are given by a1=λ−−√ $a_{1} = \sqrt{λ}$ and a2=−λ−−√ $a_{2} = - \sqrt{λ}$ in which case we find

Xλ<0(x)=c1eλ√x+c2e−λ√x $X_{λ < 0} (x) = c_{1} e^{\sqrt{λ} x} + c_{2} e^{- \sqrt{λ} x}$

Here, the two solutions are

X1,λ<0(x)X2,λ<0(x)==c1eλ√xc2sin(−λ−−√x) $\begin{array}{l} X_{1, λ < 0} (x) & = & c_{1} e^{\sqrt{λ} x} \\ X_{2, λ < 0} (x) & = & c_{2} sin (- \sqrt{λ} x) \end{array}$

si196_e

8.3.3.6 Preliminary Solution

We now have derived the solution that we can reassemble according to Eq. 8.67. The solution is given by

f(x,t)=⎧⎩⎨⎪⎪⎪⎪⎪⎪λ=0λ>0λ<0c1x+c2e−λc2t(c1cos(λ−−√x)+c2sin(λ−−√x))e−λc2t(c1eλ√x+c2e−λ√x) $f (x, t) = \{\begin{cases} λ = 0 & c_{1} x + c_{2} \\ λ > 0 & e^{- λ c^{2} t} (c_{1} cos (\sqrt{λ} x) + c_{2} sin (\sqrt{λ} x)) \\ λ < 0 & e^{- λ c^{2} t} (c_{1} e^{\sqrt{λ} x} + c_{2} e^{- \sqrt{λ} x}) \end{cases}$

si197_e (Eq. 8.77)

This is the solution to one-dimensional heat conduction equation. In all cases, the constants will be defined by the boundary conditions.

8.3.3.7 Principle of Superposition

You may wonder why the solution presented for λ > 0 and λ < 0 is the sum of the given solution and why (in case of λ > 0) we simply ignored the complex constant i. The reason is that a general solution to a differential equation must always include all potential solutions. Each of the solutions is multiplied by a constant. So in case a given solution does not apply to a given problem, the boundary values will determine that respective constant to be 0 and thus this particular solution disappears. However, whenever you have more than one solution you must always include all of them. This is referred to as the principle of superposition.

8.3.3.8 Initial Boundary Value Problem

In the next step we will look at an example where we apply boundary conditions and initial values to the problem. In this case, we will assume a rod that is fixed between two fixtures (see Fig. 8.5). At the time t = 0 the rod is assumed to have a given temperature distribution T (t = 0, x). For the sake of simplicity we will assume this distribution to be the constant temperature T₀. Heat will be conducted from the rod into the fixtures that are held at the temperature T = 0, which gives rise to the two boundary values we need for solving the heat equation: T (t, 0) = T (t, L) = 0. Please note that this temperature can be any temperature as long as the temperature profile of the rod is offset accordingly. Thus the temperature is given as a difference value to the temperature of the fixture.

f08-05-9781455731411 — Fig. 8.5 Analytical solution of the one-dimensional heat equation with initial values and boundary values. a) At t = 0 a rod fixed between two points whose ends are held at a constant temperature (in this case T = 0) has a given heat distribution. In this case the heat is distributed equally and the rod has the temperature T₀ at all points of x. b) At t = t₁ the temperature distribution will have changed. Heat will be conducted along the rod toward the two ends, giving rise to a temperature profile. The ambiance is assumed to be perfectly insulated.

We already know the solutions to the heat equation that is given by Eq. 8.77. We now apply the boundary conditions for the three cases we determined.