Software 6.4.1 (The heat kernel). Symbolic checks of the kernel and of the solution of Example 6.1.2:

>>> import sympy as sp
>>> x, y = sp.symbols('x y', real=True)
>>> t, k, a = sp.symbols('t kappa a', positive=True)
>>> Phi = sp.exp(-x**2 / (4*k*t)) / sp.sqrt(4*sp.pi*k*t)   # the heat kernel
>>> sp.simplify(Phi.diff(t) - k * Phi.diff(x, 2))
0
>>> sp.integrate(Phi, (x, -sp.oo, sp.oo))
1
>>> s = sp.sqrt(4*k*t)
>>> u = (sp.erf((x + a) / s) - sp.erf((x - a) / s)) / 2    # data 1 on |x| < a
>>> sp.simplify(u.diff(t) - k * u.diff(x, 2))
0
>>> sp.limit(u.subs(x, 0), t, 0), sp.limit(u.subs(x, 3*a), t, 0)
(1, 0)

Reading the output. The kernel satisfies the equation and has total mass \(1\). The error-function solution satisfies the equation for every \(a\) and \(\kappa \); as \(t \to 0\) it tends to \(1\) inside the box and to \(0\) outside.

MATLAB/Octave: 0.5*(erf((x+a)./sqrt(4*k*t)) - erf((x-a)./sqrt(4*k*t))).

Software 6.4.2 (Infinite speed, in numbers). The solution of Example 6.1.2 with \(\kappa = a = 1\) at \(x = 3\), two units outside the box, at four early times, computed with \(30\)-digit arithmetic and with the complementary error function so that the tiny values are not lost to cancellation; then the decay of the centre value:

>>> import mpmath as mp
>>> mp.mp.dps = 30
>>> def box(x, t):              # kappa = a = 1, written with erfc
...     s = mp.sqrt(4 * t)                 # so that tiny values are not lost
...     return (mp.erfc((x - 1) / s) - mp.erfc((x + 1) / s)) / 2
...
>>> for t in (1e-3, 1e-2, 1e-1, 1.0):      # at x = 3, outside the box
...     print(t, mp.nstr(box(3, t), 5))
...
0.001 4.5258e-437
0.01 1.0442e-45
0.1 3.8721e-6
1.0 0.076311
>>> [mp.nstr(box(0, t) * mp.sqrt(mp.pi * t), 6) for t in (10, 100, 1000)]
['0.991729', '0.999167', '0.999917']

Reading the output. At \(t = 0.001\) the temperature at \(x = 3\) is positive, but \(4.5 \times 10^{-437}\); at \(t = 0.01\) it is \(10^{-45}\). Only when the diffusion length \(\sqrt {4t}\) becomes comparable with the distance \(2\) does the heat arrive in measurable amounts: \(0.076\) at \(t = 1\). The last line shows \(u(0, t)\sqrt {\pi t} \to 1\), the asymptotic \(u(0, t) \approx 2a/\sqrt {4\pi t}\) of Example 6.1.2.

MATLAB/Octave: 0.5*(erfc((x-1)./sqrt(4*t)) - erfc((x+1)./sqrt(4*t))) (double precision underflows below \(10^{-308}\)).

Software 6.4.3 (Temperature waves in the ground). For moist soil with \(\kappa = 5 \times 10^{-7}\) m\(^2\)/s, the depths at which the daily and the annual temperature waves are reduced to \(10\%\) of their surface amplitude (\(x = \ln 10/q\)) and are exactly out of phase with the surface (\(x = \pi /q\)), in metres:

>>> import numpy as np
>>> kappa = 5e-7                           # m^2/s, moist soil
>>> for name, period in (('day', 86400), ('year', 365.25 * 86400)):
...     w = 2 * np.pi / period
...     q = np.sqrt(w / (2 * kappa))       # u = A exp(-q x) cos(w t - q x)
...     print(name, round(np.log(10) / q, 3), round(np.pi / q, 3))
...
day 0.27 0.368
year 5.16 7.041

Reading the output. The daily cycle is reduced to \(10\%\) at \(27\) cm and reversed at \(37\) cm. The annual cycle, \(\sqrt {365.25} = 19.1\) times slower in \(q\), reaches \(5.2\) m before falling to \(10\%\), and at \(7.0\) m the ground is warmest in winter. At a depth of \(1\) m the daily cycle has vanished (a factor \(e^{-8.5}\)), while the annual one keeps \(e^{-0.446} = 64\%\) of its amplitude and lags the surface by \(0.446\) radians, about \(26\) days.

MATLAB/Octave: q = sqrt(2*pi/86400/(2*5e-7)); [log(10) pi]/q.