using System;
using System.Collections.Generic;
using System.Linq;
using System.Text;
namespace FunctionalProgramming
{
//generally inspired by
//http://www.codeproject.com/Articles/375166/Functional-programming-in-Csharp
//and reaction to
//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
public static class Extensions
{
/// <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<double, double> 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)
{
//*
double k1 = lastStepSize * currentdf; //it is same as k1 = deltax * df(currentx);
double k2 = lastStepSize * df(currentx + 0.5 * lastStepSize);
double k3 = k2; //it is same as k3 = deltax * df(currentx + 0.5 * deltax);
currentdf = df(currentx + lastStepSize);
double k4 = lastStepSize * currentdf;
currentf = currentf + 1 / 6.0 * (k1 + 2 * k2 + 2 * k3 + k4);
//*/
//currentf = currentf + lastStepSize * df(currentx);
currentx = x;
}
//and return value, store it for future use
return currentf;
};
}
/// <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<double, double> 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)
{
//*
double k1 = lastStepSize * currentdf; //it is same as k1 = deltax * df(currentx);
double k2 = lastStepSize * df(currentx + 0.5 * lastStepSize);
double k3 = k2; //it is same as k3 = deltax * df(currentx + 0.5 * deltax);
currentdf = df(currentx + lastStepSize);
double k4 = lastStepSize * currentdf;
currentf = currentf + 1 / 6.0 * (k1 + 2 * k2 + 2 * k3 + k4);
//*/
//currentf = currentf + lastStepSize * df(currentx);
currentx = x;
}
//and return value, store it for future use
return currentf;
};
}
/// <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<double, double> 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;
}
};
}
/// <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<double, double> 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;
}
};
}
/// <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 f over interval where a is low bound and b is high bound
/// <see cref="http://en.wikipedia.org/wiki/Integral"/>
/// </summary>
/// <param name="f">function to integrate</param>
/// <param name="a">low bound</param>
/// <param name="b">high bound</param>
/// <param name="deltax">step size</param>
/// <returns>Integral value</returns>
public static double Integral(this Func<double, double> f, double a, double b, double deltax)
{
return f.Antiderivative(a, 0, deltax).IntervalEval(a, b);
}
/// <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);
}
/// <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);
}
/// <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>
/// 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);
}
}
}
Submit an article or tip