using System;
using System.Linq;
namespace Computationals
{
// Delegates
delegate double Function(double x);
delegate double Functionxy(double x, double y);
delegate double Integral(Function f, double a, double b, int n);
// for integration of function f over n subintervals from a to b
delegate double Differentiation(Function f, double a);
// for differentiation and double differentiation of function f at x=a in the domain of f
//___________________________________________________________________________________________________//
public class LambdaNumericals
{
// simple test functions
public static double p(double x, double y)
{
return x * y;
}
public static double q(double x, double y)
{
return x * x * y * y * y;
}
public static double r(double x, double y)
{
return 5;
}
public static double hemiSphere(double x, double y) // hemisphere radius 3
{
return Math.Sqrt(9 - x * x - y * y);
}
public static double one(double x)
{
return 1;
}
public static double minusOne(double x)
{
return -1;
}
public static double f1(double x)
{
return x;
}
public static double f2(double x)
{
return 2 * x;
}
public static double sine(double x)
{
return Math.Sin(x);
}
public static double semiCirc(double x)
{
return Math.Sqrt(1 - x * x);
}
//___________________________________________________________________________________________________//
// Differentiation methods - extend to include more 2 and 3 point recipes with foward & backward interpolations & variable step
static Differentiation D2pointCen()
{
double h = 0.000001; // h could vary in an adaptive algorithm
Differentiation differentiation = (f, x) =>
(f(x + h) - f(x - h)) / (2 * h);
return differentiation;
}
static Differentiation D5pointCen()
{
double h = 0.0005;
Differentiation differentiation = (f, x) =>
(-f(x + 2 * h) + 8 * f(x + h) - 8 * f(x - h) + f(x - 2 * h)) / (12 * h);
return differentiation;
}
static Differentiation D5pointFwd()
{
double h = 0.0005;
Differentiation differentiation = (f, x) =>
(-3 * f(x + 4 * h) + 16 * f(x + 3 * h) - 36 * f(x + 2 * h) + 48 * f(x + h) - 25 * f(x)) / (12 * h);
return differentiation;
}
// a double differentiation method -
static Differentiation DD5pointCen()
{
double h = 0.002;
Differentiation ddifferentiation = (f, x) =>
(-f(x + 2 * h) + 16 * f(x + h) - 30 * f(x) + 16 * f(x - h) - f(x - 2 * h)) / (12 * h * h);
return ddifferentiation;
}
//___________________________________________________________________________________________________//
// Integral methods - Romberg and adaptive methods could be added
static Integral SimpInt()
{
// Simpson's rule integration of function f, over n (even) subintervals in the interval [a,b]
Integral integral = (f, a, b, n) =>
{
// n, not even, throw exception code should go here
double sum = 0;
double h = (b - a) / n;
for (int i = 0; i < n; i = i + 2)
sum = sum + (f(a + i * h) + 4 * f(a + (i + 1) * h) + f(a + (i + 2) * h)) * h / 3;
return sum;
};
return integral;
}
static Integral GaussInt()
{
// 4-point Gaussian rule integration of the function f, over n subintervals in the interval (a,b)
Integral integral = (f, a, b, n) =>
{
double[] weight = new double[4];
weight[2] = 0.652145154862546;
weight[1] = weight[2];
weight[3] = 0.347854845137454;
weight[0] = weight[3];
double[] point = new double[4];
point[3] = 0.861136311594053;
point[0] = -point[3];
point[2] = 0.339981043584856;
point[1] = -point[2];
double sum = 0;
double h = (b - a) / n;
double xr, xl, m, c, s;
for (int i = 0; i < n; i++)
{
s = 0;
xl = a + i * h;
xr = xl + h;
m = (xr - xl) / 2;
c = (xr + xl) / 2;
for (int j = 0; j < 4; j++)
s = s + weight[j] * f(m * point[j] + c) * m;
sum = sum + s;
}
return sum;
};
return integral;
}
//___________________________________________________________________________________________________//
// Double integration using nested lambda functions
// method below could be extended to 3D or higher integration
// 2D integration of a function of 2 variables. (( f(x,y)dxdy
// ))
// Inner integral dx, from x=g(y) to x=h(y), over m subintervals, using algorithm inInt
// Outer integral dy, from y=a to y=b, over n subintervals, using algorithm outInt
static double dubInt(Functionxy fxy, Function g, Function h, double a, double b, int m, int n, Integral inInt, Integral outInt)
{
return outInt((double y) => inInt((double x) => fxy(x, y), g(y), h(y), m), a, b, n);
}
//___________________________________________________________________________________________________//
// partial derivatives using nested lambda functions
static double PartialDiffx(Functionxy fxy, double a, double b)
{
Differentiation diff = D2pointCen(); // could be parameter
return diff((double x) => fxy(x, b), a);
}
static double PartialDiffy(Functionxy fxy, double a, double b)
{
Differentiation diff = D2pointCen(); // could be parameter
return diff((double y) => fxy(a, y), b);
}
//___________________________________________________________________________________________________//
// Length of curve f(x) from x=a to x=b
// Uses a nested lambda function
static double curveLength(Function f, double a, double b, Integral intMethod, int n, Differentiation diffMethod)
{
return intMethod((double x) => Math.Sqrt(1 + Math.Pow(diffMethod(f, x), 2)), a, b, n);
}
//___________________________________________________________________________________________________//
// Area under z = fxy(x,y) from y=f(x) to y=g(x) and x=a to x=b
// Uses nested lambda function
static double surfaceArea(Functionxy fxy, Function f, Function g, double a, double b, int m, int n)
{
Integral intMethod = GaussInt();
return dubInt((double x, double y) =>
Math.Sqrt(1 + Math.Pow(PartialDiffx(fxy, x, y), 2) + Math.Pow(PartialDiffy(fxy, x, y), 2)), f, g, a, b, m, n, intMethod, intMethod);
}
//___________________________________________________________________________________________________//
// Curvature
static double curvature(Function f, double a)
{
Differentiation ddiff = DD5pointCen(); // could be parameters
Differentiation diff = D5pointCen();
return ddiff(f, a) / Math.Pow(1 + Math.Pow(diff(f, a), 2), 1.5);
}
//___________________________________________________________________________________________________//
public static void Main()
{
Integral integrate1, integrate2;
integrate1 = SimpInt();
Console.WriteLine(integrate1(sine, 0, Math.PI, 60));
// result 2.00000008353986
integrate2 = GaussInt();
Console.WriteLine(integrate2(sine, 0, Math.PI, 15));
// result 2
Console.WriteLine(integrate2(Math.Sin, 0, Math.PI, 15));
// result 2
// analytically, exactly 2 for the above integrals - note the superior convergence of Gauss 4-point
Console.WriteLine(dubInt(p, f1, f2, 1, 2, 4, 4, integrate1, integrate2));
// result 5.625
// analytically 5 5/8
Console.WriteLine(D5pointCen()(Math.Sin, Math.PI));
Console.WriteLine(D5pointFwd()(Math.Sin, Math.PI));
Console.WriteLine(DD5pointCen()(Math.Cos, Math.PI));
Console.WriteLine(curveLength(f1, 0, Math.Sqrt(2), GaussInt(), 2, D5pointCen()));
Console.WriteLine(2 * curveLength(semiCirc, -Math.Sqrt(2) / 2, Math.Sqrt(2) / 2, GaussInt(), 10, D5pointCen()));
Console.WriteLine(PartialDiffx(q, 5, 3));
Console.WriteLine(PartialDiffy(q, 5, 3));
Console.WriteLine(surfaceArea(p, f1, f2, 1, 4, 4, 4));
Console.WriteLine(surfaceArea(r, f1, f2, 1, 4, 4, 4));
// 7.5 as expected
Console.WriteLine(surfaceArea((double x,double y) => Math.Sqrt(9 - x * x - y * y), (double x) => -1, (double x) => 1, -1, 1, 20, 20));
// curved surface area of hemisphere, radius 3, above 2 X 2 square; area = 4.16...
Console.WriteLine(curvature((double x) => Math.Sqrt(4 - x * x), 1));
Console.WriteLine(curvature((double x) => (x + 2) * x * (x - 2), -1));
Console.WriteLine(curvature((double x) => (x + 2) * x * (x - 2), 1));
Console.WriteLine(curvature((double x) => (x + 2) * x * (x - 2), 0));
// 2 E -15 (point of inflexion)
Console.WriteLine(curvature((double x) => x * Math.Sin(x), Math.PI * 3 / 2));
Console.Read();
}
}
}
Submit an article or tip