Software 1.7.1 (Symbolic solution with pdsolve). SymPy solves linear first-order equations with constant coefficients by the method of characteristics, and can verify a proposed general solution:

>>> import sympy as sp
>>> x, y = sp.symbols('x y'); u = sp.Function('u')
>>> sp.pdsolve(2 * u(x, y).diff(x) + 3 * u(x, y).diff(y), u(x, y))
Eq(u(x, y), F(3*x - 2*y))
>>> sp.pdsolve(u(x, y).diff(x) + u(x, y).diff(y) + u(x, y) - x, u(x, y))
Eq(u(x, y), x + F(x - y)*exp(-x/2 - y/2) - 1)
>>> f = sp.Function('f'); cand = x**2 * f(y / x)       # x u_x + y u_y = 2u
>>> sp.simplify(x * cand.diff(x) + y * cand.diff(y) - 2 * cand)
0

Reading the output. For \(2u_x + 3u_y = 0\) the characteristics are the lines \(3x - 2y = \text {const}\), and the solution is an arbitrary function \(F\) of \(3x - 2y\). For \(u_x + u_y + u = x\), the solution is a particular solution \(x - 1\) plus an arbitrary function of the characteristic variable \(x - y\) times the decaying factor \(e^{-(x + y)/2}\). The last line confirms Example 1.3.3: \(x^2F(y/x)\) satisfies \(xu_x + yu_y = 2u\) for every \(F\).

MATLAB/Octave: the Symbolic Math Toolbox has no general first-order PDE solver; the characteristic ODEs are solved with dsolve.

Software 1.7.2 (Tracing characteristics numerically). For variable speed, the solution at \((x, t)\) is the initial value at the foot of the characteristic through it, found by integrating \(dx/d\tau = -c(x)\) backwards from \(x\) for time \(t\). For \(u_t + xu_x = 0\) with \(u(x, 0) = e^{-x^2}\), compared with the exact \(f(xe^{-t})\) of Example 1.5.1:

>>> import numpy as np
>>> from scipy.integrate import solve_ivp
>>> c = lambda x: x                        # u_t + x u_x = 0
>>> f0 = lambda x: np.exp(-x**2)           # u(x, 0)
>>> def u_char(xq, t):                     # follow the characteristic back
...     s = solve_ivp(lambda tau, X: [-c(X[0])], (0, t), [xq], rtol=1e-11)
...     return f0(s.y[0, -1])
...
>>> t = 1.0
>>> for xq in (0.5, 1.5, 3.0):
...     print(xq, round(u_char(xq, t), 6), round(f0(xq * np.exp(-t)), 6))
...
0.5 0.966732 0.966732
1.5 0.737489 0.737489
3.0 0.295816 0.295816

Reading the output. The numerical and exact values agree to six decimals at all three points. The same code works for any speed \(c(x)\), including those for which the characteristics have no closed form. This is the basis of semi-Lagrangian numerical methods in weather prediction, which follow characteristics back from the grid points at each time step.

MATLAB/Octave: [tt, X] = ode45(@(tau, X) -X, [0 1], xq); exp(-X(end)^2).

Software 1.7.3 (An age-structured population). With constant births \(B = 100\) a year, death rate \(\mu (a) = 0.02 + 0.001a\), and a uniform initial profile of \(60\) at every age, the density at four ages after \(0\), \(20\) and \(60\) years, from the characteristic solution of Section 1.6:

>>> B = 100.0                              # births per year, constant
>>> H = lambda a: 0.02 * a + 0.0005 * a**2  # integral of mu = 0.02 + 0.001a
>>> u0 = lambda a: 60.0                    # uniform initial age profile
>>> def density(a, t):                     # u_t + u_a = -mu(a) u
...     if a >= t:
...         return u0(a - t) * np.exp(-(H(a) - H(a - t)))   # alive at time 0
...     return B * np.exp(-H(a))                     # born after time 0
...
>>> ages = (0, 20, 40, 80)
>>> for t in (0, 20, 60):
...     print(t, [round(density(a, t), 2) for a in ages])
...
0 [60.0, 60.0, 60.0, 60.0]
20 [100.0, 32.93, 22.07, 9.92]
60 [100.0, 54.88, 20.19, 0.9]
>>> [round(B * np.exp(-H(a)), 2) for a in ages]   # the steady state
[100.0, 54.88, 20.19, 0.82]

Reading the output. At ages below \(t\) the density already equals the stable distribution \(Be^{-M(a)}\) (the last line), because those individuals were born after time \(0\). After \(20\) years only age \(0\) is on the birth branch among the ages shown; the age-\(20\) cohort was aged \(0\) at time \(0\), so it comes from the initial population, \(60e^{-M(20)} = 32.93\) rather than \(54.88\). After \(60\) years the ages \(0\), \(20\) and \(40\) all match the stable values. The age-\(80\) density, \(0.90\) after \(60\) years, still comes from the initial population, and differs from the stable value \(0.82\); it will reach it only when the last of the initial cohort has died. Information about the initial state travels along characteristics and leaves the system through the old ages.

MATLAB/Octave: the same piecewise formula as a function handle.