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}]]