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)
