Real time (stream) integration with Simpson's rule

Turn on MathML plugin

About this method of numerical integration you can read here. To do it in streaming manner we should do little formal mathematics:

Composite Simpson's rule:

$\underset{a}{\overset{b}{\int }}f\left(x\right)dx\approx \frac{h}{3}\sum _{k=1,2}^{N-1}\left(f\left(x_{k-1}\right)+4f\left(x_k\right)+f\left(x_{k+1}\right)\right)$

where k=1,2 means from 1 step by 2, N is even intervals number.

The sum will be: $\frac{h}{3}\sum _{k=1,2}^{N-1}{S}_k$

Simpsone's method requires 3 points least so we can integrate not every point (sample) but through one, so neighboring samples are I_N-2 and I_N. Different between them is:

$I_N-I_{N-2}=\frac{h}{3}\left(\sum _{k=1,2}^{N-1}{S}_k-\sum _{k=1,2}^{N-3}{S}_k\right)=\frac{h}{3}\left(S_{N-1}+\sum _{k=1,2}^{N-3}{S}_k-\sum _{k=1,2}^{N-3}{S}_k\right)$

no S_N-2 member bcz of step k=2.

So difference is: $\frac{h}{3}{S}_{N-1}$

For example, $I_4-I_2=\frac{h}{3}\left(y_2+4y_3+y_4\right)$, where I₄ is 4^th interval but 5^th sample.

And resulting formula of I_i=

$0,i\le 1$
$\frac{h}{3}\left(y_{i-2}+4y_{i-1}+y_i\right)+I_{i-1},imod2=0$
$I_{i-1},imod2\ne 0$

where i is a number of input sample. We can not get values of I on each sample, but through one (skipping odds and 0-th, 1-th). On odds samples, value will repeat previous one.

Example in SWIG:
some.h file:

// Simpsone's rule integrator
typedef struct CSIntr {
        double I; // last value of integral
        double y_1; // y[-1]
        double y_2; // y[-2]
        double h; // step of subinterval
        double h_d; // step of subinterval/3
        long long i; // iteration
} CSIntr;

some.i file:

%include "some.h"
%extend CSIntr {
    // FIXME how is more right to do this?
    void reset() {
        $self->I = 0.;
        $self->y_1 = $self->y_2 = 0.;
        $self->i = 0;
    }

    void setup(double h) {
        $self->h = h;
        $self->h_d = h/3.;
    }

    CSIntr(double h) {
        CSIntr *obj;

        if (NULL == (obj=(CSIntr*) malloc(sizeof (CSIntr))))
            return (NULL);
        CSIntr_reset(obj);
        CSIntr_setup(obj, h);
        return (obj);
    }

    ~CSIntr() {
        if ($self) free($self);
    }

    double calculate(double y) {
        if ($self->i <= 1)
            $self->I = 0;
        else if (0 == $self->i % 2)
            $self->I += $self->h_d * ($self->y_2 + 4. * $self->y_1 + y);

        $self->y_2 = $self->y_1;
        $self->y_1 = y;
        $self->i++;
        return ($self->I);
    }
};

So, if we use Tcl, then we can build and test our class:
Makefile's fragment:

some_wrap.c: some.i some.h some.c
 swig -tcl -pkgversion 1.0 some.i

some.dll: some_wrap.c
 gcc -Ic:/tcl/include -I./ -DUSE_TCL_STUBS -c some_wrap.c -o some_wrap.o
 gcc -shared -o some.dll some/some_wrap.o -L c:/tcl/lib -ltclstub85 -lm

test on linear function with y_i = {1, 2, 3, 4, 5}:

load some
% CSIntr cs 1
_f85f7601_p_CSIntr
% cs calculate 1
0.0
% cs calculate 2
0.0
% cs calculate 3
4.0
% cs calculate 4
4.0
% cs calculate 5
12.0

I₄ is really 12.0