Skip to content

Instantly share code, notes, and snippets.

Embed
What would you like to do?
Some horribly unoptimized quadrature rules, integrating the "black swan" function from Eric Lippert's "Fixing Random" series.
using System;
using System.Collections.Generic;
using System.Linq;
namespace Quadratures
{
interface Integrator
{
double Integrate(Func<double, double> f, double a, double b);
}
public static class Functions
{
public static double WeirdFunction(double x)
{
return Math.Atan(1000 * (x - 0.45)) * 20 - 31.2;
}
static readonly double INVERSE_ROOT_PI = 1.0 / Math.Sqrt(2 * Math.PI);
public static double NormalDistributionDensity(double x, double mu, double sigma)
{
return Math.Exp(-(x - mu) * (x - mu) / (2 * sigma * sigma)) * INVERSE_ROOT_PI / sigma;
}
}
struct PointWeight
{
public double Point { get; }
public double Weight { get; }
public PointWeight(double point, double weight)
{
Point = point;
Weight = weight;
}
}
class GaussianQuadratureIntegrator : Integrator
{
static readonly PointWeight[] TwoPointsRuleData =
{
new PointWeight(-Math.Sqrt(1.0 / 3.0), 1.0),
new PointWeight( Math.Sqrt(1.0 / 3.0), 1.0),
};
static readonly PointWeight[] ThreePointsRuleData =
{
new PointWeight(-Math.Sqrt(3.0 / 5.0), 5.0 / 9.0),
new PointWeight( 0.0, 8.0 / 9.0),
new PointWeight( Math.Sqrt(3.0 / 5.0), 5.0 / 9.0),
};
static readonly PointWeight[] FivePointsRuleData =
{
new PointWeight(-0.9061798459386640, 0.2369268850561891),
new PointWeight(-0.5384693101056831, 0.4786286704993665),
new PointWeight( 0.0, 0.5688888888888889),
new PointWeight( 0.5384693101056831, 0.4786286704993665),
new PointWeight( 0.9061798459386640, 0.2369268850561891),
};
static readonly PointWeight[] SevenPointsRuleData =
{
new PointWeight(-0.9491079123427585, 0.1294849661688697),
new PointWeight(-0.7415311855993944, 0.2797053914892767),
new PointWeight(-0.4058451513773972, 0.3818300505051189),
new PointWeight( 0.0, 0.4179591836734694),
new PointWeight( 0.4058451513773972, 0.3818300505051189),
new PointWeight( 0.7415311855993944, 0.2797053914892767),
new PointWeight( 0.9491079123427595, 0.1294849661688697),
};
public static readonly IReadOnlyCollection<int> ValidPointCounts = new List<int> { 2, 3, 5, 7 }.AsReadOnly();
PointWeight[] RuleData { get; }
public GaussianQuadratureIntegrator(int pointsCount)
{
switch (pointsCount)
{
case 2: RuleData = TwoPointsRuleData; break;
case 3: RuleData = ThreePointsRuleData; break;
case 5: RuleData = FivePointsRuleData; break;
case 7: RuleData = SevenPointsRuleData; break;
default: throw new ArgumentOutOfRangeException(nameof(pointsCount));
}
}
public double Integrate(Func<double, double> f, double a, double b)
{
double center = (a + b) / 2;
double halfLength = (b - a) / 2;
var result = 0.0;
for (int i = 0; i < RuleData.Length; i++)
{
result += RuleData[i].Weight * f(center + halfLength * RuleData[i].Point);
}
result *= halfLength;
return result;
}
}
class DumbRiemannSummator : Integrator
{
int NumberOfSegments { get; }
public DumbRiemannSummator(int numberOfSegments)
{
NumberOfSegments = numberOfSegments;
}
public double Integrate(Func<double, double> f, double a, double b)
{
var result = 0.0;
for (int i = 0; i < NumberOfSegments; i++)
{
result += f(a + i * (b - a) / NumberOfSegments);
}
result /= NumberOfSegments;
return result;
}
}
class AdaptiveIntegrator : Integrator
{
Integrator PreciseIntegrator { get; }
Integrator CrudeIntegrator { get; }
double MaxError { get; }
int MaxSubsegments { get; }
public int SubsegmentCount { get; private set; }
public double ErrorEstimate { get; private set; }
public AdaptiveIntegrator(Integrator preciseIntegrator, Integrator crudeIntegrator, double maxError, int maxSubsegments)
{
PreciseIntegrator = preciseIntegrator;
CrudeIntegrator = crudeIntegrator;
MaxError = maxError;
MaxSubsegments = maxSubsegments;
SubsegmentCount = 0;
ErrorEstimate = double.NaN;
}
struct EstimationResult {
public double Result { get; }
public double Error { get; }
public EstimationResult(double result, double error)
{
Result = result;
Error = error;
}
}
struct Subsegment
{
public double Left { get; }
public double Right { get; }
public Subsegment(double left, double right)
{
Left = left;
Right = right;
}
}
EstimationResult EstimateIntegral(Func<double, double> f, Subsegment s)
{
var result = PreciseIntegrator.Integrate(f, s.Left, s.Right);
var crudeResult = CrudeIntegrator.Integrate(f, s.Left, s.Right);
var error = Math.Abs(result - crudeResult);
return new EstimationResult(result, error);
}
static int FindIndexOfMaxError(List<EstimationResult> estimations)
{
var maxError = double.NegativeInfinity;
var result = -1;
for (var index = 0; index < estimations.Count; index++)
{
if (maxError < estimations[index].Error)
{
result = index;
maxError = estimations[index].Error;
}
}
return result;
}
public double Integrate(Func<double, double> f, double a, double b)
{
var subsegment = new Subsegment(a, b);
var estimation = EstimateIntegral(f, subsegment);
SubsegmentCount = 1;
ErrorEstimate = estimation.Error;
var subsegments = new List<Subsegment> { subsegment };
var estimations = new List<EstimationResult> { estimation };
if (estimation.Error < MaxError)
{
return estimation.Result;
}
// Loop invariant: there are exactly SubsegmentCount subsegments and estimations.
// Loop exit condition: there are exactly MaxSubsegments subsegments and estimations.
while (SubsegmentCount < MaxSubsegments)
{
var maxErrorIndex = FindIndexOfMaxError(estimations);
var dividedSubsegment = subsegments[maxErrorIndex];
var subsegmentCenter = (dividedSubsegment.Left + dividedSubsegment.Right) / 2;
var leftSubsegment = new Subsegment(dividedSubsegment.Left, subsegmentCenter);
var rightSubsegment = new Subsegment(subsegmentCenter, dividedSubsegment.Right);
var leftEstimation = EstimateIntegral(f, leftSubsegment);
var rightEstimation = EstimateIntegral(f, rightSubsegment);
subsegments[maxErrorIndex] = leftSubsegment;
estimations[maxErrorIndex] = leftEstimation;
subsegments.Add(rightSubsegment);
estimations.Add(rightEstimation);
var totalError = estimations.Sum(x => x.Error);
SubsegmentCount++;
ErrorEstimate = totalError;
if (totalError < MaxError)
{
var totalResult = estimations.Sum(x => x.Result);
return totalResult;
}
}
throw new Exception("Exceeded subsegments limit");
}
}
class InvocationCounter<TArg, TResult>
{
private Func<TArg, TResult> Func { get; }
public int Counter { get; private set; }
public InvocationCounter(Func<TArg, TResult> func)
{
Func = func;
Counter = 0;
}
public TResult Invoke(TArg arg)
{
Counter++;
return Func(arg);
}
}
class Program
{
static double F(double x)
{
return Functions.NormalDistributionDensity(x, 0.75, 0.09) * Functions.WeirdFunction(x);
}
static void TryIntegrator(string name, Integrator integrator)
{
const double a = 0.0;
const double b = 1.0;
var invoker = new InvocationCounter<double, double>(F);
var result = integrator.Integrate(invoker.Invoke, a, b);
Console.Write("Method '{0}': result = {1}, after {2} evaluations of F", name, result, invoker.Counter);
var adaptiveIntegrator = integrator as AdaptiveIntegrator;
if (adaptiveIntegrator != null)
{
Console.Write(", with {0} subsegments, error estimation is {1}",
adaptiveIntegrator.SubsegmentCount, adaptiveIntegrator.ErrorEstimate);
}
Console.WriteLine();
}
static void Main(string[] args)
{
const int RIEMANN_SEGMENTS_COUNT = 1000;
const double MAX_ERROR = 1e-5;
const int MAX_SUBSEGMENTS = 50;
TryIntegrator($"Dumb Riemann summator ({RIEMANN_SEGMENTS_COUNT})", new DumbRiemannSummator(RIEMANN_SEGMENTS_COUNT));
foreach (var pointCount in GaussianQuadratureIntegrator.ValidPointCounts)
{
TryIntegrator($"{pointCount}-point Gaussian quadrature", new GaussianQuadratureIntegrator(pointCount));
}
TryIntegrator($"Adaptive Gaussian quadrature with 5- and 7-points rules",
new AdaptiveIntegrator(
new GaussianQuadratureIntegrator(7),
new GaussianQuadratureIntegrator(5),
MAX_ERROR, MAX_SUBSEGMENTS));
Console.ReadLine();
}
}
}
Sign up for free to join this conversation on GitHub. Already have an account? Sign in to comment
You can’t perform that action at this time.