Illustrates numerical integration of functions in higher dimensions using classes in the Extreme.Mathematics.Calculus namespace in 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