calculus

• Introduction
• The Problem
• The Model
• Implementation
• Examples
• Accuracy
• Usage
• Further Work

Introduction

In this model we'll be applying trad4 to the problem of ordinary differential equations, again using numerical methods in a similar vein to heat_equation_1d. We'll define a function, and get the first and second order differentials of that function.

The Problem

The diagram below shows the curve f and two lines a and b.

First, we're going to find the first differential at point x=2. This can be approximated using the centred difference method which gives us the slope of line a. This is given by:

df = ( f(x3) - f(x1) ) / ( x3 - x1)

The first differential at point x=0 can be approximated by the forwards difference method which gives us the slope of line b. The forward difference method is less accurate than the centred difference method for reasons discussed below, but for now, this is given by:

df = ( f(x1) - f(x0) ) / ( x1 - x0)

The first differential at the last point on the curve (not shown) is given by the backwards difference method. This is not show on the above diagram as it's essentially the same as the forwards difference method.

The Model

Now we're going to walk through how a 10-node implementation works. First we need to define our 10 f-objects. We start in the calculate_f function where we first derive our x-values from their object id, and these run from 0.1 to 1.0. Next the y-values are derived, and for this we're using the default y = sin(5x).

Next, we move on to the df objects. The published values of these objects are the first derivative of f (given analytically as y = 5cos(5x)), and are given by the foreward and centred difference methods described above:

Finally we add our third tier d2f. This applies the same forward and centred difference methods, but applied to df rather than f. This gives us our second derivative (given analytically by y = -5*5sin(5x)):

We can see this working using spot-checks. Take the full data set of this run:

First we'll look at the green cell. According to the centred difference formula, this is given by:

	df = ( 0.909297 - 0.841471 ) / ( 0.4 - 0.2 )
	df = 0.339130

Next,the yellow cell which uses the forward difference method, which means that d2f is given by:

	d2f = ( 2.59035 - 3.62045 ) / ( 0.2 - 0.1 )
	d2f = -10.301000

Lastly we'll check the blue cell, which uses the centred difference method again:

	d2f = ( -1.01061 - -4.48961 ) / ( 0.9 - 0.7 )
	d2f = 17.395000

The small differences you see are because the tabular output above is the result of 32-bit floating point calculations whereas we've taken the rounded output to plug into the three worked examples above.

Implementation

The abstract diagram

The concrete diagram

The monitor is not shown.

The t4 files

f.t4:

static
    int mode

pub
    double x
    double y

• static:
  • mode is used to easily switch between if ( mode == 0 ) y = sin(5x) else y = x*x
• pub:
  • x is the x-value of each node. As in heat_equation_1d this is derived from the object's id, and for all the models discussed in this document runs from 0.0-1.0, with various step-sizes.
  • y is the value of the function we're differentiating. In our stock example this is sin(5x)

The mode switch is just to allow you to switch between the two demo functions without having to recompile. Arguably it shouldn't exist on every f-object but on a singleton data_server. See the Usage section below for how to use this.

df.t4:

sub
    f f_down
    f f_this
    f f_up

pub
    double df
    double dx

• sub:
  • f_down is the f-object to the below left of this object
  • f_this is the f-object below this object
  • f_up is the f-object below right of this object
• pub:
  • df is the first order differential of the f_this object
  • dx is the step-size of the f_this object. This is a pass-through, discussed below.

The pub element dx is described as a pass-through above. This is required as d2f needs to know the step-size of the f of which it's finding the second differential. While it's theoretically possible for d2x to read x directly from f (by using another level of indirection of the kind used in the sub macros) this is not formally supported - any object only has access to the data of objects to which it directly subscribes. For this reason dx has to copy x and publish it to d2f unmolested. This is known as a pass-through.

d2f.t4:

sub
    df df_down
    df df_this
    df df_up

pub
    double d2f

• sub:
  • df_down is the df-object to the below left of this object
  • df_this is the df-object below this object
  • df_up is the df-object below right of this object
• pub:
  • d2f is the second differential of the f_this object

monitor.t4:

sub
    f my_f[NUM_NODES]
    df my_df[NUM_NODES]
    d2f my_d2f[NUM_NODES]

• sub:
  • my_f[NUM_NODES] holds the ids of all the f objects
  • my_df[NUM_NODES] holds the ids of all the df objects
  • my_d2f[NUM_NODES] holds the ids of all the d2f objects

Examples

For this section we're switching to a 100-node graph. This will provide smoother curves, but not too smooth as we'll see in the next section.

sin(5x)

The full 100-node data set is available here.

You'll notice a glitch at the end of the numeric d2f graph. This is down to the fact the forwards/backwards difference methods are not as accurate as the central difference method. This will be discussed in the next section.

x squared

The full 100-node data set is available here.

You'll see the same glitch at the start and end of the d2f graph - discussed further below.

Accuracy

The Edges

We're now going to look at where those glitches at the edges come from, why the forwards/backwards difference methods produce these artifacts and why there's nothing we can do about it.

Firstly, you might wonder why only d2f is affected. The answer is it's not - df has this artifact too, it's just not as visible when looking at the graph. Take a look at the right-hand end of the 100-node data for the sin(5x) function. You can see the df error in yellow, and the d2f errors in green and blue:

Take a look at the curve below. This looks at the start of an arbitrary graph but as already explained the problems at the start and the ends are identical. As above we have a curve f, and two lines a and b. We're trying to find the slope of the curve at x=0.

The actual method used - the forward difference method - corresponds to the slope of line a. If the imaginary x=-1 point existed we could use the centred difference method which would give us the slope of line b. Looking at this we can see that b gives us a better approximation of df(0), but the problem is that the point x=-1 doesn't exist as we're only given the curve from 0 onwards - for all we know the function could be discontinuous for x < 0. So short of any complex curve extrapolation we are forced to use the forward difference method which we can see gives us a less good approximation of df(0).

I hope that explains why we're forced to use the less-accurate forwards/backwards method at the ends of the graph. The next thing we need to understand is a) why this error propagates inwards for d2f (i.e. the blue and green cells in the table above) and b) why the d2f curve for x_squared shows such a marked jump at both ends.

The reason this error propagates inward is because the green cell subscribes to the yellow one and the yellow one has the error, the error is propagated. The blue cell has two problems - it also subscribes to the yellow cell but it also suffers from the backwards difference error, which explains why it's error is greater than the green cell.

The reason the d2f error is so marked for the graph x_squared is that the underlying curve (df) for which we're trying to find the derivative is (or should be, analytically) a linear slope. As it's not, due to the 1st order problem described above, the d2f curve which says "there is no rate of change in df" actually does see a rate of change for the first and last elements due to the forwards/backwards error, and for this reason, plus the propagation issue, it shows a marked jump.

Apologies if you think I've laboured this point but it took me a long time to figure out what was happening here and convince myself that these errors weren't due to bugs in my code, and I wanted to save you the trouble.

Step Size

The error across the first derivative is a known quantity of the finite difference method and is proportional to the step-size squared, so the smaller your step size (i.e. the greater NUM_NODES) the better approximation you get. Below, an abridged table showing the sin(5x) data for a 1000-node graph where you can see that the diffs are smaller compared to the 100-node graph above. The full 1000-node outputs for sin(5x) and x_squared are on-line.

Below, the graph for the 1000-node data set. You can see the backwards error still exists but for a much smaller part of the graph.

Usage

The default distribution is set up with a 100-node graph and y = sin(5x). You can switch between sin(5x) and x_squared by toggling the mode var on f:

$ t4db
SQL> update f set mode=1;

And then send the reload signal (there's no need to restart):

$ send_reload.sh

If you want to try a different function, modify the calculate_f function in f.c as appropriate, recompile and restart:

calculus$ make
calculus$ calculus

If you want to change the x-range, again modify the calculate_f as appropriate, recompile and restart.

Further Work

Big graphs and the monitor

The need_refresh function in the trad4 internals explicitly lists each object to which any type subscribes. This has worked well until now, but presents a problem with the monitor and big graphs. As the monitor subscribes to 3*NUM_NODES objects (as there are 3 tiers), big graphs, for example NUM_NODES=10k, mean a massive need_refresh function as 30k objects are listed in this function by the precompiler. This means for a long compile time and a big binary - monitor.o is 3MB for 10k nodes whereas the rest of the object files weigh in at a few KB.

This doesn't actually affect run-time performance as because every object to which the monitor subscribes will have changed, lazy evaluation kicks in so only the first object needs to be tested. It's almost certainly the compiler resolving the implication of the lazy evaluation that takes it so long and makes the binary so big.

Nth-differentials

Should you need the 3rd differential or above it's pretty trivial to add another tier.

Switching between data sets

Due to a limitation in the current trad4 infrastructure, if you want to run against the 10 or 1000 node data sets you'll have to set the value of NUM_NODES in calculus/src/constants.t4s as appropriate, run the precompiler and recompile:

heat_equation_1d$ t4p
heat_equation_1d$ make all