Software 3.4.1 (A canonical form). The discriminant and characteristic slopes of \(u_{xx} + 4u_{xy} + 3u_{yy} = 0\), and a check of the general solution of Example 3.2.2:

>>> import sympy as sp
>>> x, y = sp.symbols('x y')
>>> a, b, c = 1, 2, 3                      # u_xx + 4 u_xy + 3 u_yy = 0
>>> b**2 - a*c                             # discriminant: hyperbolic
1
>>> m = sp.symbols('m')
>>> sp.solve(a*m**2 - 2*b*m + c, m)        # characteristic slopes dy/dx
[1, 3]
>>> F, G = sp.Function('F'), sp.Function('G')
>>> u = F(y - x) + G(y - 3*x)              # constant on each family
>>> sp.simplify(u.diff(x, 2) + 4*u.diff(x, y) + 3*u.diff(y, 2))
0

Reading the output. The discriminant is \(1 > 0\), the slopes are \(1\) and \(3\), and \(F(y - x) + G(y - 3x)\) satisfies the equation identically for arbitrary twice-differentiable \(F\) and \(G\).

MATLAB/Octave: roots([1 -4 3]).

Software 3.4.2 (Hadamard’s example). The size of the Cauchy data and of the solution at \(y = 1\) in Example 3.3.1:

>>> import numpy as np
>>> for n in (1, 5, 10, 20, 40):           # Hadamard: u_xx + u_yy = 0
...     data = 1 / n**2                    # max |u_y(x, 0)|, with u(x, 0) = 0
...     sol = np.sinh(n) / n**3            # max |u(x, 1)|
...     print(n, f"{data:.2e}", f"{sol:.2e}")
...
1 1.00e+00 1.18e+00
5 4.00e-02 5.94e-01
10 1.00e-02 1.10e+01
20 2.50e-03 3.03e+04
40 6.25e-04 1.84e+12

Reading the output. As \(n\) goes from \(1\) to \(40\) the data shrink by a factor of \(1600\), while the solution at \(y = 1\) grows from \(1.18\) to \(1.84 \times 10^{12}\). A perturbation of the data smaller than any measuring instrument could detect changes the solution beyond recognition.

MATLAB/Octave: n = [1 5 10 20 40]; [1./n.^2; sinh(n)./n.^3].

Software 3.4.3 (The backward heat equation). An error of \(10^{-10}\) in each of the first \(40\) sine coefficients of a final temperature profile, propagated backwards; the root sum of squares of the coefficient errors at three earlier times:

>>> k = np.arange(1, 41)                   # backward heat: u_t = -u_xx
>>> eps = 1e-10                            # error in each sine coefficient
>>> for t in (0.01, 0.05, 0.1):
...     print(t, f"{np.sqrt(np.sum((eps * np.exp(k**2 * t))**2)):.3e}")
...
0.01 9.987e-04
0.05 5.542e+24
0.1 3.070e+59
>>> eps * np.exp(-k[-1]**2 * 0.1)          # the forward problem damps it
3.257488532207521e-80

Reading the output. Going back only \(0.01\) time units turns errors of \(10^{-10}\) into errors of \(10^{-3}\), and going back \(0.1\) produces \(3 \times 10^{59}\): the computed initial state is meaningless. The same error run forwards for \(0.1\) units is multiplied by \(e^{-160}\) in the \(40\)th mode; the forward problem is as stable as the backward one is unstable.

MATLAB/Octave: k = 1:40; norm(1e-10*exp(k.^2*0.05)).