Higher Dimensional Numerical Integration in IronPython QuickStart Sample

Illustrates numerical integration of functions in higher dimensions using classes in the Numerics.NET.Calculus namespace in IronPython.

This sample is also available in: C#, Visual Basic, F#.

Overview

This QuickStart sample demonstrates how to perform numerical integration of functions in two dimensions using Numerics.NET’s integration capabilities.

The sample shows several key aspects of multi-dimensional integration:

  • Using the AdaptiveIntegrator2D class for efficient integration over rectangular regions
  • Working with the Repeated1DIntegrator2D class that uses repeated one-dimensional integration
  • Setting up integrand functions using delegates
  • Integrating over arbitrary non-rectangular regions by specifying boundary functions
  • Monitoring integration progress through properties like Status, EstimatedError, and function evaluation counts
  • Comparing numerical results with known analytical solutions

The example includes integration of several test functions including:

  • A rational function: 4/(1 + 2x + 2y) over the unit square
  • A trigonometric function: (π²/4)sin(πx)sin(πy) over the unit square
  • A polynomial function: x²y² over the unit disk

For each integration, the sample demonstrates how to set up the problem, execute the integration, and analyze the results including error estimates and computational effort required.

The code

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