Basic Integration in IronPython QuickStart Sample

Illustrates the basic numerical integration classes in IronPython.

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

Overview

This QuickStart sample demonstrates how to use the basic numerical integration classes in Numerics.NET.

The sample shows how to use three different non-adaptive numerical integration algorithms:

• Simpson’s Rule - A simple and widely used integration method
• Gauss-Kronrod Integration - An optimized fixed-point scheme using carefully chosen evaluation points
• Romberg Integration - A method that combines Simpson’s Rule with convergence acceleration

For each integration method, the sample demonstrates:

• How to configure the integrator (tolerances, convergence criteria)
• How to perform the actual integration
• How to access results and diagnostic information (error estimates, iteration counts)
• The strengths and limitations of each method

The sample concludes by showing how simpler integration methods can break down when dealing with difficult integrands that have singularities or discontinuities, setting up the motivation for using adaptive integration methods covered in other samples.

The code

import numerics

from math import *

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

# Illustrates the basic use of the numerical integration
# classes in the Extreme.Mathematics.Calculus namespace of the Extreme
# Optimization Mathematics Library for .NET.

# Numerical integration algorithms fall into two
# main categories: adaptive and non-adaptive.
# This QuickStart Sample illustrates the use of
# the non-adaptive numerical integrators.
#
# All numerical integration classes derive from
# NumericalIntegrator. This abstract base class
# defines properties and methods that are shared
# by all numerical integration classes.

#
# The integrand
#

# The function we are integrating must be
# provided as a Func<double, double>. 
f = sin
			
#
# SimpsonIntegrator
# 

# The simplest numerical integration algorithm
# is Simpson's rule. 
simpson = SimpsonIntegrator()
# You can set the relative or absolute tolerance
# to which to evaluate the integral.
simpson.RelativeTolerance = 1e-5
# You can select the type of tolerance using the
# ConvergenceCriterion property:
simpson.ConvergenceCriterion = ConvergenceCriterion.WithinRelativeTolerance
# The Integrate method performs the actual 
# integration:
result = simpson.Integrate(sin, 0, 2)
print "sin(x) on [0,2]"
print "Simpson integrator:"
# The result is also available in the Result 
# property:
print "  Value:", simpson.Result
# To see whether the algorithm ended normally, # inspect the Status property:
print "  Status:", simpson.Status
# You can find out the estimated error of the result
# through the EstimatedError property:
print "  Estimated error:", simpson.EstimatedError
# The number of iterations to achieve the result
# is available through the IterationsNeeded property.
print "  Iterations:", simpson.IterationsNeeded
# The number of function evaluations is available 
# through the EvaluationsNeeded property.
print "  Function evaluations:", simpson.EvaluationsNeeded

#
# Gauss-Kronrod Integration
#

# Gauss-Kronrod integrators also use a fixed point 
# scheme, but with certain optimizations in the 
# choice of points where the integrand is evaluated.

# The NonAdaptiveGaussKronrodIntegrator uses a
# succession of 10, 21, 43, and 87 point rules
# to approximate the integral.
nagk = NonAdaptiveGaussKronrodIntegrator()
nagk.Integrate(sin, 0, 2)
print "Non-adaptive Gauss-Kronrod rule:"
print "  Value:", nagk.Result
print "  Status:", nagk.Status
print "  Estimated error:", nagk.EstimatedError
print "  Iterations:", nagk.IterationsNeeded
print "  Function evaluations:", nagk.EvaluationsNeeded
			
#
# Romberg Integration
#

# Romberg integration combines Simpson's Rule
# with a scheme to accelerate convergence.
# This algorithm is useful for smooth integrands.
romberg = RombergIntegrator()
result = romberg.Integrate(sin, 0, 2)
print "Romberg integration:"
print "  Value:", romberg.Result
print "  Status:", romberg.Status
print "  Estimated error:", romberg.EstimatedError
print "  Iterations:", romberg.IterationsNeeded
print "  Function evaluations:", romberg.EvaluationsNeeded
			
# However, it breaks down if the integration
# algorithm contains singularities or 
# discontinuities.
#
# The AdaptiveIntegrator can handle this type
# of integrand, and many other difficult cases.
# See the AdvancedIntegration QuickStart sample
# for details.
result = romberg.Integrate(lambda x: 0.0 if x <= 0.0 else x ** -0.9 * log(1/x), 0.0, 1.0)
print "Romberg on hard integrand:"
print "  Value:", romberg.Result
print "  Actual value: 100"
print "  Status:", romberg.Status
print "  Estimated error:", romberg.EstimatedError
print "  Iterations:", romberg.IterationsNeeded
print "  Function evaluations:", romberg.EvaluationsNeeded

# Function that will cause difficulties to the
# simplistic integration algorithms.
# </summary>
def HardIntegrand(x):
	# This is put in because some integration rules
	# evaluate the function at x=0.
	if x <= 0:
		return 0
	return x**-0.9 * log(1/x)