|
Posted on May 8, 2007 @ 10:29:09 PM by Paul Meagher
Understanding Simpson's rule can be difficult if you are laboring under the assumption that there is only one version of it. In fact, there are at least 4 basic versions of Simpson's rule: quadratic, cubic, adaptive, and composite. I implemented each of these variants as methods in a class file called Math_Integration_Simpson.php.
<?php
/**
* @package Math_Integration
* @filepath Math/Integration/Simpson.php
*
* A numerical intergration class implementing 4 variants of
* Simpson's rule.
*
* Rule of thumb:
*
* If the function is smooth, Simpson's rule is better.
* If the function has abrupt changes, then Trapezoid is better.
*
* @author Paul Meagher
* @version 0.1
* @license LGPL
* @modified May 10, 2007
*
* @see http://math.fullerton.edu/mathews/n2003/NewtonCotesMod.html
* @see http://en.wikipedia.org/wiki/Adaptive_Simpson's_method
* @see http://www.cs.queensu.ca/~daver/271/Notes.pdf
*/
class Math_Integration_Simpson {
/**
* Approximation of f by a quadratic polynomial (i.e., uses three
* weighted points). x0 as a, x2 as b and x1 as 1/2 distance
* between a and b. This is the most basic form of Simpson's
* rule and a favorite among textbook authors - also known as
* Simpson's one-third rule.
*/
public function quadratic($f, $a, $b) {
$c = ($a+$b)/2.0;
$h = abs($b-$a)/6.0;
return $h*($f->valueAt($a) + 4.0*$f->valueAt($c) + $f->valueAt($b));
}
/**
* Approximation f by a cubic polynomial (i.e., uses four
* weighted points). x0 as a, x3 as b and x1 and x2 as 1/3 and
* 2/3 distances between a and b. Also known as Simpson's
* three-eights rule.
*/
public function cubic($f, $a, $b) {
$h = abs($b-$a)*(1.0/3.0);
$c = $a + $h;
$d = $b - $h;
return (3.0*$h/8.0)*($f->valueAt($a) + 3.0*$f->valueAt($c)+ 3.0*$f->valueAt($d)+ $f->valueAt($b));
}
/**
* Recursive implementation of adaptive Simpson's rule.
*/
private function recursive($f, $a, $b, $eps, $sum) {
$c = ($a+$b)/2.0;
$left = $this->quadratic($f, $a, $c);
$right = $this->quadratic($f, $c, $b);
if (abs($left + $right - $sum) <= (15 * $eps))
return $left + $right + ($left + $right - $sum)/15;
return $this->recursive($f, $a, $c, $eps/2, $left) + $this->recursive($f, $c, $b, $eps/2, $right);
}
/**
* Calculate integral of f from a to b using adaptive Simpson's rule
* with max error of eps.
*/
public function adaptive($f, $a, $b, $eps=0.000001) {
return $this->recursive($f, $a, $b, $eps, $this->quadratic($f, $a, $b));
}
/**
* Calculate integral of f from a to b using a piecewise approach. This
* form is commonly used to estimate the integral with $n used to control
* the accuracy of the estimate.
*/
public function composite($f, $a, $b, $n=10) {
$h = ($b-$a)/(2.0*$n);
$s1 = ($f->valueAt($a)+ $f->valueAt($b));
for($i=1; $i<2*$n; $i=$i+2)
$s2 = $s2 + $f->valueAt($a + $i*$h);
for($i=2; $i<2*$n; $i=$i+2)
$s3 = $s3 + $f->valueAt($a + $i*$h);
return ($h/3.0)*($s1+ 4.0*$s2 + 2.0*$s3);
}
}
?>
The following simpson_methods.php script executes each of the methods defined in this class.
<?php
include "Math/Integration/Simpson.php";
class F {
function valueAt($x) {
return 4.0/(1.0+$x*$x);
}
function toString() {
return "4.0/(1.0+x^2)";
}
function toHTML() {
return "4.0/(1.0+x<sup>2</sup>)";
}
}
$f = new F;
$a = 0;
$b = 1;
$simpson = new Math_Integration_Simpson;
$quadratic = $simpson->quadratic($f, $a, $b);
$cubic = $simpson->cubic($f, $a, $b);
$adaptive = $simpson->adaptive($f, $a, $b);
$composite = $simpson->composite($f, $a, $b);
echo "quadratic method: $quadratic <br />";
echo "cubic method: $cubic <br />";
echo "adaptive method: $adaptive <br />";
echo "composite method: $composite <br />";
?>
Each method in the script generates a π estimate:
quadratic method: 3.13333333333
cubic method: 3.13846153846
adaptive method: 3.14159265371
composite method: 3.14159265297
The method producing the output closest to an exact π value of 3.1415926535 is the adaptive method, beating out the composite method by a tiny margin.
Other classes that might be found in a Math_Integration package include:
Math_Integration_Trapezoid.phpMath_Integration_Gaussian.php (i.e., gaussian quadrature)Math_Integration_Romberg.php
All would have basic, adaptive, and composite variants.
Exercises
1. Play around with the $eps parameter in the adaptive Simpson's rule and the $n parameter in the composite Simpson's rule. Can the composite method be made to out perform the adaptive method in computing π?
2. Apply these methods to a variety of other functions and note any
differences.
3. Implement a Math_Integration_Trapezoid.php class that contains a basic, adaptive, and composite implementation of the trapezoid method of integration.
4. Compare the relative performance of Simpson's rule to the Trapezoid method of integration on functions that are smooth or contain abrupt changes.
5. How do Lagrange polynomials relate to Simpson's rule?
6. How do the Newton-Cotes formulas relate to Simpson's rule and other methods of integration?
|
