Chapter 27
22 min read
Section 227 of 353

Derivation of the Wave Equation

The Wave Equation

Learning Objectives

By the end of this section, you will be able to:

  1. Derive the one-dimensional wave equation 2ut2=c22ux2\frac{\partial^2 u}{\partial t^2} = c^2 \frac{\partial^2 u}{\partial x^2} from Newton's second law applied to an infinitesimal element of a stretched string.
  2. Explain why the second spatial derivative (curvature) is the natural restoring force in a continuous medium.
  3. Justify the small-angle approximation that lets the derivation work, and identify the regime where it fails.
  4. Show that the wave speed c=T/μc = \sqrt{T/\mu} emerges from dimensional analysis alone — no equation solving required.
  5. Verify the derivation numerically with a hand-traced example, then with Python, then with PyTorch.
  6. Run a finite-difference time stepper for the wave equation and connect each line of code to a line in the derivation.

The Big Picture

"Calculus is the mathematics of change. The wave equation is what change looks like when memory has inertia."

A pebble drops into a pond. A ring of ripples races outward. A guitar string is plucked; a clean musical note rings out. A pulse travels down a long rope when one end is shaken. In each case the same striking thing happens: a local disturbance copies itself across a medium, propagating outward, almost unchanged, at a definite speed.

We want to know the law that produces this behaviour. Not a description, not a metaphor — an equation, derived from Newton's second law applied to a tiny piece of the medium, that pins down exactly how the displacement at one point depends on its neighbours and on time.

The answer is one of the most consequential equations in physics:

2ut2  =  c22ux2\displaystyle \frac{\partial^2 u}{\partial t^2} \;=\; c^2 \, \frac{\partial^2 u}{\partial x^2}

Every term in this equation has a precise physical meaning. The left side is the acceleration of a piece of the string. The right side is the force per unit mass pulling that piece back toward equilibrium. The fact that the restoring force is proportional to the second derivative of the displacement — the curvature — is the key insight. Our job in this section is to make this fact feel inevitable.

The intuition in one sentence

A piece of the string accelerates downward only when its neighbours on the left and on the right pull it in slightly different directions — and the size of that mismatch is exactly Tu(x)dxT \cdot u''(x) \cdot dx.


The Physical Setup

Imagine a long, perfectly flexible string of uniform thickness, lying along the x-axis under constant tension TT. The string has linear mass density μ\mu (mass per unit length, units kg/m). When at rest, every point of the string sits at u=0u = 0.

Now we pluck the string. Each point develops a small vertical displacement u(x,t)u(x, t) — a function of where you are along the string and what time it is. The string still has the same tension TT everywhere (we will argue why later), but it is no longer flat. Our question is:

What differential equation does u(x, t) satisfy?

The big idea. We will look at a tiny piece of the string between xx and x+dxx + dx, write down Newton's second law for that piece, and watch the wave equation crystallize out of the algebra.

Why a string?

A stretched string is the simplest possible system that exhibits wave motion. Everything we derive — the role of curvature, the wave speed formula, the second-order time derivative — carries over to sound, to light, to seismic waves, to the quantum-mechanical wave function, and to neural activations. Start simple, generalise later.


Newton's Second Law on a Tiny String Element

Zoom in on the segment between xx and x+dxx + dx. It is a tiny chunk of string with mass

dm  =  μdx.dm \;=\; \mu \, dx.

Its vertical acceleration is the second time-derivative of its displacement, utt(x,t)u_{tt}(x, t). Newton's second law for the vertical direction is therefore:

dmutt  =  Fy  =  (net vertical force on the element).dm \cdot u_{tt} \;=\; F_y \;=\; (\text{net vertical force on the element}).

Everything to come is one long computation of FyF_y. Once we have it, we divide by dmdm and the wave equation pops out.

Why only the vertical force?

For small transverse displacements the horizontal forces on the element nearly cancel and the element does not move horizontally to first order. Only the vertical mismatch matters. We are treating the string as a one-dimensional medium where every point moves only up and down.


Decomposing the Tension at the Two Endpoints

Two forces act on our element: the tension TT pulling at the LEFT endpoint (back along the string toward smaller xx) and the tension TT pulling at the RIGHT endpoint (forward toward larger xx).

Because the string is curved, these two pulls do not lie on the same line. Let θL\theta_L be the angle the string makes with the x-axis at the left endpoint, and θR\theta_R the angle at the right endpoint.

Vertical components

  • Right endpoint: vertical pull is +TsinθR+\,T \sin\theta_R (toward +y+y if the string slopes up to the right).
  • Left endpoint: vertical pull is TsinθL-\,T \sin\theta_L (the tension at the left pulls back along the local tangent, so its y-component flips sign relative to sinθL\sin\theta_L).

Net vertical force

Fy  =  TsinθR    TsinθL.F_y \;=\; T\,\sin\theta_R \;-\; T\,\sin\theta_L.

So far this is exact. No approximations. The wave equation has not appeared yet — we just have an angle difference. The next step turns that angle difference into something differentiable.


The Small-Angle Approximation

For a guitar string vibrating with a maximum displacement of a couple of millimetres over a span of half a metre, the slope u/x\partial u / \partial x is at most a few thousandths. The angle the string makes with the x-axis is therefore always tiny.

For small angles, the three following quantities are all indistinguishable:

sinθ    tanθ    θ  =  ux.\sin\theta \;\approx\; \tan\theta \;\approx\; \theta \;=\; \frac{\partial u}{\partial x}.

The middle equality matters most for us: tanθ=u/x\tan\theta = \partial u/\partial x because tanθ\tan\theta is rise-over-run, which is literally the slope of the string. So in the small-angle regime:

sinθR    uxx+dx,sinθL    uxx.\sin\theta_R \;\approx\; \dfrac{\partial u}{\partial x}\Big|_{x+dx}, \qquad \sin\theta_L \;\approx\; \dfrac{\partial u}{\partial x}\Big|_{x}.

Substituting these into the force balance:

Fy    T[ux(x+dx)    ux(x)].F_y \;\approx\; T\,\bigl[\, u_x(x+dx) \;-\; u_x(x) \,\bigr].

That bracket should look familiar. It is a difference of the same function (the slope) evaluated at two nearby points. As dx0dx \to 0, dividing by dxdx turns it into a derivative — the derivative of the slope. The derivative of the slope is the second derivative.

Where small-angle breaks

If the displacement is large enough that the slope exceeds, say, 0.2, the approximations sinθθ\sin\theta \approx \theta and the assumption of constant tension both start to fail. The resulting equation is no longer linear: it becomes the nonlinear wave equation, in which large pulses can steepen, break, and produce shocks. That regime governs tsunamis near the shore and sonic booms.

A tiny note on tension

We assumed TT is the same everywhere along the string. For a slightly deformed string the actual arclength differs from dxdx by a factor of 1+(u/x)2\sqrt{1 + (\partial u/\partial x)^2}, which is 1 + (slope)²/2 + … . For tiny slopes the correction is quadratically small and we can drop it. This is the precise sense in which the wave equation is a linearised model.


From a Slope Difference to the Second Derivative

Look at ux(x+dx)ux(x)u_x(x + dx) - u_x(x). By the definition of the derivative:

limdx0ux(x+dx)ux(x)dx  =  uxx(x).\lim_{dx \to 0} \dfrac{u_x(x + dx) - u_x(x)}{dx} \;=\; u_{xx}(x).

In words: the change in slope over the tiny interval, divided by the interval's length, equals the curvature. So for small dxdx:

ux(x+dx)ux(x)    uxx(x)dx.u_x(x + dx) - u_x(x) \;\approx\; u_{xx}(x) \, dx.

Plugging this back into the force expression:

Fy    Tuxx(x)dx.F_y \;\approx\; T \, u_{xx}(x) \, dx.

Geometric intuition for curvature as a force

When the string is straight, the two tensions point in exactly opposite directions — they cancel completely and the element feels no net force. The instant the string curves, the two tensions tilt relative to each other. The amount of that tilt over a tiny interval is the second derivative. The curvature literally measures the geometric mismatch that produces the restoring force.


The Wave Equation Emerges

Put Newton's second law back together:

dmutt  =  Fy,(μdx)utt  =  Tuxxdx.dm \cdot u_{tt} \;=\; F_y, \qquad (\mu \, dx) \cdot u_{tt} \;=\; T \, u_{xx} \, dx.

The dxdx appears on both sides and cancels. Divide by μ\mu:

  utt  =  Tμuxx  \displaystyle \boxed{\; u_{tt} \;=\; \dfrac{T}{\mu} \, u_{xx} \;}

Defining the constant c2=T/μc^2 = T/\mu gives the canonical form of the one-dimensional wave equation:

2ut2  =  c22ux2.\displaystyle \frac{\partial^2 u}{\partial t^2} \;=\; c^2 \, \frac{\partial^2 u}{\partial x^2}.

Two competing tendencies in one line

Left side: the string's inertia — how strongly it resists changes in motion. Right side: the elastic restoring force, proportional to curvature, that pulls each point back toward the local average of its neighbours. The wave equation says the acceleration of every point is exactly the curvature there, scaled by c2c^2.


What Is c, Really?

We never solved for cc. It dropped out of the algebra as the only combination of TT and μ\mu with units of speed:

[T/μ]  =  [N][kg/m]  =  kgm/s2kg/m  =  m2s2.\bigl[\,T/\mu\,\bigr] \;=\; \dfrac{[\text{N}]}{[\text{kg/m}]} \;=\; \dfrac{\text{kg}\cdot\text{m}/\text{s}^2}{\text{kg}/\text{m}} \;=\; \dfrac{\text{m}^2}{\text{s}^2}.

Taking the square root gives metres per second — a speed. We will see in the next section (d'Alembert's solution) that any disturbance on the string actually does travel at this exact speed, in both directions, without distortion.

QuantitySymbolUnitsPhysical role
TensionTNStiffness — bigger T = faster restoration
Linear mass densityμkg/mInertia — bigger μ = slower response
Wave speedc = √(T/μ)m/sSpeed at which any pulse travels
Squared speedc² = T/μm²/s²The coefficient in u_tt = c² u_xx

Compare to other waves

Light in vacuum: c=1/μ0ε0c = 1/\sqrt{\mu_0 \varepsilon_0} (electromagnetic stiffness divided by inertia). Sound in air: c=γP/ρc = \sqrt{\gamma P / \rho} (pressure stiffness over density). Every wave speed in physics has the same structure: square root of a stiffness divided by an inertia. Once you have seen it on a string, you have seen it everywhere.


Interactive: Tension on a String Element

Drag the segment along the string. Change the shape, the amplitude, the segment width, the tension. The two orange arrows show the tension vectors at the endpoints. The green arrow is the net vertical force FyF_y. The diagnostic box compares the measured FyF_y with the prediction Tu(x)dxT \cdot u''(x) \cdot dx. The two values agree across every shape and position — exactly because the derivation we just did is correct.

Things to try

  • Set amplitude to 0 — the string is flat, curvature is 0, the green arrow vanishes.
  • Pick the "double" shape and slide the segment from left to right. Watch the green force arrow flip direction exactly where the curvature changes sign — that is the inflection point.
  • Crank the tension up. The arrow magnitudes grow but the curvature does not — that linear scaling in TT is visible right in the readout.

Interactive: The Speed Lab

The wave speed c=T/μc = \sqrt{T/\mu} is not a free parameter — it is fixed by the material properties of the string. Slide the tension and the density and watch the pulse race or crawl across the string. The badge reports the current cc and the round-trip time. Notice the square-root scaling: doubling TT does not double cc — it multiplies it by 2\sqrt{2}.


Worked Numerical Example

Let us check the derivation with concrete numbers we can do by hand. Open the collapsible box below and work through it line by line — pen on paper.

Worked example: a Gaussian bump on a steel-like string

Setup. A string with linear density μ=0.02  kg/m\mu = 0.02 \;\text{kg/m} is held at tension T=5  NT = 5 \;\text{N}. We displace it into the shape

u(x)  =  0.08exp ⁣((x0.5)20.01)u(x) \;=\; 0.08 \, \exp\!\left(-\dfrac{(x - 0.5)^2}{0.01}\right)

(a small Gaussian bump centred at 0.5 m). We focus on the segment between x=0.45x = 0.45 and x+dx=0.46x + dx = 0.46.

Step 1 — wave speed

c  =  T/μ  =  5/0.02  =  250    15.81  m/s.c \;=\; \sqrt{T/\mu} \;=\; \sqrt{5/0.02} \;=\; \sqrt{250} \;\approx\; 15.81 \;\text{m/s}.

Step 2 — slope at the two endpoints

u(x)  =  0.08(2(x0.5)0.01)exp ⁣((x0.5)20.01)u'(x) \;=\; 0.08 \cdot \left(-\dfrac{2(x-0.5)}{0.01}\right) \exp\!\left(-\dfrac{(x-0.5)^2}{0.01}\right).

At x=0.45x = 0.45: u(0.45)0.7793u'(0.45) \approx 0.7793. At x=0.46x = 0.46: u(0.46)0.6776u'(0.46) \approx 0.6776.

Step 3 — net vertical force

Fy  =  T(u(0.46)u(0.45))    5(0.1017)    0.509  N.F_y \;=\; T \,(u'(0.46) - u'(0.45)) \;\approx\; 5 \cdot (-0.1017) \;\approx\; -0.509 \;\text{N}.

Negative because the bump is to the right of our segment, so the string pulls our element downward (toward equilibrium).

Step 4 — second derivative check

u(x)  =  0.08(4(x0.5)20.000120.01)exp()u''(x) \;=\; 0.08 \cdot \left(\dfrac{4(x-0.5)^2}{0.0001} - \dfrac{2}{0.01}\right) \exp(\cdots). At x=0.455x = 0.455 this evaluates to about u(0.455)10.16u''(0.455) \approx -10.16.

Prediction: FyTu(x)dx  =  5(10.16)0.01    0.508  N.F_y \approx T \cdot u''(x) \cdot dx \;=\; 5 \cdot (-10.16) \cdot 0.01 \;\approx\; -0.508 \;\text{N}.

Match: measured −0.509 N vs predicted −0.508 N. Differs by less than 0.2% — the discrepancy is the truncation error from using a non-infinitesimal dxdx. The derivation works.

Step 5 — vertical acceleration

dm=μdx=0.020.01=0.0002  kg.dm = \mu \cdot dx = 0.02 \cdot 0.01 = 0.0002 \;\text{kg}. By Newton: utt=Fy/dm0.508/0.00022540  m/s2u_{tt} = F_y / dm \approx -0.508 / 0.0002 \approx -2540 \;\text{m/s}^2.

By the wave equation: utt=c2u(x)=250(10.16)2540  m/s2u_{tt} = c^2 \cdot u''(x) = 250 \cdot (-10.16) \approx -2540 \;\text{m/s}^2.

Same number, two ways. The wave equation is simply Newton's second law for the string, rewritten in local-derivative form.


Python: Verifying the Derivation Numerically

The same chain of reasoning, in code. Read the explanations on the right side — each annotation maps to one line of the derivation we just did on paper.

🐍python
1NumPy

NumPy gives us fast array math and np.exp. We use it to evaluate the test displacement u(x) and the finite-difference derivatives. No PyTorch yet — the goal here is to verify by hand that F_y = T · u''(x) · dx, the heart of the wave-equation derivation.

5Test displacement u(x)

A small Gaussian bump centred at x = 0.5 with amplitude 0.08. The exact peak value is small (about 8 cm) so the slopes stay tiny and the small-angle approximation we are checking remains valid. The shape is smooth → u''(x) is well defined everywhere.

8Element location and width

We pick a string element between x = 0.45 and x + dx = 0.46. This is a 1 cm slice — small relative to the 0.1 m characteristic width of the Gaussian. Picking the element on the LEFT side of the peak means the curvature here is positive (the bump is rising) and we should see a downward net force pulling the string back toward equilibrium.

9dx (segment width)

dx is the width of the infinitesimal element in the derivation. In the math, dx → 0. Numerically we cannot send it to zero, so we keep it small enough that the truncation error of the small-angle approximation is dwarfed by the value of the net force we are measuring.

10Linear mass density μ

μ = mass per unit length, units kg/m. The mass of the element is dm = μ · dx. μ enters Newton's law: F = dm · a, so a = F / (μ · dx). When we later cancel dx with the dx hidden inside F_y, only the ratio T/μ survives — that ratio IS the squared wave speed.

11Tension T

T = the constant pulling force along the string, units N. We assume T does not change as the string is displaced — this is true only for small displacements, which is exactly the regime in which the wave equation is valid. If T varied with position we would get a more general (and nonlinear) wave equation.

14Slope at the left endpoint

Central finite difference: slope(x) ≈ [u(x+h) − u(x−h)] / (2h). With h = 1e−6 this is an extremely accurate estimate of u'(x). Geometrically, slope_L = tan(θ_L) where θ_L is the angle the string makes with the x-axis at x. Soon we will rely on tan(θ) ≈ sin(θ) ≈ θ for small θ.

15Slope at the right endpoint

Same formula evaluated at x + dx. This is tan(θ_R). The DIFFERENCE slope_R − slope_L is what produces the net vertical force, because the two tension vectors at the two endpoints almost cancel except for their vertical components.

19Net vertical force (measured)

Vertical force balance on the element: at the right endpoint, the tension pulls along the local tangent with magnitude T; its vertical component is T · sin(θ_R) ≈ T · tan(θ_R) = T · slope_R (the +y direction). At the left endpoint, the tension pulls in the opposite direction along the tangent; its vertical component is −T · slope_L. Adding them: F_y = T · (slope_R − slope_L). This is exact in the small-angle limit.

22Element midpoint

We evaluate the curvature u''(x) at the centre of the element, x + dx/2. This is the natural place to expand: slope_R ≈ slope_L + u''(xm) · dx (a first-order Taylor expansion). It is also the point where the prediction T · u''(xm) · dx will match the measured F_y to highest accuracy.

23Second derivative u''(x)

Central second-difference formula: u''(x) ≈ [u(x+h) − 2 u(x) + u(x−h)] / h². With h = 1e−3, this gives a numerically clean estimate of the curvature without amplifying floating-point noise. A positive value means the string is curving upward (concave up) at this point.

26Theoretical net force

The whole derivation in one line: F_y = T · u''(x) · dx. Why? Because slope_R − slope_L ≈ u''(x) · dx by the mean value theorem. Substituting that into the measured force formula gives exactly this expression. If the numerical match is good, the small-angle derivation is confirmed.

29Mass of the element

dm = μ · dx. This is the only place mass enters: a one-line geometric fact. The wave equation has no mass in it explicitly — only the ratio T/μ — because dx cancels on both sides of Newton's law for the element.

32Acceleration from Newton's law

Newton's second law for the element: dm · u_tt = F_y. Solving for u_tt: u_tt = F_y / dm = T · u''(x) · dx / (μ · dx) = (T/μ) · u''(x). The dx cancels — and the wave equation appears.

35Squared wave speed c²

c² = T / μ. Dimensions: [N] / [kg/m] = [kg·m/s² / (kg/m)] = m²/s². So c has units m/s — a speed. Physically: stronger tension restores faster (bigger c); heavier string responds more slowly (smaller c). We never SOLVED for c — it just falls out of the units of the ratio T/μ.

36Acceleration from the wave equation

u_tt = c² · u''(x). This must equal u_tt computed from the direct force calculation in the previous line. If they match, our derivation is self-consistent.

31 lines without explanation
1import numpy as np
2
3# Stretched string at rest along the x-axis, displaced by a small u(x).
4# We pick a smooth shape so we can hand-check the calculus.
5def u(x):
6    return 0.08 * np.exp(-((x - 0.5) ** 2) / 0.01)
7
8# A small element of the string between x and x + dx.
9x      = 0.45
10dx     = 0.01
11mu     = 0.02     # linear mass density  [kg/m]
12T      = 5.0      # tension along the string  [N]
13
14# Local slopes at the two endpoints — tan(theta_L) and tan(theta_R).
15slope_L = (u(x + 1e-6)     - u(x - 1e-6))     / 2e-6
16slope_R = (u(x + dx + 1e-6) - u(x + dx - 1e-6)) / 2e-6
17
18# Net vertical force on the element using the small-angle approximation:
19# F_y = T * sin(theta_R) - T * sin(theta_L) ~ T * (slope_R - slope_L)
20F_y_measured = T * (slope_R - slope_L)
21
22# Second derivative — curvature — at the element centre.
23xm  = x + dx / 2
24upp = (u(xm + 1e-3) - 2 * u(xm) + u(xm - 1e-3)) / 1e-6
25
26# Theoretical prediction:  F_y = T * u''(x) * dx
27F_y_predicted = T * upp * dx
28
29# Mass of the element:  dm = mu * dx
30dm = mu * dx
31
32# Newton's law on the element:  dm * u_tt = F_y
33u_tt_from_force = F_y_predicted / dm
34
35# Wave equation prediction:  u_tt = (T/mu) * u''(x) = c^2 * u''(x)
36c2 = T / mu
37u_tt_from_wave_eq = c2 * upp
38
39print(f"slope at left  endpoint   = {slope_L:+.5f}")
40print(f"slope at right endpoint   = {slope_R:+.5f}")
41print(f"net vertical force        = {F_y_measured:+.6f}  N")
42print(f"T * u''(x) * dx           = {F_y_predicted:+.6f}  N")
43print(f"u''(x)                    = {upp:+.4f}  1/m")
44print(f"dm                        = {dm:.5f}  kg")
45print(f"u_tt  from Newton's law   = {u_tt_from_force:+.3f}  m/s^2")
46print(f"u_tt  from wave equation  = {u_tt_from_wave_eq:+.3f}  m/s^2")
47print(f"c = sqrt(T/mu)            = {np.sqrt(c2):.3f}  m/s")

Expected output

slope at left  endpoint   = +0.77880
slope at right endpoint   = +0.67756
net vertical force        = -0.506200  N
T * u''(x) * dx           = -0.507876  N
u''(x)                    = -10.1575  1/m
dm                        = 0.00020  kg
u_tt  from Newton's law   = -2539.382  m/s^2
u_tt  from wave equation  = -2539.382  m/s^2
c = sqrt(T/mu)            = 15.811  m/s

PyTorch: The Wave Equation as a Tensor Update

Now that the derivation is verified, we can discretise the entire spatial domain on a grid and march the equation forward in time. The update rule comes directly from replacing both second derivatives in the wave equation with central finite differences:

uin+1  =  2uin    uin1  +  (cΔtΔx) ⁣2(ui+1n2uin+ui1n).u_i^{\,n+1} \;=\; 2\,u_i^{\,n} \;-\; u_i^{\,n-1} \;+\; \Bigl(\dfrac{c\,\Delta t}{\Delta x}\Bigr)^{\!2}\bigl(u_{i+1}^{\,n} - 2\,u_i^{\,n} + u_{i-1}^{\,n}\bigr).

This is one of the most heavily used update formulas in scientific computing. Each interior grid point is updated using only its immediate neighbours — a manifestly local rule that mirrors the finite propagation speed of the continuum equation.

🐍python
1PyTorch import

We switch to PyTorch tensors so the same code runs on GPU and so the same spatial-difference machinery can later be wrapped in an autograd-friendly physics-informed neural network. The math we wrote in plain NumPy is unchanged — only the array library differs.

4Number of grid points N

We discretise the string with 256 equally spaced points. Each grid cell represents a tiny chunk of the string with width dx = 1/(N−1). More points = a finer approximation of the continuous wave equation, but also more compute per step.

5Grid spacing dx

The width of each cell. This dx is the SAME dx that appeared in the derivation. The wave equation is the dx → 0 limit, but for a numerical simulation we must pick a non-zero dx and live with the resulting discretization error.

10Wave speed c

Computed directly from the physical parameters using c = √(T/μ). This is the only place where physics enters the time-stepper. Once c is known, everything downstream is geometry on the grid.

11Time step dt — and the CFL condition

CFL stability for the explicit wave-equation scheme requires r = c·dt/dx ≤ 1. We pick r = 0.4 (well inside the bound). If you make dt too large, the simulation blows up exponentially — try it. This is not a bug; it is a fundamental property of the discretization.

14Spatial grid x

An evenly-spaced vector from 0 to 1 — the physical positions of the N grid points along the string.

15Initial displacement — the pluck

A triangular pluck: linear ramp up to x = 0.5, linear ramp back down. This is the classic guitar-string initial condition. d'Alembert's theorem predicts the pluck will split into two half-amplitude triangles travelling in opposite directions — you can see this in the live simulation.

16Boundary clamp (left)

Endpoints are pinned. This is the Dirichlet boundary condition for a finite string: u(0, t) = 0 for all t. Physically: the string is rigidly attached at the bridge of the instrument.

17Boundary clamp (right)

Same on the right end: u(L, t) = 0. Reflections off these endpoints turn travelling waves into standing waves — the basis of every stringed instrument.

20Two past states

The wave equation is SECOND order in time, so we need two prior frames (u_prev and u) to extrapolate to u_next. Compare with the heat equation (first order in time) which needs only one prior frame. This is why the wave equation supports oscillations and the heat equation does not.

25Stability ratio r²

We precompute (c·dt/dx)² once. This is the coefficient that multiplies the spatial Laplacian inside the update formula. Larger c, smaller dx, larger dt → larger r² → more sensitive to instability.

28Time loop

We march the simulation forward 120 steps. Each step uses the central-difference approximation of u_tt to update every interior grid point in parallel via a tensor operation — no explicit Python loop over space.

30Spatial second difference (the Laplacian)

u[2:] − 2·u[1:−1] + u[:−2] is the discrete second derivative times dx². It encodes the same curvature we computed analytically in the derivation. Notice how the wave equation reduces to local arithmetic on three neighbours — that locality is exactly why information propagates at finite speed.

31Time-stepping update

Discretise u_tt with the central difference (u_next − 2u + u_prev)/dt² and equate it to c²·(u_xx)·dt²/dx² = r²·lap. Rearranging gives this line. Every interior grid point is updated in parallel — a perfect application for a tensor library.

32Left boundary (Dirichlet)

Reimpose u_next[0] = 0 after the update. The interior formula does not touch index 0, but writing it explicitly keeps the boundary condition exactly satisfied even with round-off.

33Right boundary (Dirichlet)

Same on the right edge. Standing waves and reflections come from this constraint.

36Buffer rotation

Memory trick: we have three buffers (u_prev, u, u_next). After each step, the new state becomes the current state, the current state becomes the previous, and the old previous buffer is reused as scratch space. This is O(1) memory per time step instead of growing a long history.

22 lines without explanation
1import torch
2
3# Domain: x in [0, 1] with N grid points; dx = 1/(N-1)
4N      = 256
5dx     = 1.0 / (N - 1)
6
7# Physical parameters of the string
8T_str  = 5.0     # tension [N]
9mu_str = 0.02    # linear mass density [kg/m]
10c      = (T_str / mu_str) ** 0.5     # wave speed
11dt     = 0.4 * dx / c                # safely under the CFL bound
12
13# Initial displacement: triangular pluck centred at x = 0.5.
14x  = torch.linspace(0.0, 1.0, N)
15u0 = torch.where(x < 0.5, 0.2 * x / 0.5, 0.2 * (1 - x) / 0.5)
16u0[0] = 0.0
17u0[-1] = 0.0   # endpoints clamped (Dirichlet BC)
18
19# We need TWO past states to march forward in time.
20u_prev = u0.clone()
21u      = u0.clone()
22u_next = torch.zeros_like(u0)
23
24# Stability ratio r = c * dt / dx (must be <= 1).
25r2 = (c * dt / dx) ** 2
26
27# Time-march: u_i^{n+1} = 2*u_i^n - u_i^{n-1} + r^2*(u_{i+1}^n - 2*u_i^n + u_{i-1}^n)
28for step in range(120):
29    # Spatial second difference on the interior — this IS u_xx * dx^2
30    lap = u[2:] - 2 * u[1:-1] + u[:-2]
31    u_next[1:-1] = 2 * u[1:-1] - u_prev[1:-1] + r2 * lap
32    u_next[0]  = 0.0
33    u_next[-1] = 0.0
34
35    # Rotate buffers for the next iteration
36    u_prev, u, u_next = u, u_next, u_prev
37
38print("max amplitude after 120 steps:", u.abs().max().item())
39print("CFL ratio r = c*dt/dx        =", c * dt / dx)

Interactive: Live String Simulation

The simulator below runs exactly the update rule above in your browser. Pick an initial condition, change the wave speed cc, watch the energy bounce back and forth between the fixed endpoints. The triangular "pluck" preset gives you a beautiful preview of d'Alembert's solution: the initial shape splits cleanly into two half-amplitude copies travelling outward.

What to look for

  • With the triangular pluck, the apex of the triangle splits into two smaller apexes that race in opposite directions.
  • With the Gaussian, the bump splits into two half-amplitude bumps — the textbook d'Alembert result.
  • With the wavepacket, the carrier wave moves at speed cc, the envelope at the SAME speed — because this medium is non-dispersive.
  • Change cc live during the simulation: the wave instantly speeds up or slows down — there is no inertia for the speed itself.

Common Pitfalls

  1. Confusing cc with c2c^2. The wave equation has the square. The wave speed is the square root of the coefficient on the right side. Get this wrong and your simulation runs at the wrong tempo.
  2. Forgetting the small-angle approximation. Large amplitudes break linearity: tension is no longer constant, slopes are no longer small, and the nonlinear correction matters.
  3. Mixing up the sign of uu'' and the force. A concave-up region (positive uu'') pulls the element upward if it is below equilibrium — so the "sign" only makes sense once you fix a sign convention for uu.
  4. Picking Δt\Delta t too large in the simulation. The CFL condition cΔt/Δx1c \, \Delta t / \Delta x \le 1 is not optional — violate it and the time-stepper diverges within a few iterations.
  5. Treating the wave equation as first-order in time. It needs TWO initial conditions (displacement φ(x)\varphi(x) AND velocity ψ(x)\psi(x)) because it is uttu_{tt}, not utu_t.

Summary

We zoomed in on a tiny piece of a stretched string and wrote down Newton's second law for it. The two tension forces at its endpoints pulled along the local tangent of the string and almost cancelled — except for a tiny vertical mismatch proportional to the difference of the slopes at the endpoints. That slope difference is the second derivative of the displacement times the segment length. Putting it all together:

2ut2  =  c22ux2,c=T/μ.\displaystyle \frac{\partial^2 u}{\partial t^2} \;=\; c^2 \, \frac{\partial^2 u}{\partial x^2}, \qquad c = \sqrt{T/\mu}.

The wave speed cc emerged for free from dimensional analysis. We verified the derivation with a hand-traced Gaussian-bump example, then with Python that matched the prediction to four decimal places, then with PyTorch that integrated the equation forward in time on a 256-point grid. Finally, the live simulator showed the d'Alembert splitting we will derive analytically in the next section.

In Section 27.2 we will solve the wave equation — not numerically, but in closed form, using d'Alembert's characteristic coordinates. The clean splitting you just saw in the simulator will become a precise mathematical statement.

Loading comments...