C# Two Variable Integration Example

using System;

using CenterSpace.NMath.Core;
using CenterSpace.NMath.Analysis;

namespace CenterSpace.NMath.Analysis.Examples.CSharp
{
  class TwoVariableIntegrationExample
  {
    /// <summary>
    /// A .NET example in C# showing how to integrate over two variables.
    /// </summary>
    static void Main(string[] args)
    {
      Console.WriteLine();

      // Integrate the function F over the area where 0 < x < 1 and 0 < y < 3

      // First create the integrand
      MultiVariableFunction integrand = new MultiVariableFunction(new NMathFunctions.DoubleVectorDoubleFunction(F));

      // Create the integrator with defaults
      TwoVariableIntegrator integrator = new TwoVariableIntegrator();

      Console.WriteLine("Calculating the integral of ( x + 10 ) * y^2");
      Console.WriteLine("where 0 < x < 1 and 0 < y < 3");
      Console.WriteLine();

      double integral = integrator.Integrate(integrand, 0, 1, 0, 3);
      Console.WriteLine("integral is... " + integral);
      Console.WriteLine();

      // What if one of the bounds is not constant but a function? 
      // Here the upper bound of y is 3x + 5.
      OneVariableFunction yUpper = new OneVariableFunction(new NMathFunctions.DoubleUnaryFunction(Y));

      // Integrate the function, F, over the area where 0 < x < 1 and 0 < y < ( 3x + 5 )
      Console.WriteLine("where 0 < x < 1 and 0 < y < 3x + 5");
      Console.WriteLine();

      integral = integrator.Integrate(integrand, 0, 1, 0, yUpper);
      Console.WriteLine("integral is... " + integral);
      Console.WriteLine();

      Console.WriteLine("switch to Romberg");
      Console.WriteLine();
      integrator.DxIntegrator = new RombergIntegrator();
      integrator.DyIntegrator = new RombergIntegrator();
      integral = integrator.Integrate(integrand, 0, 1, 0, yUpper);
      Console.WriteLine("integral is... " + integral);
      Console.WriteLine();

      Console.WriteLine();
      Console.WriteLine("Press Enter Key");
      Console.Read();
    }

    private static double F(DoubleVector v)
    {
      return (v[0] + 10) * (v[1] * v[1]);
    }

    private static double Y(double x)
    {
      return (3 * x) + 5;
    }

  }  // class

} // namespace