Functional Programming in Numerical Integration

Alexandr Stefek

4.67/5 (4 votes)

Jun 4, 2014

CPOL

5 min read

20943

318

The article depicts usage of functional programming for creating indefinite integral in form Func from delegate Func.

Introduction
Theory
First Part of Code
Antiderivative Method - Other Implementations
Testing of Different Implementations of Antiderivate Method
Usage of Antiderivative Method
Multiple Integration
2D Integrals Numerical Evaluation
3D Integrals Numerical Evaluation
Results
References
History

Introduction

The functional programming is widely used and discussed. My first experience with it is still fresh. During my studies, I read a famous article http://www.codeproject.com/Articles/375166/Functional-programming-in-Csharp which made some "clicks" in my mind.

Because I did some experiments with numerical integration and I found articles http://www.codeproject.com/Articles/95073/Using-Nested-Lambda-Functions-to-Implement-Numeric and http://www.codeproject.com/Articles/334308/Numerical-Integration-with-Simpsons-Rule, I decided to try something similar but in a slightly different way.

Even if the presented functions / methods are fully functional, they are not optimal and probably you will be able to find faster implementations. But still here remains some usable concept.

Because I really like the extension, the whole code is using it. For more explanation on extension methods, you can look at http://www.codeproject.com/Articles/261639/Extension-Methods-in-NET or http://msdn.microsoft.com/en-us/library/bb383977.aspx.

Theory

As mentioned above, the article is focused on numerical integration. Theory about integrals can be easily learned at wiki pages http://en.wikipedia.org/wiki/Integral and http://en.wikipedia.org/wiki/Antiderivative.

The equation (image) was taken from http://en.wikipedia.org/wiki/Integral.

In C# language, we want to develop a function which will be formally declared as:

public static Func<double, double> Antiderivative(this Func<double, double> df)

so the user will be able to call something like:

Func<double, double> cos = Math.Cos;
Func<double, double> sin = cos.Antiderivative();

First Part of Code

The first method which I implemented was Antiderivative function in form which is very close to form introduced in code strip above. For implementation, I used Runge-Kutta method which is the general equivalent of Simpson rule.

/// <summary>
/// limit for numerical computation
/// </summary>
public static double eps = 1e-16;

/// <summary>
/// Numerical antiderivation (integral computation)
/// <para/><see cref="http://en.wikipedia.org/wiki/Antiderivative"/>
/// <para/><see cref="https://www.khanacademy.org/math/integral-calculus/
/// indefinite-definite-integrals/indefinite_integrals/v/antiderivatives-and-indefinite-integrals"/>
/// <para/><see cref="http://www-math.mit.edu/~djk/calculus_beginners/chapter16/section01.html"/>
/// </summary>
/// <example>
/// <code>Func&lt;double, double&gt; sin = Antiderivative(cos, 0, 0, 0.001);</code>
/// </example>
/// <param name="df">derivation</param>
/// <param name="x0">start point</param>
/// <param name="f0">start value f0 = f(x0)</param>
/// <param name="stepSize">step size</param>
/// <returns>integral of input function</returns>
public static Func<double, double> Antiderivative
    (this Func<double, double> df, double x0, double f0, double stepSize)
{
    //save values - state for function which will be returned
    double currentx = x0;
    double currentf = f0;
    double currentdf = df(x0);

    //create function
    return (x) =>
    {
        //not implemented for left values, can be done in future
        if (x < x0)
            throw new ArgumentException("Cannot compute");

        //test relation between currentx (state value) and wanted x | f(x)
        if (currentx > x)
        {
            //reset state if needed
            //this can be improved by values caching this implementation is simplified
            currentx = x0;
            currentf = f0;
            currentdf = df(x0);
        }

        double newLimit = x - stepSize;
        //make loop till the distance is less than deltax
        while (currentx < newLimit)
        {
            //*
            double k1 = stepSize * currentdf; //it is same as k1 = deltax * df(currentx);
            double k2 = stepSize * df(currentx + 0.5 * stepSize);
            double k3 = k2; //it is same as k3 = deltax * df(currentx + 0.5 * deltax);
            currentdf = df(currentx + stepSize);
            double k4 = stepSize * currentdf;
            currentf = currentf + 1 / 6.0 * (k1 + 2 * k2 + 2 * k3 + k4);
            //*/

            currentx = currentx + stepSize;
        }

        //now compute last step, if it is needed
        double lastStepSize = x - currentx;
        if (lastStepSize > eps)
        {
            currentf = currentf + lastStepSize * df(currentx);
            currentx = x;
        }
        //and return value, store it for future use
        return currentf;
    };
}

Antiderivative Method - Other Implementations

Thanks to Marc Clifton (see discussion part bellow article) I included other possible implementations of Antiderivate method.

Firstly the implementation where I avoid use currentx, currentf and currentdf. Full code of AntiderivativePure method:

/// <summary>
/// Numerical antiderivation (integral computation)
/// <para/><see cref="http://en.wikipedia.org/wiki/Antiderivative"/>
/// <para/><see cref="
///https://www.khanacademy.org/math/integral-calculus/indefinite-definite-integrals/
///indefinite_integrals/v/antiderivatives-and-indefinite-integrals"/>
/// <para/><see cref="http://www-math.mit.edu/~djk/calculus_beginners/chapter16/section01.html"/>
/// </summary>
/// <example>
/// <code>Func&lt;double, double&gt; sin = Antiderivative(cos, 0, 0, 0.001);</code>
/// </example>
/// <param name="df">derivation</param>
/// <param name="x0">start point</param>
/// <param name="f0">start value f0 = f(x0)</param>
/// <param name="stepSize">step size</param>
/// <returns>integral of input function</returns>
public static Func<double, double> AntiderivativePure
    (this Func<double, double> df, double x0, double f0, double stepSize)
{
    //create function
    return (x) =>
    {
        //not implemented for left values, can be done in future
        if (x < x0)
            throw new ArgumentException("Cannot compute");

        double currentx = x0;
        double currentf = f0;
        double currentdf = df(x0);

        double newLimit = x - stepSize;
        //make loop till the distance is less than deltax
        while (currentx < newLimit)
        {
            //*
            double k1 = stepSize * currentdf; //it is same as k1 = deltax * df(currentx);
            double k2 = stepSize * df(currentx + 0.5 * stepSize);
            double k3 = k2; //it is same as k3 = deltax * df(currentx + 0.5 * deltax);
            currentdf = df(currentx + stepSize);
            double k4 = stepSize * currentdf;
            currentf = currentf + 1 / 6.0 * (k1 + 2 * k2 + 2 * k3 + k4);
            //*/

            currentx = currentx + stepSize;
        }

        //now compute last step, if it is needed
        double lastStepSize = x - currentx;
        if (lastStepSize > eps)
            currentf = currentf + lastStepSize * df(currentx);

        //and return value, store it for future use
        return currentf;
    };
}

As an second implementation the recursive implementation will be shown. Notice the expression

return df.AntiderivativeRecursive(x0, f0, stepSize)(x - stepSize) + delta;

Full code of AntiderivativeRecursive method:

/// <summary>
/// Numerical antiderivation (integral computation)
/// <para/><see cref="http://en.wikipedia.org/wiki/Antiderivative"/>
/// <para/><see cref="https://www.khanacademy.org/math/integral-calculus/
/// indefinite-definite-integrals/indefinite_integrals/v/antiderivatives-and-indefinite-integrals"/>
/// <para/><see cref="http://www-math.mit.edu/~djk/calculus_beginners/chapter16/section01.html"/>
/// </summary>
/// <example>
/// <code>Func&lt;double, double&gt; sin = Antiderivative(cos, 0, 0, 0.001);</code>
/// </example>
/// <param name="df">derivation</param>
/// <param name="x0">start point</param>
/// <param name="f0">start value f0 = f(x0)</param>
/// <param name="stepSize">step size</param>
/// <returns>integral of input function</returns>
public static Func<double, double> AntiderivativeRecursive
    (this Func<double, double> df, double x0, double f0, double stepSize)
{
    //create function
    return (x) =>
    {
        //not implemented for left values, can be done in future
        if (x < x0)
            throw new ArgumentException("Cannot compute");

        if (x - x0 < stepSize)
        {
            double cstep = x - x0;
            double k1 = cstep * df(x0);
            double k2 = cstep * df(x0 + 0.5 * cstep);
            double k3 = k2; //it is same as k3 = cstep * df(x0 + 0.5 * cstep);
            double k4 = cstep * df(x0 + cstep);
            double currentf = f0 + 1 / 6.0 * (k1 + 2 * k2 + 2 * k3 + k4);

            return currentf;
        }
        else
        {
            double k1 = stepSize * df(x - stepSize);
            double k2 = stepSize * df(x - 0.5 * stepSize);
            double k3 = k2; //it is same as k3 = stepSize * df(x - 0.5 * stepSize);
            double k4 = stepSize * df(x);
            double delta = 1 / 6.0 * (k1 + 2 * k2 + 2 * k3 + k4);

            return df.AntiderivativeRecursive(x0, f0, stepSize)(x - stepSize) + delta;
        }
    };
}

And as last implementation (at least at this stage) AntiderivativeEx method:

/// <summary>
/// Numerical antiderivation (integral computation)
/// <para/><see cref="http://en.wikipedia.org/wiki/Antiderivative"/>
/// <para/><see cref="https://www.khanacademy.org/math/
/// integral-calculus/indefinite-definite-integrals/indefinite_integrals/v/antiderivatives-and-indefinite-integrals"/>
/// <para/><see cref="http://www-math.mit.edu/~djk/calculus_beginners/chapter16/section01.html"/>
/// </summary>
/// <example>
/// <code>Func&lt;double, double&gt; sin = Antiderivative(cos, 0, 0, 0.001);</code>
/// </example>
/// <param name="df">derivation</param>
/// <param name="x0">start point</param>
/// <param name="f0">start value f0 = f(x0)</param>
/// <param name="stepSize">step size</param>
/// <returns>integral of input function</returns>
public static Func<double, double> AntiderivativeEx
    (this Func<double, double> df, double x0, double f0, double stepSize)
{
    return (x) => df.AntiderivativeExH(x0, f0, stepSize)(x, df(x));
}

private static Func<double, double, double> AntiderivativeExH
    (this Func<double, double> df, double x0, double f0, double stepSize)
{
    //create function
    return (x, dfx) =>
    {
        //not implemented for left values, can be done in future
        if (x < x0)
            throw new ArgumentException("Cannot compute");

        if (x - x0 < stepSize)
        {
            double cstep = x - x0;
            double k1 = cstep * df(x0);
            double k2 = cstep * df(x0 + 0.5 * cstep);
            double k3 = k2; //it is same as k3 = cstep * df(x0 + 0.5 * cstep);
            double k4 = cstep * dfx; //it is same as k4 = cstep * df(x0 + cstep);
            double currentf = f0 + 1 / 6.0 * (k1 + 2 * k2 + 2 * k3 + k4);

            return currentf;
        }
        else
        {
            double dfxmh = df(x - stepSize);
            double k1 = stepSize * dfxmh;
            double k2 = stepSize * df(x - 0.5 * stepSize);
            double k3 = k2; //it is same as k3 = stepSize * df(x - 0.5 * stepSize);
            double k4 = stepSize * dfx;
            double delta = 1 / 6.0 * (k1 + 2 * k2 + 2 * k3 + k4);

            return df.AntiderivativeExH(x0, f0, stepSize)(x - stepSize, dfxmh) + delta;
        }
    };
}

This last implementation uses function Func<double, double, double> which can accommodate already calculated value of dfxm = df(x - stepSize). Once the value is computed, it is shared with other evaluations (see last expression df.AntiderivativeExH(x0, f0, stepSize)(x - stepSize, dfxmh) in code above).

Testing of Different Implementations of Antiderivate Method

Of course there are differences in implementation and a good question is how big an impact on performance those differences will have? To be able to answer this question, I prepared a test where two aspects - call count and time consumption - will be observed.

The test code is as follows:

private void AntiderivativePerformanceTest()
{
    int calls = 0;
    Func<double, double> cos = (x) =>
    {
        calls++;
        return Math.Cos(x);
    };

    Func<double, double> sin = cos.Antiderivative(0, 0, 0.01);
    Func<double, double> sinB = cos.AntiderivativePure(0, 0, 0.01);
    Func<double, double> sinC = cos.AntiderivativeRecursive(0, 0, 0.01);
    Func<double, double> sinD = cos.AntiderivativeEx(0, 0, 0.01);

    int cycles = 1000;
    calls = 0;
    int start = Environment.TickCount;
    for (int i = 0; i < cycles; i++)
        for (double x = 0; x < Math.PI; x = x + 0.1)
        {
            double sinValueI = sin(x);
        }
    int elapsed = Environment.TickCount - start;
    Msg(string.Format("Antiderivative: {0}ms, calls: {1}", elapsed, calls));

    calls = 0;
    start = Environment.TickCount;
    for (int i = 0; i < cycles; i++)
        for (double x = 0; x < Math.PI; x = x + 0.1)
        {
            double sinValueI = sinB(x);
        }

    elapsed = Environment.TickCount - start;
    Msg(string.Format("AntiderivativePure: {0}ms, calls: {1}", elapsed, calls));

    calls = 0;
    start = Environment.TickCount;
    for (int i = 0; i < cycles; i++)
        for (double x = 0; x < Math.PI; x = x + 0.1)
        {
            double sinValueI = sinC(x);
        }
    elapsed = Environment.TickCount - start;
    Msg(string.Format("AntiderivativeRecursive: {0}ms, calls: {1}", elapsed, calls));

    calls = 0;
    start = Environment.TickCount;
    for (int i = 0; i < cycles; i++)
        for (double x = 0; x < Math.PI; x = x + 0.1)
        {
            double sinValueI = sinD(x);
        }
    elapsed = Environment.TickCount - start;
    Msg(string.Format("AntiderivativeEx: {0}ms, calls: {1}", elapsed, calls));
}

All implementations are used for evaluation of 3200 values of integral function. Results of this tests are given below:

Antiderivative: 46ms, calls: 642999
AntiderivativePure: 750ms, calls: 9974000
AntiderivativeRecursive: 1625ms, calls: 14919000
AntiderivativeEx: 1235ms, calls: 9978000

As we can read, Antiderivative method is fastes (46ms) and does least calls of df function (643k). If we avoid recursive implementation and use of state variables (currentx, currentf and currentdf), we get 750ms and nearly one million of calls (0.997M).

There are two implementations with use of recursive calls AntiderivativeRecursive and AntiderivativeEx. AntiderivativeEx is faster (1235ms) and does less call (0.998M) than AntiderivativeRecursive. Both aspects are better because of optimization of df function calls. This optimization is done via sharing of derivation values between recursive calls as you can see in appropriate code (AntiderivativeExH method).

Usage of Antiderivative Method

For a test of Antiderivative method, I choose test with sin(x) and cos(x).

private void AntiderivativeTest()
{
    Func<double, double> cos = Math.Cos;
    Func<double, double> sin = cos.Antiderivative(0, 0, 0.01);
    for (double x = 0; x < Math.PI; x = x + 0.1)
    {
        double sinValueI = sin(x);
        double sinValue = Math.Sin(x);
        double delta = sinValue - sinValueI;
        Msg(string.Format("x = {0:F6}, Math.Sin(x) = {1:F6}, sin(x) = {2:F6}, delta = {3:F6}", x, sinValue, sinValueI, delta));
    }
}

The result of method is:

x = 0.000000, Math.Sin(x) = 0.000000, sin(x) = 0.000000, delta = 0.000000
x = 0.100000, Math.Sin(x) = 0.099833, sin(x) = 0.099833, delta = 0.000000
x = 0.200000, Math.Sin(x) = 0.198669, sin(x) = 0.198669, delta = 0.000000
x = 0.300000, Math.Sin(x) = 0.295520, sin(x) = 0.295520, delta = 0.000000
x = 0.400000, Math.Sin(x) = 0.389418, sin(x) = 0.389418, delta = 0.000000
x = 0.500000, Math.Sin(x) = 0.479426, sin(x) = 0.479426, delta = 0.000000
x = 0.600000, Math.Sin(x) = 0.564642, sin(x) = 0.564642, delta = 0.000000
x = 0.700000, Math.Sin(x) = 0.644218, sin(x) = 0.644218, delta = 0.000000
x = 0.800000, Math.Sin(x) = 0.717356, sin(x) = 0.717356, delta = 0.000000
x = 0.900000, Math.Sin(x) = 0.783327, sin(x) = 0.783327, delta = 0.000000
x = 1.000000, Math.Sin(x) = 0.841471, sin(x) = 0.841471, delta = 0.000000
x = 1.100000, Math.Sin(x) = 0.891207, sin(x) = 0.891207, delta = 0.000000
x = 1.200000, Math.Sin(x) = 0.932039, sin(x) = 0.932039, delta = 0.000000
x = 1.300000, Math.Sin(x) = 0.963558, sin(x) = 0.963558, delta = 0.000000
x = 1.400000, Math.Sin(x) = 0.985450, sin(x) = 0.985450, delta = 0.000000
x = 1.500000, Math.Sin(x) = 0.997495, sin(x) = 0.997495, delta = 0.000000
x = 1.600000, Math.Sin(x) = 0.999574, sin(x) = 0.999574, delta = 0.000000
x = 1.700000, Math.Sin(x) = 0.991665, sin(x) = 0.991665, delta = 0.000000
x = 1.800000, Math.Sin(x) = 0.973848, sin(x) = 0.973848, delta = 0.000000
x = 1.900000, Math.Sin(x) = 0.946300, sin(x) = 0.946300, delta = 0.000000
x = 2.000000, Math.Sin(x) = 0.909297, sin(x) = 0.909297, delta = 0.000000
x = 2.100000, Math.Sin(x) = 0.863209, sin(x) = 0.863209, delta = 0.000000
x = 2.200000, Math.Sin(x) = 0.808496, sin(x) = 0.808496, delta = 0.000000
x = 2.300000, Math.Sin(x) = 0.745705, sin(x) = 0.745705, delta = 0.000000
x = 2.400000, Math.Sin(x) = 0.675463, sin(x) = 0.675463, delta = 0.000000
x = 2.500000, Math.Sin(x) = 0.598472, sin(x) = 0.598472, delta = 0.000000
x = 2.600000, Math.Sin(x) = 0.515501, sin(x) = 0.515501, delta = 0.000000
x = 2.700000, Math.Sin(x) = 0.427380, sin(x) = 0.427380, delta = 0.000000
x = 2.800000, Math.Sin(x) = 0.334988, sin(x) = 0.334988, delta = 0.000000
x = 2.900000, Math.Sin(x) = 0.239249, sin(x) = 0.239249, delta = 0.000000
x = 3.000000, Math.Sin(x) = 0.141120, sin(x) = 0.141120, delta = 0.000000
x = 3.100000, Math.Sin(x) = 0.041581, sin(x) = 0.041581, delta = 0.000000

Once we get a method for creating antiderivate function, we can use it for numerical solution of definite integrals:

The equation (image) was taken from http://en.wikipedia.org/wiki/Integral.

So the method in C# can look like:

/// <summary>
/// Takes function and evals f(to) - f(from)
/// </summary>
/// <param name="f"></param>
/// <param name="from"></param>
/// <param name="to"></param>
/// <returns></returns>
public static double IntervalEval(this Func<double, double> f, double from, double to)
{
    double l = f(from);
    double h = f(to);
    return h - l;
    //abit better than
    //return f(to) - f(from);
    //especially in usage with Antiderivative function
}

/// <summary>
/// Calculates integral of fx over interval where a is low bound and b is high bound
/// <see cref="http://en.wikipedia.org/wiki/Integral"/>
/// </summary>
/// <param name="fx">function to integrate</param>
/// <param name="a">low bound</param>
/// <param name="b">high bound</param>
/// <param name="stepSize">step size</param>
/// <returns>Integral value</returns>
public static double Integral(this Func<double, double> fx, double a, double b, double stepSize)
{
    return fx.Antiderivative(a, 0, stepSize).IntervalEval(a, b);
}

Interval method is here for simplification of interval evaluation of equation right side. Notice that body of Integral method consists of just one line.

At this stage of implementation, the question "What about 2D integrals?" will appear.

Multiple Integration

Multiple integration is discussed on Wikipedia, see http://en.wikipedia.org/wiki/Integral#Extensions - Multiple integration or http://en.wikipedia.org/wiki/Multiple_integral. A very common form of multiple integration is 2D and 3D integrals which can be used for area or volume calculations.

2D Integrals Numerical Evaluation

2D integrals can be computed by usage of the next equation:

The equation (image) was taken from http://en.wikipedia.org/wiki/Multiple_integral.

/// <summary>
/// Create function f(y) from f(x, y) in form f(const x, y)
/// </summary>
/// <param name="fxy">function f(x, y)</param>
/// <param name="x">const x</param>
/// <returns>function f(y)</returns>
public static Func<double, double> FixateX(this Func<double, double, double> fxy, double x)
{
    return (y) => fxy(x, y);
}

/// <summary>
/// Calculates 2D integral of f over shape defined by lowx, highx, lowy(x) and highy(x) bounds.
/// </summary>
/// <param name="fxy">function to integrate</param>
/// <param name="lowX">x low bound</param>
/// <param name="highX">x high bound</param>
/// <param name="lowY">y(x) func low bound</param>
/// <param name="highY">y(x) func high bound</param>
/// <param name="stepSize">step size</param>
/// <returns>Integral value</returns>
public static double Integral(this Func<double, double, double> fxy, 
double lowX, double highX, Func<double, double> lowY, Func<double, double> highY, double stepSize)
{
    Func<double, double> helperFunc = 
        (cx) => fxy.FixateX(cx).Integral(lowY(cx), highY(cx), stepSize);
    return helperFunc.Integral(lowX, highX, stepSize);
}

FixateX method is here for simplification. The implementation can be shortened if this general method is omitted. Even if the body of Integral method can be typed on one line, the two expression implementation is much easier to understand (maybe one more line can be optimal).

Test for this method can look like:

//calculate circle with radius r area (Pi*r*r)
private double Circle(double r)
{
    Func<double, double, double> fxy = (x, y) => 1;
    Func<double, double> yboundLow = (x) => -Math.Sqrt(r * r - x * x);
    Func<double, double> yboundHigh = (x) => Math.Sqrt(r * r - x * x);

    double result = 0;
    result = fxy.Integral(-r, r, yboundLow, yboundHigh, 0.01);
    return result;
}

private void I2DTest()
{
    Msg("Circle = " + Circle(1));
    Msg("Expected : " + Math.PI);
}

Result of I2DTest is:

Circle = 3.14143024919302
Expected : 3.14159265358979

What now 3D integration? Of course! Let's do it!

3D Integrals Numerical Evaluation

For 3D integrals, we can use a very similar equation (compare with 2D):

The equation was taken from http://en.wikipedia.org/wiki/Multiple_integral.

C# implementation can be this:

/// <summary>
/// Create function f(y, z) from f(x, y, z) in form f(const x, y, z)
/// </summary>
/// <param name="fxy">function f(x, y, z)</param>
/// <param name="x">const x</param>
/// <returns>function f(y, z)</returns>
public static Func<double, double, double> FixateX
    (this Func<double, double, double, double> fxyz, double x)
{
    return (y, z) => fxyz(x, y, z);
}

/// <summary>
/// Calculates 3D integral of f over volume defined 
/// by lowx, highx, lowy(x), highy(x), lowz(x, y) and highz(x, y) bounds.
/// </summary>
/// <param name="fxyz">function to integrate</param>
/// <param name="lowX">x low bound</param>
/// <param name="highX">x high bound</param>
/// <param name="lowY">y(x) func low bound</param>
/// <param name="highY">y(x) func high bound</param>
/// <param name="lowZ">z(x, y) func low bound</param>
/// <param name="highZ">z(x, y) func high bound</param>
/// <param name="stepSize">step size</param>
/// <returns>Integral value</returns>
public static double Integral(this Func<double, double, double, double> fxyz,
    double lowX, double highX, Func<double, double> lowY, Func<double, 
    double> highY, Func<double, double, double> lowZ, Func<double, 
    double, double> highZ, double stepSize)
{
    Func<double, double> helperFunc = 
        (cx) => fxyz.FixateX(cx).Integral(lowY(cx), 
        highY(cx), lowZ.FixateX(cx), highZ.FixateX(cx), stepSize);
    return helperFunc.Integral(lowX, highX, stepSize);
}

FixateX method if here for simplification. The implementation can be shorted if this general method is omitted. Even if the body of Integral method can be typed on one line, the two expression implementation is much easier to understand (maybe one more line can be optimal). Really the same words like in 2D...