first.c

Approximation of the first derivative of a function (1st version)

Approximation of the first derivative of a function (1st version)

Given a function f(x)=\cos(2\pi x), compute an approximation of its first derivative.

We will discretise the function at regular intervals on the segment x\in [0:1].

#include "grid/cartesian1D.h"

The main() function is the starting point of any C program.

int main() {

We declare a new (integer) variable which holds the number of grid points.

  int n = 10;

We create the grid.

  init_grid (n);

We allocate a new scalar field f…

  scalar f[];

initialise it…

  foreach()
    f[] = cos(2.*pi*x);

…and apply boundary conditions.

  boundary ({f});

Our discretisation now looks like this

set term pngcairo size 600,200
set xtics 0,0.1,10
set grid xtics
f(x)=cos(2.*pi*x)
unset key
plot [-0.05:1.05] f(x), "< seq 0.05 0.1 0.95" u 1:(f($1)) pt 7, \
                        "< seq -0.05 1.1 1.05" u 1:(f($1)) pt 7

(script)

where the green values (at the center of the discretisation cell) were set by the foreach() loop and the blue “ghost” values were set by the call to boundary(). By default, the boundary conditions are “symmetry” i.e. f'(0) = f'(1) = 0.

We will store the first derivative in field df.

  scalar df[];

To compute the first derivative, we will use a local stencil, for example

set label 'f[]' at 0.25,f(0.25)+0.2 center
set label 'f[1,0]' at 0.37,f(0.35)+0.2 center
set label 'f[-1,0]' at 0.13,f(0.15)-0.2 center
plot [-0.05:1.05] f(x), "< seq 0.05 0.1 0.95" u 1:(f($1)) pt 7, \
                        "< seq -0.05 1.1 1.05" u 1:(f($1)) pt 7

(script)

We can then write an approximation of the first derivative for each cell as

  foreach()
    df[] = (f[] - f[-1,0])/Delta;

where \Delta is the grid spacing (1/10 for the moment). Of course this means that the derivative is defined at the mid-point between f[] and f[-1,0] (0.2 in the example above). Given that the foreach() loop only iterates over the green points, this also means that the derivative is not defined for x=1.

Note that we could also have defined df[] using f[] and f[1,0], it is just a matter of choice.

To check whether this works, we can write the values to standard output which is redirected to a file called ‘out’.

  foreach()
    printf ("%g %g\n", x - Delta/2., df[]);

Note that the x coordinate in the foreach() loop is the position of the center of the cell (the green dots), so that we need to substract \Delta/2 to get the position of the derivative.

After the program has finished we can plot both the points in ‘out’ and the exact value of the derivative -2\pi\sin(2\pi x).

unset label
plot [0:1]-2.*pi*sin(2*pi*x), 'out' pt 7 lt 4

(script)

Although it looks like the error is small, we need to measure it and do a proper convergence study. This will be the next exercise.