## Contents

- 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.

public static double eps = 1e-16;
public static Func<double, double> Antiderivative
(this Func<double, double> df, double x0, double f0, double stepSize)
{
double currentx = x0;
double currentf = f0;
double currentdf = df(x0);
return (x) =>
{
if (x < x0)
throw new ArgumentException("Cannot compute");
if (currentx > x)
{
currentx = x0;
currentf = f0;
currentdf = df(x0);
}
double newLimit = x - stepSize;
while (currentx < newLimit)
{
double k1 = stepSize * currentdf;
double k2 = stepSize * df(currentx + 0.5 * stepSize);
double k3 = k2;
currentdf = df(currentx + stepSize);
double k4 = stepSize * currentdf;
currentf = currentf + 1 / 6.0 * (k1 + 2 * k2 + 2 * k3 + k4);
currentx = currentx + stepSize;
}
double lastStepSize = x - currentx;
if (lastStepSize > eps)
{
currentf = currentf + lastStepSize * df(currentx);
currentx = x;
}
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:

public static Func<double, double> AntiderivativePure
(this Func<double, double> df, double x0, double f0, double stepSize)
{
return (x) =>
{
if (x < x0)
throw new ArgumentException("Cannot compute");
double currentx = x0;
double currentf = f0;
double currentdf = df(x0);
double newLimit = x - stepSize;
while (currentx < newLimit)
{
double k1 = stepSize * currentdf;
double k2 = stepSize * df(currentx + 0.5 * stepSize);
double k3 = k2;
currentdf = df(currentx + stepSize);
double k4 = stepSize * currentdf;
currentf = currentf + 1 / 6.0 * (k1 + 2 * k2 + 2 * k3 + k4);
currentx = currentx + stepSize;
}
double lastStepSize = x - currentx;
if (lastStepSize > eps)
currentf = currentf + lastStepSize * df(currentx);
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:

public static Func<double, double> AntiderivativeRecursive
(this Func<double, double> df, double x0, double f0, double stepSize)
{
return (x) =>
{
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;
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;
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:

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)
{
return (x, dfx) =>
{
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;
double k4 = cstep * dfx;
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;
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:

public static double IntervalEval(this Func<double, double> f, double from, double to)
{
double l = f(from);
double h = f(to);
return h - l;
}
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.

public static Func<double, double> FixateX(this Func<double, double, double> fxy, double x)
{
return (y) => fxy(x, y);
}
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:

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:

public static Func<double, double, double> FixateX
(this Func<double, double, double, double> fxyz, double x)
{
return (y, z) => fxyz(x, y, z);
}
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...

Test for this method can look like:

private double Sphere(double r)
{
Func<double, double, double, double> fxyz = (x, y, z) => 1;
Func<double, double> yboundLow = (x) => -Math.Sqrt(r * r - x * x);
Func<double, double> yboundHigh = (x) => Math.Sqrt(r * r - x * x);
Func<double, double, double> zboundLow = (x, y) => -Math.Sqrt(r * r - x * x - y * y);
Func<double, double, double> zboundHigh = (x, y) => Math.Sqrt(r * r - x * x - y * y);
double result = 0;
result = fxyz.Integral(-r, r, yboundLow, yboundHigh, zboundLow, zboundHigh, 0.01);
return result;
}
private void I3DTest()
{
Msg("Sphere = " + Sphere(1).ToString());
Msg("Expected : " + 4 * Math.PI / 3);
}

Result of `I3DTest`

is:

Sphere = 4.18857816178346
Expected : 4.18879020478639

## Results

In the Introduction part, there was a warning about optimality of implemented algorithms. Yes, this implementation is not optimal in case of evaluation speed. But why not sometime create something "smooth". This work helped me to become more familiar with functional programming.

## References

## History

- 04-Jun-2014 First implementation (from scratch)
- 11-Jun-2014 Update according discussion (and fixed bug - improved accuracy)