# MATLAB tutorial 2.6: Wave Equations

Wave Equations

The wave equation is a typical example of more general class of partial differential equations called hyperbolic equations. They occur in classical physics, geology, acoustics, electromagnetics, and fluid dynamics. Wave equations usually describe wave propagations in different media. Historically, the problem of a vibrating string such as that of a musical instrument was first studied by the French mathematician, mechanician, physicist, philosopher, and music theorist Jean le Rond d'Alembert. Because he was the first who found a solution of one-dimensional wave equation in 1746, the latter is usually referred to as d'Alembert's equation. Many others contributed to study of the wave equation, among first of them we mention Leonhard Euler (who discovered the wave equation in three space dimensions), Daniel Bernoulli ( the Euler–Bernoulli beam equation), and Joseph-Louis Lagrange (classical and celestial mechanics).

The wave equation for real-valued function $$u(x_1, x_2, \ldots , x_n , t)$$ of n spatial variables and a time variable t is

$\frac{\partial^2 u}{\partial t^2} = c^2 \nabla^2 u , \qquad \mbox{or} \qquad \square u =0 ,$
where c is a positive constant (having dimensions of speed) and
$\nabla^2 u \equiv \Delta u = \frac{\partial^2 u}{\partial x_1^2} + \frac{\partial^2 u}{\partial x_2^2} + \cdots + \frac{\partial^2 u}{\partial x_n^2} \qquad\mbox{and} \qquad \square u \equiv \square_c u = \frac{\partial^2 u}{\partial t^2} - c^2 \Delta u .$
The differential operator □ is called the d'Alembertian and △ is the Laplacian.

Suppose we have a medium whose displacement may be described by a scalar function u(x,t), where $${\bf x} \in \mathbb{R}^n , \quad t\in\mathbb{R} .$$ Suppose that the system is conservative and it has the Lagrangian $${\cal L} = \mbox{K} - \Pi ,$$ where the kinetic energy K and potential energy Π of the medium are

$\mbox{K} \left( u_t \right) = \frac{1}{2} \int \rho\,u_t\,{\text d}{\bf x} , \qquad \Pi \left( u \right) = \frac{1}{2} \int k \left\vert \nabla u \right\vert^2 {\text d}{\bf x} .$
Here ρ(x) is a mass-density and k(x) is a stiffness, both assumed positive, and $$u_t = \partial u/\partial t .$$ The corresponding action becomes
$S \left( u \right) = \int {\text d}t \, {\cal L} \left( u, u_t \right) = \int {\text d}t \int {\text d}{\bf x} \,\frac{1}{2} \left\{ \rho\,u_t^2 - k\left\vert \nabla u \right\vert^2 \right\} .$
Note that the kinetic and potential energies and the Lagrangian are functions of the spacial field and velocity at each fixed time, whereas the action is a functional of the space-time field u(x, t), obtained by integrating the Lagrangian with respect to time.

The Euler--Lagrange equation is satisfied by a stationary point (which is a function u(x, t)) of this action becomes

$\rho\, u_{tt} - \nabla \cdot \left( k\,\nabla u \right) =0 .$
If ρ and k are constants, then we get the wave equation
$\square_c u \equiv u_{tt} - c^2 \Delta u =0 \qquad\mbox{or}\qquad \frac{\partial^2 u}{\partial t^2} - c^2 \nabla^2 u =0 .$
The action functional for the wave equation is not positive definite. therefore, we cannot expect a solution of the wave equation to be a minimizer of the action, in general, only a critical point.   ■

We derive the wave equation in one space dimension that models the transverse vibrations of an elastic string. If such string is placed horizontally between end points x=0 and x=ℓ, it can freely vibrate within a vertical plane. Generally speaking it is not true; however, if displacements u(x,t) are small, we can assume that spring motion occur only within a plane perpendicular to its equilibrium horizontal position.

Perhaps the easiest case is observed with the investigation of mechanical vibrations. Suppose that an elastic string of length ℓ is tightly stretched between two supports at the same horizontal level, which we identify with x-axis. Then its end points may be taken as x=0 and x=ℓ. The elastic string may be thought of as a guitar or violin string, a guy wire, or possibly an electric power line. The positions of points on the string can be described by the displacement, which we denote by u(x,t), from the equilibrium horizontal position. If damping effects, such as air resistance, are neglected, and if the magnitude of the motion is not too large, then the displacement function satisfies the partial differential equation (called one dimensional wave equation)

$\square_c u =0 \qquad\mbox{or} \qquad u_{tt} = c^2 u_{xx}$
in the domain 0 < x < ℓ    0 < t < ∞. The constant coefficient c² is given by
$c^2 = T/\rho ,$
where T is the tension (force) in the string, and ρ is the mass per unit length of the string material (density). To describe the motion of the string completely, we need to impose some auxiliary conditions. Of these, we need to specify the initial displacement and its initial velocity
$u(x,0) = d(x), \qquad \dot{u} (x,0) \equiv \left. \frac{\partial u}{\partial t} \right\vert_{t=0} = v(x) ,$
where d and v are known functions. If we consider a ideal (and not realistic) case that the string has an infinite length, we arrive at so called the initial value problem:
$u_{tt} = c^2 u_{xx}, \qquad u(x,0) = d(x), \qquad \dot{u} (x,0) = v(x) , \qquad -\infty < x < \infty .$
This problem was first solved by Jean d'Alembert, which is presented below.

The one dimensional d'Alembertian operator can be recomposed into the product of the first order differential operators:

$\square_c \equiv \frac{\partial^2 }{\partial t^2} - c^2 \frac{\partial^2 }{\partial x^2} = \left( \frac{\partial}{\partial t} - c\,\frac{\partial}{\partial x} \right) \left( \frac{\partial}{\partial t} + c\,\frac{\partial}{\partial x} \right)$
This allows one to reduce the second order partial differential equation utt - c²uxx = 0 to two first order equations
$\frac{\partial v}{\partial t} - c\,\frac{\partial v}{\partial x} =0 \qquad\mbox{and}\qquad \frac{\partial v}{\partial t} + c\,\frac{\partial v}{\partial x} =0.$
By introducing a new variable ξ = x ±ct, each of the above equations is reduced to a simple ordinary differential equation
$\frac{{\text d} v}{{\text d} \xi} = \frac{\partial v}{\partial x}\,\frac{\partial x}{\partial \xi} + \frac{\partial v}{\partial t}\,\frac{\partial t}{\partial \xi} = \frac{\partial v}{\partial x} \pm \frac{1}{c}\,\frac{\partial v}{\partial t} =0,$
which can be integrated directly. This leads to the conclusion that a solution of the wave equation utt - c²uxx = 0 is the sum
$u(x,t) = f(x+ct) + g(x-ct)$
of two functions f(ξ) and g(ξ) of one variable. This formula represents a superposition of two waves, one traveling to the right and another traveling to the left, each with velocity x. However, in practice, traveling waves are excited by the initial disturbance
$u(x,0) = d(x) \qquad \mbox{and}\qquad \left. \frac{\partial u}{\partial t} \right\vert_{t=0} = v(x) ,$
where d(x) is the initial displacement (initial configuration) and v(x) is the initial velocity of the string. Upon substituting the general solution into the initial condition, we get two equations
$\begin{split} d(x) &= f(x) + g(x) , \\ v(x) &= c\, f' (x) - c\, g' (x) . \end{split}$
Integrating the latter, we obtain
$\int_0^x v(\xi )\,{text d}\xi = c\,f(x) - c\,g(x) .$
This enables us to express the general solution in terms of the initial displacement and the initial velocity (called the d'Alembert's formula)
$u(x,t) = \frac{d(x+ct) + d(x-ct)}{2} + \frac{1}{2c} \,\int_{x-ct}^{x+ct} v(\xi )\,{\text d}\xi .$

Let us consider the wave equation in semi-infinite domain 0 < x < ∞ For simplicity, we first assume that the boundary condition at left end x = 0 are of first type (Dirichlet):

$\left\{ \begin{array}{ll} \ddot{u} = c^2 u_{xx} , & 0 < x \infty ,\quad 0<t<t^{\ast} <\infty ; \\ u(0,t) =0, & 0 < t< t^{\ast} < \infty ; \\ u(x,0) = d(x) , \quad \dot{u} (x,0) = v(x) ,\quad & 0<x<\infty . \end{array} \right.$
To understand the solution, we assume temporally that input is a snakey. (A "snakey" is a slinky-like device that consists of a large concentration of small-diameter metal coils.) If an upward displaced pulse is introduced at the left end of the snakey, it will travel rightward across the snakey until it reaches the fixed end on the right side of the snakey. Upon reaching the fixed end, the single pulse will reflect and undergo inversion. That is, the upward displaced pulse will become a downward displaced pulse. Now suppose that a second upward displaced pulse is introduced into the snakey at the precise moment that the first crest undergoes its fixed end reflection. If this is done with perfect timing, a rightward moving, upward displaced pulse will meet up with a leftward moving, downward displaced pulse in the exact middle of the snakey. As the two pulses pass through each other, they will undergo destructive interference. The animation below shows several snapshots of the meeting of the two pulses at various stages in their interference.

We are now in a position to solve the general initial boundary value problem for the wave equation subject to homogeneous boundary conditions of the first type (Dirichlet's conditions are chosen for simplicity):

$\left\{ \begin{array}{ll} \ddot{u} = c^2 u_{xx} , & 0 < x < \ell , \quad 0<t<t^{\ast} <\infty ; \\ u(0,t) = u(\ell ,t) =0, & 0 < t< t^{\ast} < \infty ; \\ u(x,0) = d(x) , \quad \dot{u} (x,0) = v(x) ,\quad & 0<x<\ell . \end{array} \right.$
As usual, we consider an auxiliary problem that consists of the wave equation and homogeneous boundary conditions. We seek its partial nontrivial solutions in the form $$u(x,t) = X(x)\,T(t) .$$ Substituting this function into the wave equation gives
$\frac{X''}{X} = \frac{\ddot{T}(t)}{c^2 T(t)} = -\lambda .$
Separation of variables yields
$\begin{split} \ddot{T} + c^2 \lambda\, T(t) &=0 , \\ X'' (x) + \lambda \,X(x) &=0 . \end{split}$
Similarly to the previous sections, we substitute $$u(x,t) = X(x)\,T(t)$$ into the boundary conditions to obtain X(0) = 0 and X(ℓ) = 0. The differential equation together with these boundary conditions constitute the Sturm--Liouville problem:
$X'' (x) + \lambda \,X(x) =0, \qquad X(0) =0, \quad X(\ell ) =0 .$
Of the usual three possibilities for the parameter, λ = 0, λ = -α² < 0, and λ = α² > 0, only the last choice leads to nontrivial solutions. Corresponding to λ = α² > 0, with α > 0, the general solution to the differential equation $$X'' (x) + \alpha^2 \,X(x) =0$$ is
$X(x) = c_1 \cos \alpha x + c_2 \sin \alpha x .$
The boundary conditions X(0) = 0 and X(ℓ) = 0 indicate that c1 = 0 and $$c_2 \sin \alpha \ell =0 .$$ The last equation implies that $$\sin \alpha \ell =0$$ because otherwise c2 = 0 forces to obtain a trivial solution. Since we cannot effort vanishing c2, we have to choose α to be
$\alpha = \sqrt{\lambda} = \frac{n\pi}{\ell} , \qquad n=1,2,3,\ldots .$
Therefore, we obtain the eigenvalues and corresponding eigenfunctions:
$\lambda_n = \left( \frac{n\pi}{\ell} \right)^2 , \qquad X_n (x) = \sin \frac{n\pi x}{\ell}, \qquad n=1,2,3,\ldots .$
The general solution of the second order differential equation $$\ddot{T} + c^2 \lambda_n T(t) =0$$ becomes
$T_n (t) = A_n \cos \frac{cn\pi t}{\ell} + B_n \sin \frac{cn\pi t}{\ell} , \qquad n=1,2,3,\ldots .$
Since the eigenfunction $$X_n (x) = \sin \frac{n\pi x}{\ell}$$ can be multiplied by any nonzero constant, we find partial nontrivial solutions of our auxiliary problem:
$u_n (x,t) = X_n (x) \, T_n (t) = \sin \frac{n\pi x}{\ell} \left( A_n \cos \frac{cn\pi t}{\ell} + B_n \sin \frac{cn\pi t}{\ell} \right) , \qquad n=1,2,3,\ldots .$
Finally, we represent the required solution of the given initial boundary value problem as the sum of all partial nontrivial solutions
$u(x,t) = \sum_{n\ge 1} u_n (x,t) = \sum_{n\ge 1} X_n (x) \, T_n (t) = \sum_{n\ge 1} \sin \frac{n\pi x}{\ell} \left( A_n \cos \frac{cn\pi t}{\ell} + B_n \sin \frac{cn\pi t}{\ell} \right) .$
Setting t = 0, we get from the initial condition u(x,0) = d(x) that
$u(x,0) = d(x) = \sum_{n\ge 1} u_n (x,0) = \sum_{n\ge 1} A_n X_n (x) = \sum_{n\ge 1} A_n \sin \frac{n\pi x}{\ell} .$
Since the last series is just Fourier sine series for d(x), we can write
$A_n = \frac{2}{\ell} \int_0^{\ell} d(x)\, \sin \frac{n\pi x}{\ell} \,{\text d} x , \qquad n=1,2,3,\ldots .$
To determine Bn, we differentiate the general solution with respect to t and the set t = 0:
$\begin{split} \frac{\partial u}{\partial t} &= \sum_{n\ge 1} \left( -A_n \frac{cn\pi}{\ell} \, \sin \frac{cn\pi t}{\ell} + B_n \frac{cn\pi}{\ell}\, \cos \frac{cn\pi t}{\ell} \right) \sin \frac{n\pi x}{\ell} , \\ \left. \frac{\partial u}{\partial t} \right\vert_{t=0} &= v(x) = \sum_{n\ge 1} B_n \frac{cn\pi}{\ell}\, \sin \frac{n\pi x}{\ell} . \end{split}$
For this last sine Fourier series to be expansion of v(x), the total coefficient of sine must be given by
$B_n \frac{cn\pi}{\ell} = \frac{2}{\ell} \int_0^{\ell} v(x)\, \sin \frac{n\pi x}{\ell} \,{\text d} x ,$
from which we obtain
$B_n = \frac{2}{n\pi c} \int_0^{\ell} v(x)\, \sin \frac{n\pi x}{\ell} ,{\text d} x , \qquad n=1,2,3,\ldots .$
Recall that the constant c appearing in the solution of the initial boundary value problem is given by $$c= \sqrt{T/\rho} ,$$ where ρ is mass per unit length and T is the magnitude of the tension in the spring. When T is large enough, the vibrating string produces a musical sound. This sound, according to the solution
$u(x,t) = \sum_{n\ge 1} u_n (x,t) ,$
is the result of sum (superposition) of standing waves or normal modes.
Each normal mode can be written as
$u_n (x,t) = C_n \sin \left( \frac{n\pi c}{\ell} \,t + \phi_n \right) \sin \frac{n\pi x}{\ell} ,$
where
$C_n = \sqrt{A_n^2 + b_n^2} , \qquad \sin \phi_n = A_n /C_n \quad\mbox{or}\quad \cos \phi_n = B_n / C_n .$
For $$n= 1,2,3,\ldots$$ the standing waves are essentially the graphs of sine function $$\sin \frac{n\pi x}{\ell}$$ with a time-varying amplitude given by
$C_n \sin \left( \frac{n\pi c}{\ell} \,t + \phi_n \right) .$
Alternatively, we see from the formula for the normal mode un(x,t) that at a fixed value of x this function represents simple harmonic motion with amplitude $$C_n \sin \frac{n\pi x}{\ell}$$ and frequency fn = nc/(2ℓ). In other words, each point on a standing wave vibrates with a different amplitude but with the same frequency. The mode corresponding n = 1,
$u_1 (x,t) = C_1 \sin \left( \frac{c\pi}{\ell}\,t + \phi_1 \right) \sin \frac{\pi x}{\ell} ,$
is called the first standing wave, the first normal mode, or the fundamental mode of vibration, and its frequency is called the fundamental frequency:
$f_1 = \frac{c}{2\ell} = \frac{1}{2\ell} \sqrt{\frac{T}{\rho}} .$
The fundamental frequency or first harmonic is directly related to the pitch (or note) produced by a string instrument. It is apparent that the greater the tension on the string, the higher the pitch of the sound. The frequencies fn of the other normal modes, which are integer multiples of the fundamental frequency, are called overtones.

Standing wave patterns are wave patterns produced in a medium when two waves of identical frequencies interfere in such a manner to produce points along the medium that always appear to be standing still. These points that have the appearance of standing still are referred to as nodes. There are several frequencies with which the snakey can be vibrated to produce the patterns. Each frequency is associated with a different standing wave pattern. These frequencies and their associated wave patterns are referred to as harmonics. The two individual waves are drawn in blue and green and the resulting shape of the medium is drawn in black.

Solve the initial boundary value problem:

$\left\{ \begin{array}{ll} \ddot{u} = 9 \,u_{xx} , & 0 < x < 50, \quad 0<t<t^{\ast} <\infty ; \\ u(0,t) = u(50,t) =0, & 0 < t< t^{\ast} < \infty ; \\ u(x,0) = 3\,\sin (2\pi x) , \quad \dot{u} (x,0) = 5\,\sin (3\pi x) ,\quad & 0<x<50. \end{array} \right.$

Suppose we pluck a string by pulling it upward and release it from rest (see Fig.\;\ref{F555.1}). If the point of the pluck is in the third of a string of length $\ell$ (which is usually the case when playing guitar), we can model the vibration of the string by solving the following initial boundary value problem:

$\begin{array}{ll} \ddot{u} = c^2 \,u_{xx} , & 0< x< \ell , \quad 0< t <t^{\ast} <\infty ; \\ u(0,t) = u(\ell ,t) =0, & 0<t<t^{\ast} <\infty ; \\ u(x,0) = f(x) = \begin{cases} \dfrac{3x}{\ell} , & 0 < x < \frac{\ell}{3} , \\ \dfrac{3}{2\ell} \, (\ell -x) , & \frac{\ell}{3} < x < \ell , \end{cases} \\ \dot{u} (x,0) = v(x) \equiv 0 ,\quad & 0<x< \ell . \end{array}$

syms x n t
y1 = (x./50).*sin((n.*pi.*x)./120);
y2 = (120-x./70).*sin((n.*pi.*x)./120);
I1 = int(y1,x,0,50);
I2 = int(y2,x,50,120);
Cn = ((1/60).*I1)+((1/60).*I2);
utemp =
subs(Cn.*sin(n.*pi.*x.*(1/120)).*cos(n.*pi.*t.*(1/40)),n,1:10);
U = sum(utemp);
U = subs(U,t,1);
U = subs(U,x,0:120);
plot(0:120,U)

picture:

syms x n t
y1 = (x./50).*sin((n.*pi.*x)./120);
y2 = (120-x./70).*sin((n.*pi.*x)./120);
I1 = int(y1,x,0,50);
I2 = int(y2,x,50,120);
Cn = ((1/60).*I1)+((1/60).*I2);
utemp =
subs(Cn.*sin(n.*pi.*x.*(1/120)).*cos(n.*pi.*t.*(1/40)),n,1:10);
U = sum(utemp);
U = subs(U,t,5);
U = subs(U,x,0:120);
plot(0:120,U)

picture:

syms x n t
y1 = (x./50).*sin((n.*pi.*x)./120);
y2 = (120-x./70).*sin((n.*pi.*x)./120);
I1 = int(y1,x,0,50);
I2 = int(y2,x,50,120);
Cn = ((1/60).*I1)+((1/60).*I2);
utemp =
subs(Cn.*sin(n.*pi.*x.*(1/120)).*cos(n.*pi.*t.*(1/40)),n,1:10);
U = sum(utemp);
U = subs(U,t,15);
U = subs(U,x,0:120);
plot(0:120,U)

picture:

syms x n t
y1 = (x./50).*sin((n.*pi.*x)./120);
y2 = (120-x./70).*sin((n.*pi.*x)./120);
I1 = int(y1,x,0,50);
I2 = int(y2,x,50,120);
Cn = ((1/60).*I1)+((1/60).*I2);
utemp =
subs(Cn.*sin(n.*pi.*x.*(1/120)).*cos(n.*pi.*t.*(1/40)),n,1:10);
U = sum(utemp);
U = subs(U,t,25);
U = subs(U,x,0:120);
plot(0:120,U)

picture:

syms x n t
y1 = (x./50).*sin((n.*pi.*x)./120);
y2 = (120-x./70).*sin((n.*pi.*x)./120);
I1 = int(y1,x,0,50);
I2 = int(y2,x,50,120);
Cn = ((1/60).*I1)+((1/60).*I2);
utemp =
subs(Cn.*sin(n.*pi.*x.*(1/120)).*cos(n.*pi.*t.*(1/40)),n,1:10);
U = sum(utemp);
U = subs(U,t,0:35);
U = subs(U,x,15);
plot(0:35,U)

picture:

syms x n t
y1 = (x./50).*sin((n.*pi.*x)./120);
y2 = (120-x./70).*sin((n.*pi.*x)./120);
I1 = int(y1,x,0,50);
I2 = int(y2,x,50,120);
Cn = ((1/60).*I1)+((1/60).*I2);
utemp =
subs(Cn.*sin(n.*pi.*x.*(1/120)).*cos(n.*pi.*t.*(1/40)),n,1:10);
U = sum(utemp);
U = subs(U,t,0:35);
U = subs(U,x,25);
plot(0:35,U)

picture:

syms x n t
y1 = (x./50).*sin((n.*pi.*x)./120);
y2 = (120-x./70).*sin((n.*pi.*x)./120);
I1 = int(y1,x,0,50);
I2 = int(y2,x,50,120);
Cn = ((1/60).*I1)+((1/60).*I2);
utemp =
subs(Cn.*sin(n.*pi.*x.*(1/120)).*cos(n.*pi.*t.*(1/40)),n,1:10);
U = sum(utemp);
U = subs(U,t,0:35);
U = subs(U,x,47);
plot(0:35,U)

picture:

syms x n t
y1 = sin((n.*pi.*x)./90);
I1 = int(y1,x,10,15);
Kn = ((2.*n.*pi).*I1);
utemp = subs(Kn.*sin(n.*pi.*x.*(1/90)).*sin(n.*pi.*t.*(1/10)),n,1:10);
U = sum(utemp);
U = subs(U,t,1);
U = subs(U,x,0:90);
figure(1)
plot(0:90,U)
syms x n t
y1 = sin((n.*pi.*x)./90);
I1 = int(y1,x,10,15);
Kn = ((2.*n.*pi).*I1);
utemp = subs(Kn.*sin(n.*pi.*x.*(1/90)).*sin(n.*pi.*t.*(1/10)),n,1:10);
U = sum(utemp);
U = subs(U,t,3.5);
U = subs(U,x,0:90);
figure(1)
plot(0:90,U)

picture:

syms x n t
y1 = sin((n.*pi.*x)./90);
I1 = int(y1,x,10,15);
Kn = ((2.*n.*pi).*I1);
utemp = subs(Kn.*sin(n.*pi.*x.*(1/90)).*sin(n.*pi.*t.*(1/10)),n,1:10);
U = sum(utemp);
U = subs(U,t,11);
U = subs(U,x,0:90);
figure(1)
plot(0:90,U)

picture:

syms x n t
y1 = sin((n.*pi.*x)./90);
I1 = int(y1,x,10,15);
Kn = ((2.*n.*pi).*I1);
utemp = subs(Kn.*sin(n.*pi.*x.*(1/90)).*sin(n.*pi.*t.*(1/10)),n,1:10);
U = sum(utemp);
U = subs(U,t,18);
U = subs(U,x,0:90);
figure(1)
plot(0:90,U)

picture: ■

Sounds from a piano, unlike the guitar, are put into effect by striking strings. When a player presses a key, it causes a hammer to strike the strings. The corresponding IBVP is

\begin{equation*} \begin{split} &\text{(PDE)} \quad u_{tt} = c^2 u_{xx} , \quad 0<x<\ell, \ 0< t\le t^{\ast} < \infty , \\ &\text{(BC)} \hspace{6.2mm} u(0,t) =0 , \quad u (\ell ,t) =0 , \quad 0< t\le t^{\ast} < \infty , \\ &\text{(IC)} \qquad u(x,0) =0 , \quad u_t (x,0) = v(x) , \end{split} \end{equation*}
where the initial velocity is the step function
$v(x) = \begin{cases} 1 , & \mbox{ when } s< x< s+h , \\ 0, & \mbox{ otherwise. } \end{cases}$
Here s is a position of the left hammer's end and h is the width of the hammer. It is assumed that both s and s+h are within the string length ℓ. ■
Note that $$a_0 = u(x,t) .$$