Software 10.4.1 (Transform pairs). Fourier transforms computed symbolically from the definition, with the convention \(\hat f(\xi ) = \int f(x)e^{-ix\xi }\,dx\):

>>> import sympy as sp
>>> x, xi = sp.symbols('x xi', real=True)
>>> a, t, k, s = sp.symbols('a t kappa s', positive=True)
>>> def FT(f):                     # f-hat(xi) = int f(x) e^{-i x xi} dx
...     return sp.simplify(sp.integrate(f * sp.exp(-sp.I*x*xi),
...                                     (x, -sp.oo, sp.oo)))
...
>>> FT(sp.exp(-a*x**2))                    # a Gaussian goes to a Gaussian
sqrt(pi)*exp(-xi**2/(4*a))/sqrt(a)
>>> sp.simplify(2 * sp.integrate(sp.exp(-a*x) * sp.cos(x*xi), (x, 0, sp.oo),
...                              conds='none'))   # e^{-a|x|}, even
2*a/(a**2 + xi**2)
>>> heat = sp.exp(-x**2 / (4*k*t)) / sp.sqrt(4*sp.pi*k*t)
>>> FT(heat)                       # the heat kernel -> e^{-kappa xi^2 t}
exp(-kappa*t*xi**2)

Reading the output. The Gaussian \(e^{-ax^2}\) goes to \(\sqrt {\pi /a}\,e^{-\xi ^2/4a}\), the two-sided exponential to \(2a/(a^2 + \xi ^2)\), and the heat kernel to \(e^{-\kappa \xi ^2t}\), the multiplier that solves the heat equation in frequency.

MATLAB/Octave: with the Symbolic Math Toolbox, fourier(exp(-a*x^2)), which uses the same sign convention.

Software 10.4.2 (Spectral solution by the FFT). On a large periodic box the fast Fourier transform computes \(\hat g\) at the wavenumbers of the box, and each mode is advanced exactly by its multiplier. First the heat equation from \(e^{-x^2}\) to \(t = 2\), against the exact solution \(e^{-x^2/9}/3\); then \(u_t + u_{xxx} = 0\), where the box is \(4000\) long so that no part of the wave train wraps around in the time computed; the columns are \(t\), \(\int u^2\), the position and value of the maximum, and the exact value there from the Airy-kernel convolution:

>>> import numpy as np
>>> def grid(L, N):                        # periodic box and its wavenumbers
...     xs = (np.arange(N) - N//2) * L / N
...     return xs, 2 * np.pi * np.fft.fftfreq(N, d=L/N)
...
>>> def evolve(g, q, symbol, t):   # multiply g-hat by exp(symbol(xi) t)
...     return np.fft.ifft(np.fft.fft(g) * np.exp(symbol(q) * t)).real
...
>>> xs, q = grid(80.0, 4096)
>>> u = evolve(np.exp(-xs**2), q, lambda q: -q**2, 2.0)    # heat, kappa = 1
>>> float(np.abs(u - np.exp(-xs**2 / 9) / 3).max())        # exact at t = 2
8.326672684688674e-17
>>> from scipy.special import airy
>>> from scipy.integrate import quad
>>> def airy_exact(X, T):          # kernel (3T)^(-1/3) Ai(x / (3T)^(1/3))
...     c = (3 * T)**(1/3)
...     K = lambda y: airy((X - y) / c)[0] / c * np.exp(-y**2)
...     return quad(K, -8, 8, limit=200)[0]
...
>>> xs, q = grid(4000.0, 2**16)            # u_t + u_xxx = 0: symbol i xi^3
>>> for T in (1.0, 8.0, 27.0):
...     v = evolve(np.exp(-xs**2), q, lambda q: 1j * q**3, T)
...     i = np.argmax(v); norm = np.sum(v**2) * (xs[1] - xs[0])
...     print(T, round(norm, 6), round(xs[i], 3), round(v[i], 6),
...           round(airy_exact(xs[i], T), 6))
...
1.0 1.253314 -1.282 0.586298 0.586298
8.0 1.253314 -2.869 0.319361 0.319361
27.0 1.253314 -4.333 0.216484 0.216484
>>> np.sqrt(np.pi) * airy(-1.0188)[0] / 3**(1/3)   # limit of max(u) T^(1/3)
0.6582957084826466

Reading the output. For the heat equation the spectral solution is exact to rounding (\(8 \times 10^{-17}\)), because the multiplier is exact for every mode and the Gaussian is negligible at the edges of the box. For the dispersive equation, \(\int u^2 = 1.253314 = \sqrt {\pi /2}\) at every time, the maximum moves left, and its value agrees with the quadrature of the Airy convolution to six digits. Times \(t^{1/3}\), the maxima are \(0.586\), \(0.639\) and \(0.649\), increasing towards the limit \(\sqrt \pi \,\max \operatorname {Ai}/3^{1/3} = 0.658\), which they cannot exceed since \(g \geq 0\).

MATLAB/Octave: real(ifft(fft(g).*exp(1i*q.^3*T))) with q the wavenumbers.

Software 10.4.3 (Inverting a Laplace transform). The transform \(e^{-x\sqrt s}/s\) of Example 10.3.1 (with \(\kappa = 1\)) inverted symbolically, and then numerically by Talbot’s contour at \(x = 1\) and three times, against \(\operatorname {erfc}(1/2\sqrt t)\):

>>> x0 = sp.symbols('x', positive=True)
>>> sp.inverse_laplace_transform(sp.exp(-x0*sp.sqrt(s)) / s, s, t)
1 - erf(x/(2*sqrt(t)))
>>> import mpmath as mp
>>> f = lambda S: mp.exp(-mp.sqrt(S)) / S  # u(1, t) for kappa = 1
>>> for T in (0.1, 1, 10):
...     print(T, mp.nstr(mp.invertlaplace(f, T, method='talbot'), 10),
...           mp.nstr(mp.erfc(1 / (2*mp.sqrt(T))), 10))
...
0.1 0.02534731868 0.02534731868
1 0.4795001222 0.4795001222
10 0.8230632738 0.8230632738

Reading the output. SymPy returns \(1 - \operatorname {erf}(x/2\sqrt t) = \operatorname {erfc}(x/2\sqrt t)\). The numerical inversion agrees with it to all ten digits shown: at \(x = 1\) the temperature is \(0.025\) after \(0.1\) time units, \(0.48\) after \(1\), and \(0.82\) after \(10\).

MATLAB/Octave: ilaplace(exp(-x*sqrt(s))/s) with the Symbolic Math Toolbox.