Higher Dimensional Numerical Integration QuickStart Sample (IronPython)

import numerics

from math import *

# The numerical integration classes reside in the
# Extreme.Mathematics.Calculus namespace.
from Extreme.Mathematics.Calculus import *
# Function delegates reside in the Extreme.Mathematics
# namespace.
from Extreme.Mathematics import *

#/ Illustrates numerical integration in higher dimensions using
#/ classes in the Extreme.Mathematics.Calculus namespace of the Extreme
#/ Optimization Numerical Libraries for .NET.

#
# Two-dimensional integration
#

# The function we are integrating must be
# provided as a Func<double, double, double> delegate. 

# The AdaptiveIntegrator2D class is the most efficient
# 2D integrator in most cases. It uses an adaptive algorithm.

# Construct an instance of the integrator class:
integrator1 = AdaptiveIntegrator2D()

# An example of setting the integrand and bounds through properties
# is given below. Here, we put the integrand and the bounds 
# of the integration region directly in the call to Integrate, # which performs the calculation:
integrator1.Integrate(lambda x, y: 4 / (1 + 2 * x + 2 * y), 0, 1, 0, 1)
print "4 / (1 + 2x + 2y) on [0,1] * [0,1]"
print "  Value:       {0:.15f}".format(integrator1.Result)
print "  Exact value: {0:.15f} = Ln(3125 / 729)".format(log(3125.0 / 729.0))
# To see whether the algorithm ended normally, # inspect the Status property:
print "  Status:", integrator1.Status
print "  Estimated error:", integrator1.EstimatedError
print "  Iterations:", integrator1.IterationsNeeded
print "  Function evaluations:", integrator1.EvaluationsNeeded

# Another integrator uses repeated 1-dimensional
# integration:
integrator2 = Repeated1DIntegrator2D()

# You can set the order of integration, as well as
# the integration rules for the X and the Y direction:
integrator2.InitialDirection = Repeated1DIntegratorDirection.X

# You can set the integrand and the bounds of the integration region
# by setting properties of the integrator object:
integrator2.Integrand = lambda x, y: pi ** 2 / 4.0 * sin(pi * x) * sin(pi * y)
integrator2.XLowerBound = 0
integrator2.XUpperBound = 1
integrator2.YLowerBound = 0
integrator2.YUpperBound = 1

result = integrator2.Integrate()
print "Pi^2 / 4 Sin(Pi x) Sin(Pi y)   on [0,1] * [0,1]"
print "  Value:       {0:.15f}".format(integrator2.Result)
print "  Exact value: {0:.15f}".format(1.0)
# To see whether the algorithm ended normally, # inspect the Status property:
print "  Status:", integrator2.Status
print "  Estimated error:", integrator2.EstimatedError
print "  Iterations:", integrator2.IterationsNeeded
print "  Function evaluations:", integrator2.EvaluationsNeeded

#
# Integration over arbitrary regions
#

# The repeated 1D integrator can also be used to compute
# integrals over arbitrary regions. To do this, you need to
# supply function that return the lower bound and upper bound 
# of the region as a function of x.

# Here, we integrate x^2 * y^2 over the unit disk.
integrator2.LowerBoundFunction = lambda x: 0.0 if abs(x) >= 1.0 else -sqrt(1.0 - x*x)
integrator2.UpperBoundFunction = lambda x: 0.0 if abs(x) >= 1.0 else sqrt(1.0 - x*x)
integrator2.XLowerBound = -1
integrator2.XUpperBound = 1

integrator2.Integrand = lambda x, y: x ** 2 * y ** 2

result = integrator2.Integrate()
print "x^2 * y^2 on the unit disk"
print "  Value:       {0:.15f}".format(integrator2.Result)
print "  Exact value: {0:.15f} = Pi / 24".format(pi / 24)
# To see whether the algorithm ended normally, # inspect the Status property:
print "  Status:", integrator2.Status
print "  Estimated error:", integrator2.EstimatedError
print "  Iterations:", integrator2.IterationsNeeded
print "  Function evaluations:", integrator2.EvaluationsNeeded