Prophecy

I am disappointed at the outcome of the EU referendum, as well as the way the campaign was conducted, for several reasons, which need not be detailed here.

But I have been cheered up, after my attention was drawn to this speech by Nick Clegg in 2014:

That’s what we are fighting for. So what are we fighting against? Imagine again what it will be like in 2020, but this time with the Conservatives in Government on their own. Britain, diminished and divided after a botched attempt to renegotiate our relationship with Europe and a vote to withdraw from the European Union. Companies pulling out of the UK left, right and centre, the markets losing confidence, hiking up our borrowing costs and halting the recovery in its tracks. Workers fearing for their jobs, not just because the companies they work for are plunged into uncertainty but because their bosses can fire them at will, no questions asked. The young and the working poor hit time and time again as George Osborne takes his axe to the welfare budget with no regard for the impact on people’s lives. Schools run in the interests of profit for shareholders rather than the life chances of their pupils. A Home Office state snooping on your emails and social media. Opportunity reserved for a few at the top and everyone else told to make do with what they’ve got. A Tory party leadership in hock to their right wing, desperately running after and pandering to UKIP’s ugly nationalism. A Prime Minister trapped between being a poor man’s Margaret Thatcher and a rich man’s Nigel Farage. “Compassionate Conservatism” just a sound bite from a bygone age.

How prescient.

Solving ODEs in C

I set my students an exercise: to solve the nonlinear differential equation

\ddot{x} + \sin(x) = A\sin(\omega t).

What follows is an outline of a model solution (but by no means complete, or optimal). I take this opportunity to introduce several other features of C that I haven’t had the chance to talk about in the classes.

This equation is, of course, the (nondimensionalised) equation of motion of a forced pendulum, where x represents the angle that the pendulum makes with the downward vertical. If A = 0 or \omega = 0, then the pendulum is unforced, and executes regular oscillations. The period of these oscillations is approximately equal to 2\pi for small-amplitude oscillations (where the motion is approximately simple harmonic motion), but increases as the amplitude increases. If forcing is turned on, then a second period—that of the forcing—is introduced. The pendulum may still oscillate with its ‘intrinsic’ period, or if the forcing is strong enough then the pendulum may start oscillate with the period of the forcing. If the period of the forcing is close to the intrinsic period of the pendulum, then resonance may occur.

We first write the equation as a system of two equations,

\dot{x} = y \text{ and } \dot{y} = - \sin(x) + \sin(\omega t)

which is first-order. We then solve this system numerically using the Runge-Kutta method. The program diffsolve is meant to do this using a Runge-Kutta solver. Its full source code can be downloaded here. The code, as it is written, is quite general; it can solve any system of two ODEs of the form

\dot{x}=f(t,x,y) \text{ and } \dot{y}=g(t,x,y).

Annotation of the code

After the initial comments and the #includes, we first come across a typedef struct:

typedef struct {
    double amp;
    double omega;
} problem_t;

I declare a type, problem_t, which will store the parameters of the problem, amp and omega. I usually wrap all parameters into a struct like this, because it lets me pass a single variable, of type problem_t, to a function, instead of passing each of the individual parameters. It also means that, if I want to change my program to allow more parameters, I can simply change the declaration of problem_t, rather than changing the definition of every function that needs these parameters. This pays off when we work with problems with ten or twenty parameters.

The first function defined makes use of this:

double forcing(double t, problem_t params) 
{
    return params.amp * sin( params.omega * t );
}

This is used by the function rhs_y, which calculates \dot{y} . (You could do away with the function forcing and have the above inside rhs_y, if you prefer.)

The function rhs_x, which gives the derivative \dot{x} , appears to require arguments that it does not use:

double rhs_x(double t, double x, double y, problem_t params)
{
    return y;
}

Although the arguments t, x and params are not actually used, I have chosen to include them anyway, because in the future I might want to move from solving \dot{x} = y to something more general, of the form \dot{x} = f(t,x,y) .

Now we come to the main function. The first thing we notice is that the function prototype says

int main(int argc, char* argv[]) 

instead of simply

int main()

We want to be able to run our program without having to recompile the program each time we want to change our parameters such as initial conditions. We will type something like

$ ./diffsolve 1 2 100 0.005 1 0

into the command line, in order to specify the different parameters to diffsolve. These are command line arguments, and they are treated by C as arguments to main. The argument argc is the argument count, the number of arguments; and argv is the argument vector, an array of strings. Note that the program’s name, ./diffsolve, is treated as the 0th argument, and so argc is always at least equal to 1.

The arguments are actually read in the following block:

if (argc > 6) 
{
    params.amp   = atof(argv[1]);
    params.omega = atof(argv[2]);
    tmax  = atof(argv[3]);
    dt    = atof(argv[4]);
    x     = atof(argv[5]);
    y     = atof(argv[6]);
}
else
{
    fprintf(stderr, "Usage: %s amp omega tmax dt x0 y0\n", argv[0]);
    exit(-1);
}

If the user has specified at least seven arguments (where the 0th argument is the name of the program), then these arguments are taken. The atof function converts strings to floats (although our variables are declared as doubles; this is fine, because floats can be implicitly cast into doubles).

If the user has not specified enough arguments, then the program gives an error which tells the user how to use it, and exits with a return code of -1, which (by convention) indicates an error. (The exit function is provided in stdlib.h.) The error message is printed using fprintf, which stands for ‘printf to a file’, or ‘print formatted text to a file’. Here the ‘file’ is stderr, the standard error stream. In contrast, the function printf writes to the standard output stream, represented by stdout. By default, both stdout and stderr write to the screen, but the two streams may be separated from each other. It is a good idea to use stderr for errors, warnings or diagnostic messages, and stdout for actual output that you want to process.

The rest of the program is fairly straightforward: the values of t, x and y are updated in a while loop, and at each time step their values are printed (to standard output).

Usage (on Unix-like systems)

After compiling, the command

$ ./diffsolve 1 2 100 0.005 1 0

will print plenty of lines, giving the solution of the system at each timestep.

You probably want to save this information to another file, which you can do by writing

$ ./diffsolve 1 2 100 0.005 1 0 > out.txt

This will write the results to the file out.txt, which you can then process somewhere else. For example, you can use gnuplot to plot the results:

$ gnuplot
gnuplot> set terminal png
gnuplot> set output "out.png"
gnuplot> plot "out.txt" using 1:2 with lines

This produces the image out.png:
out.png

Mathematica modules

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

u_t + 6uu_x + u_{xxx} = 0

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 k, the angular wavenumber, whereas this code uses \nu, which is equal to k / (2\pi) .

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

\phi_{tt} - \phi_{xx} + \sin\phi = 0.

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:

\phi_t = \psi \text{ and } \psi_t = \phi_{xx} - \sin\phi.

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

The discrete Fourier transform in Mathematica

The Fourier transform is a very useful and widely used tool in many branches of maths. The continuous Fourier transform, available in Mathematica as FourierTransform, tries to perform the Fourier transform of a symbolic expression. This is useful for solving many problems if the Fourier transform can also be expressed in a nice symbolic form, but such a form is not necessarily available (and indeed the Fourier transform might not exist at all).

In practice, the discrete Fourier transform, which works with a finite set of data points and approximates the continuous Fourier transform, is more important. The DFT is available in Mathematica as Fourier, and the inverse DFT as InverseFourier. Note that Mathematica uses the 1/\sqrt{n} normalisation by default for both the DFT and the IDFT.

Example

Here is an example usage of the DFT in Mathematica:

xmax = 1;
nx = 128;
dx = xmax / nx;
xs = Range[0, xmax - dx, dx];
f[x_] := Exp[Sin[6 Pi x]] + Sin[40 Pi x];
data = Map[f, xs];
ListLinePlot[Transpose[{xs, data}]]

dft = Fourier[data];
numax = 1/(2*dx);
dnu = 1/xmax;
nus = Join[Range[0, numax - dnu, dnu], Range[-numax, -dnu, dnu]];
ListPlot[
Transpose[{nus, Abs[dft]}],
PlotRange -> All]

The first block generates a list data which consists of sampled values of the function f at the values xs, which are 128 values between 0 and 1 (not including the endpoint at 1). These sampled points are plotted.

The second plot calculates the discrete Fourier transform of data (according to the default normalisation). Then the list nus is constructed: these are the wavenumbers that correspond to each element of dft. (With the default normalisation, this is not the angular wavenumber: the wavenumber of \sin(40\pi x) is said to be 20, while the angular wavenumber is 40\pi. The symbol k is often used for the angular wavenumber; I am using \nu for the wavenumber.) We then plot the absolute values of dft against the wavenumbers nus.

Aliasing

Note how the values of \nu_{\max} and \Delta\nu are related to the values of \Delta x and x_{\max}:

\nu_{\max} = \frac{1}{2\Delta x} \text{ and } \Delta\nu = \frac{1}{x_{\max}}.

What this means is the following:

  • If you sample more densely (so \Delta x is smaller), then you can pick up higher frequencies.
  • If you sample for a longer domain (so x_{\max} is larger), then you can pick up lower frequencies.

The highest frequency that you can pick up, $1/2\Delta x$, is called the Nyquist frequency and is equal to half the sampling rate of the data. If your function in fact oscillates at a higher frequency than this, then your sampled data will not pick up that component properly, and many of your results will be inaccurate. This effect is known as aliasing.

Applications

The Fourier transform can be used to solve many different problems, and is particularly useful for solving linear partial differential equations. Numerical PDE solvers that make use of the DFT are part of a class of numerical schemes known as spectral methods. For linear problems where the solution is oscillatory or where a ‘conservation law’ (usually analogous to conservation of mass or momentum) applies, spectral methods tend to be better than finite differences. The DFT can also be useful for quasilinear PDEs, using a family of methods known as pseudospectral methods.

Resources for the C course

For the benefit of my students, here are the codes that I have used in the lessons. The contents of that directory will be updated periodically.

A useful set of tutorials can be found at Programiz or at W3schools.

HackerRank provides a number of challenges or exercises for your programming skills, but you may need to get a bit more familiarity with the language before being able to complete them.

Why write in C?

Having been asked the question ‘Why would you write in C?’, I write here a few words on that. A Google search would recommend articles by people with much more expertise than me, who can go into much greater detail. The book 21st Century C may be useful as well.

Here are some of the reasons that I like C:

  • C is free, in the sense of free speech as well as free beer. The C standard is open, free implementations of C compilers (such as GCC) exist widely, and many of the libraries available are free and open source. You can distribute C code to other people, and know that they will be able to compile and run it, provided that they have a suitable platform. (In contrast, your colleague will only be able to run your Matlab programs if they too have a license for Matlab.)
  • C is a small language. The basic grammar of C consists of only a few keywords, and the standard library provides only a handful of mostly self-explanatory functions.
  • C is portable. A corollary of C’s smallness is that that a C compiler can be implemented relatively simply on many different systems. This is why C is often used for firmware, especially on small or specialised devices such as games consoles.
  • C is well-designed. For example, C has straightforward scoping rules, whereas Matlab (and other scripting languages such as PHP) tend to put almost all variables into a global scope, which increases the likelihood of ‘cross-contamination’ between scripts or different parts of a script. Matlab also has rules on naming conventions or on syntax which are unnecessarily strict; C is not so strict with these petty rules, and allows you to use the right style for each job (although it is up to you as the programmer to decide what that style is).
  • C evolves slowly. Since 1989, changes to standard C have occurred approximately once every ten years; these changes have been very incremental and almost always backwards-compatible. (That said, new libraries are written all the time, and programming styles have changed.) Contrast this with Matlab, where new versions are released twice a year and code written in older versions often fail with newer versions.
  • C is the parent of many other languages. The syntax of C has been inherited by languages such as Java, Python and PHP, and in turn by their descendants. For example, although functions for printing formatted text existed in languages before C, it is C’s printf function whose name and syntax is used by most languages.
  • C is widely used. Because of its age and its portability, C remains a standard language amongst systems programmers. Because C is relatively low-level, allowing (and often mandating) direct control over memory usage, it is used for firmware as well as things such as the Linux kernel. That said, C’s usefulness is not restricted to systems programming; the wide range of libraries available means that it is also used for applications or scientific computing.

On the other hand, C has some downsides. C is hard to learn compared to Matlab, and its low-level nature means that you have to be careful about things like memory management (which is done automatically by languages such as Java).

Creating user interfaces, graphics or plots with C is not immediately easy. All the libraries are there, but you will still have to do quite a lot of work. I usually prefer to have my C applications output data into text files, which I then give to gnuplot or Matlab if I want a plot out of them.

I wish I could give a comparison against Fortran or Python as well as Matlab, as those are other languages that are often used for scientific computing. But I don’t know either of them.

Building C: gcc and make

As discussed in Basic principles of C, code written in C needs to be built (compiled and linked) to produce executables. What follows is a short discussion of this process.

Building with GCC

The GNU Compiler Collection (GCC) is a suite of tools for the compilation process. The GCC is free software. It is pre-installed on most GNU/Linux distributions, and is also available on OS X and Windows. (On Windows you will need to install MinGW, a development environment that emulates, and gives the functionalities of, a basic Unix-like system.)

Due to its popularity and its free and open-source nature, the GCC has become, officially or unofficially, the standard compilation system on most Unix-like systems.

The GCC is designed so that the interface is the same on all systems. To compile the C source file foo.c, the basic command is
$ gcc foo.c
This will compile and link foo.c and produce the executable file a.out. You can specify a different name for the output file using the -o flag:
$ gcc -o foo foo.c

Multiple source files

If your program is written across several different source files, then the linking process can be taken care of automatically:
$ gcc -o output foo.c bar.c baz.c

Options

There are many options that can be passed to gcc (and to its sisters, described below). Here are some of the more important ones:

  • -I/path/to/dir: Adds the specified directory to the list of places where #include commands will search for files. On Unix-like systems, directories such as /usr/include and /usr/local/include are usually searched automatically, but if you need the header file of a library that is not installed there, then you will need to specify the path.

    For example, on my laptop (running OS X), the GSL is installed in /usr/local/Cellar/gsl/1.16. To use the GSL headers, I would need to pass the flag -I/usr/local/Cellar/gsl/1.16/include/.
  • -llibrary: Indicates that library should be linked. For example, if I am using the GSL, then I would need to link using -lgsl.
  • -Wall: Tells the compiler to emit all warnings, not just fatal errors or serious warnings. This can be useful for revealing those places where your code works, but is not written as cleanly as it could be (for example, if you have declared variables which you do not actually use).
  • -g: Tells the compiler to include debugging information in the produced object code or executable. For example, line references or even the actual text of the source file might be included. This is useful because it allows a debugger such as gdb or valgrind to tell you not just that memory leaks are happening, but also tell you where they take place.
  • -O0: Tells the compiler not to optimise the compiled code. Usually, compilers do not simply translate source code directly into machine code; they may, for example, change the order in which some statements are executed. This improves the speed of the program and may reduce the size of the executable. Turning off optimisation may be desirable if you want a faster compilation, and is also useful in testing because the produced executable will better resemble the source in its structure.

Other tools in the GCC

The GCC also provides a suite of compilers for other languages, including g++ for C++ and gfortran for Fortran. All of these generate object code files (files with the .o extension) which are compatible: they may be linked together. Hence different parts of a program may be written in different languages, compiled into object code and then linked into one executable.

GNU Make

If you are working on a project with many source files and you have to compile many times during development, then typing out the above commands over and over again, and trying to remember what flags to pass to gcc, can be tedious. GNU Make is a very useful tool for automating the build process.

The GNU Make utility is not part of the C language, but it is an important part of its ‘culture’: Software written in C (or C++, or Fortran, or any other compiled language) is often built using GNU Make. Like GCC, GNU Make usually comes with GNU/Linux distributions and is available on OS X and Windows (with MinGW).

GNU Make will be described in more detail in a future post.

Other tools

I mainly stick to using GNU Make, but it is good to mention a few other tools: also not parts of the C language, but common parts of the culture.

Git

Git is a version control system that allows you to keep track of changes to source files in a project. You edit your files and then commit your changes each time you make a change. The history of these commits is recorded, and you can revert to previous versions.

Git is designed for multiple developers working on different parts of a project. Alice, having changed a few files, commits and then pushes her commits. Bob can then pull these changes: Git will try to update Bob’s copy of the code to apply the changes that Alice has made, without overwriting any changes that Bob might have made.

Alice and Bob may also create branches of the code if they want to work on different features of the project whilst keeping the master branch unchanged.

GNU Autotools

If you are working on a large project, then building can be tedious even with the aid of GNU Make: Keeping track of your source files and updating your makefiles to reflect their changes is difficult. The GNU Autotools suite is designed to help automate that process.

Integrated development environments (IDEs)

Alternatively, you may use an IDE such as Eclipse (cross-platform), Xcode (on OS X), or Microsoft Visual Studio (on Windows) to keep track of files and build processes. Most IDEs provide a compilation and linking system as well as a debugger. Their interfaces tend to be daunting, however.

Basic principles of C

History of C

C has existed since 1972. Its current standard is C11 (published in 2011), but only a few details have changed since the 1989 standard (C89, also known as ANSI C). It is an extremely standardised language; different operating systems provide different libraries with different functions, but the syntax of the language is very simple. At its core, C is a very small language; there are only a handful of keywords offering basic notions such as basic types (e.g. int and float) and control structures (e.g. for and switch). The basic language is supplemented with the standard library which provides functions for basic operations, such as input and output, string handling and mathematical functions.

C is the precursor of many modern programming languages, and has influenced many others. Hence, if you are familiar with any of C++, C#, Java, Perl, PHP or Python, then you will be familiar with many aspects of C’s syntax.

Compilation and linking

C is a compiled language. Raw C code is stored in source files or header files. In order to produce an executable application, the source files must be compiled into object files, and then the object files must be linked. (Traditionally, source files have extension .c, header files have extension .h, and the object file produced by compiling foo.c is called foo.o.) The combined process of compiling and linking is called building. Sometimes, the term compiling is used to mean both compiling and linking.

Often, the same functions will need to be used by many different applications: in particular, the standard library functions. It would be wasteful for each C programmer to implement their own exp function, but exp (and other basic mathematical functions) are provided in the standard library. When writing your own program foo, all you need to do before you can use exp is to tell the compiler that exp is a function that takes in a floating-point number and returns another floating-point number. This is a declaration. For a standard library function such as exp, the declaration is included in a standard header file, here math.h, which you include into your source file foo.c:

#include"math.h"

Somewhere in math.h is a line that looks like

double exp(double);

which is the declaration of exp. You then compile foo.c and produce an object file foo.o. At this point, foo.o does not contain any details about what exp actually does; all it knows is what arguments exp takes in and what it returns. Hence you cannot run foo.o; you must link it with another object file, probably called something like math.o, where the actual implementation of exp is given.

There are two advantages to doing things this way. One is that the implementation of a function like exp can be used by millions of different programmers. The other is that your usage of a function is independent of its implementation; different systems or machines will often provide different implementations, but your part of the code will work with any of them.

In practice

The build processes require a number of commands. When building a large project, it is common to automate the build processes using a number of tools—not actually part of the C language, but important parts of its ‘culture’. GNU Make, which will be the subject of a future post, is one such tool; many projects distribute their source code along with a ‘makefile’ which contain instructions for automating the build process.