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

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;
        }
    }
}