In Mathematica, a *module* is an extended function that may perform several steps. An advantage of a module is that it may declare and define *local variables* which can be used inside the module but not in the global scope. It is generally good programming practice to keep variables inside as local a scope as possible, because variables in global scope can conflict with each other if you are not careful to keep track of their names.

## A module for solving the KdV equation

Below is an example usage of a module. The module `solat` computes the solution to the KdV equation

at a desired time `tm`, given initial data `data0`. It uses a split-step pseudospectral integrator, whose details, as well as an equivalent code in Matlab, can be found at Wikiwaves. (Note that their Matlab code makes reference to , the angular wavenumber, whereas this code uses , which is equal to .

nx = 128; xmax = 1; dx = 1/nx; xs = Range[0, xmax - dx, dx]; numax = 1/(2*dx); dnu = 1/xmax; nus = Join[Range[0, numax - dnu, dnu], Range[-numax, -dnu, dnu]]; u0[x_] := Sin[6 Pi x]; data = u0 /@ xs; solat[tmax_, data_] := Module[{t = 0, dft = Fourier[data], dt = 0.4/nx^2, dft1 = 0, sol = 0}, While[t < tmax, dt = Min[dt, tmax - t]; dft1 = dft*Exp[I (2 Pi nus)^3 dt]; dft = dft1 - 3 I 2 Pi nus dt Fourier[InverseFourier[dft1]^2]; t += dt; ]; sol = Re[InverseFourier[dft]]; sol ]

We declare the variables `t`, `dft`, `dt`, `dft1` and `sol` as local to the module, and initialise them to sensible values. (Note that `dft1` and `sol` has been initialised to `0`; this is completely arbitrary, because this value will be overwritten a few lines later.)

The last line inside the module, `sol`, simply causes `sol` to be returned from the module.

## A module for solving the sine-Gordon equation

We now consider the sine-Gordon equation

Like KdV, the sine-Gordon equation is quasilinear, in the sense that no term contains any nonlinear function of a derivative. Moreover, in both equations there is only one nonlinear term and no inhomogeneities; the rest are linear and may be treated easily. From a computational perspective, the most important difference between the two is that the sine-Gordon equation is second-order in time. We first write it as a pair of coupled first-order equations:

This system can now be solved using a split-step pseudospectral method, much like the above, except that we must keep track of both dynamical variables and $\latex \psi$. The derivation of this is a bit messy, but follows the same principle of evolving the linear parts of the equations separately from the nonlinear part. Here is the Mathematica code:

nx = 128; xmax = 1; dx = 1/nx; xs = Range[0, xmax - dx, dx]; numax = 1/(2*dx); dnu = 1/xmax; nus = Join[Range[0, numax - dnu, dnu], Range[-numax, -dnu, dnu]]; u0[x_] := Sin[6 Pi x]; data = u0 /@ xs; sinegordonsol[tmax_, phi0_, psi0_] := Module[ {t = 0, phift = Fourier[phi0], psift = Fourier[psi0], dt = 0.4/nx^2, phift1 = 0, psift1 = 0, sol = 0}, While[t < tmax, dt = Min[dt, tmax - t]; phift1 = phift Cos[dt] + psift Sin[dt]; psift1 = psift Cos[dt] - phift Sin[dt]; phift = phift1; psift = psift1 - Fourier[Sin[InverseFourier[phift]]] dt; t += dt; ]; sol = Re[InverseFourier[phift]]; sol ]

To use it, do something like

sgs = sinegordonsol[1, data, 0*data]; ListLinePlot[Transpose[{xs, sgs}]]