math-suite/Source/GraphingCalculator/Expression/Integrator.cs

117 lines
3.7 KiB
C#
Raw Normal View History

2018-02-05 23:24:46 +00:00
using System;
using System.Collections.Generic;
using System.Linq;
using System.Text;
namespace GraphingCalculator
{
public static class Integrator
{
#region Settings
private static int maximumRecursionDepth = 16;
public static int MaximumRecursionDepth
{
get { return maximumRecursionDepth; }
set { maximumRecursionDepth = value; }
}
private static double initialEpsilon = 1e-10;
public static double InitialEpsilon
{
get { return initialEpsilon; }
set { initialEpsilon = value; }
}
#endregion
#region Some variables
private static string variable = "x";
private static Expression expression;
private static Dictionary<double, double> results = new Dictionary<double, double>();
#endregion
#region Swap function
private static void Swap(ref double a, ref double b)
{
double c = a;
a = b;
b = c;
}
#endregion
private static double EvaluateCache(double x)
{
if (!results.ContainsKey(x))
{
expression.Variables[variable] = x;
results[x] = expression.Evaluate();
}
return results[x];
}
private static double CalculateIntegralAux(double xA, double xB, double epsilon, double S, double yA, double yB, double yC, int depth)
{
double xC = (xA + xB) / 2;
double width = (xB - xA);
double xD = (xA + xC) / 2;
double xE = (xC + xB) / 2;
double yD = EvaluateCache(xD);
double yE = EvaluateCache(xE);
double Sleft = (width / 12) * (yA + 4 * yD + yC);
double Sright = (width / 12) * (yC + 4 * yE + yB);
double S2 = Sleft + Sright;
if (depth <= 0 || Math.Abs(S2 - S) <= 15 * epsilon)
return S2 + (S2 - S) / 15;
return CalculateIntegralAux(xA, xC, epsilon / 2, Sleft, yA, yC, yD, depth - 1) +
CalculateIntegralAux(xC, xB, epsilon / 2, Sright, yC, yB, yE, depth - 1);
}
private static double CalculateIntegral(double xA, double xB, double epsilon, int maxDepth)
{
double xC = (xA + xB) / 2;
double width = xB - xA;
double yA = EvaluateCache(xA);
double yC = EvaluateCache(xC);
double yB = EvaluateCache(xB);
double S = (width / 6) * (yA + 4 * yC + yB);
return CalculateIntegralAux(xA, xB, epsilon, S, yA, yB, yC, maxDepth);
}
/// <summary>
/// Calculates a definite integral (the area between the function graph and the X axis, where x ranges between beg and end.
/// </summary>
/// <param name="expr">An expression to integrate.</param>
/// <param name="beg">Beginning of the interval</param>
/// <param name="end">End of the interval</param>
/// <param name="var">The name of the variable, by default is 'x'</param>
public static double Integrate(Expression expr, double beg, double end, string var = "x")
{
double result = 0;
bool changeSign = false;
// Make sure beg < end
if (beg > end)
{
Swap(ref beg, ref end);
changeSign = true;
}
// Initialization
results.Clear();
variable = var;
expression = expr;
// Calculation
result = CalculateIntegral(beg, end, initialEpsilon, maximumRecursionDepth);
if (changeSign) return -result;
return result;
}
}
}