php/Math   
Recreational Mathematics   
   home  |  library  |  contact
 Math Notes
 Math Programming [25]
 Regression [3]
 Data Mining [17]
 Notation [6]
 Linear Algebra [9]
 Stats & Prob [15]
 Math Cognition [5]
 Space & Physics [6]
 Formulas [5]
 Fun & Games [2]
 Haskell [1]
 Bayes Theory [1]
 Site News [0]
 Math Projects [5]
 Polynomials [1]
 Calculus [9]
 Number Theory [3]
 Optimization [2]
 Financial [1]

 Math Links
 PHP/ir
 Andrew Gelman
 Chance Wiki
 Daniel Lemire
 KD Knuggets
 Social Stats
 MySQL Performance
 Hunch.net
 Matthew Hurst
 JMLR
 JSS
 Hal Daume III
 Math Notes >> Calculus

On the way to NURBS [Calculus
Posted on January 28, 2008 @ 02:03:17 AM by Paul Meagher

NURBS is an acronym for "Non-Uniform Rational B-Spline". NURBS are used extensively in CAD systems to represent surfaces. I hope to someday understand the math and programming behind NURBS better. Towards that end, I have ported some code that manipulates Bernstein polynomials. Bernstein polynomials unify the math behind Bezier approximation schemes and Spline approximation schemes. NURBS are another approximation scheme which requires that you know something about Bernstein polynomials in order to appreciate the math involved.

The Bernstein class below can be used to generate the coordinates needed to draw a curve that interpolates and approximates a set of control points. Vector-based image manipulation programs offer a Bezier tool that allows you to create curves by moving around a set of bounding control points. As you move the curve control points, the shape of the Bezier curve changes to interpolate and approximate the new position of the control points. The class below implements the math required to interpolate and approximate the position of the control points. The math looks like this:

The implementation of this math looks like this:

<?php
/**
* Provides methods to manipulate the Bernstein polynomial.
*
* An nth-degree Berstein polynomial is given by:
*
* B(i,n,u) = n! * u^i * (1-u)^(n-i) / i! * (n-i)!
* with  0 <= i <= n
*
* @author R.P. from OpaleTeam
* @see http://opale.belios.de
* @license LPGL
* @date 05/2002
*
* Ported by Paul Meagher on Jan 20, 2008.
*/
class Bernstein {

  
/**
  * Returns the evaluation, for a given value, of a Bernstein
  * polynomal defined by the given arguments.
  *
  * @param int $i 
  * @param int $n the degree of the polynomial
  * @param double $u the value for which the polynomial will be computed
  * @return double the computed value
  */
  
public function getValue($i$n$u) {
  
    
$temp[$n-$i] = 1.0;
    
$u1          1.0 $u;

    for (
$k 1$k <= $n$k++) 
      for (
$j $n$j >= $k$j--) 
        
$temp[$j] = $u1 $temp[$j] + $u $temp[$j-1];       
    
    return 
$temp[$n];

  }
  
  
/**
  * Returns all Bernstein polynomials evaluations, for a given
  * value, of a Bernstein polynomal defined by the given arguments.
  *
  * @param int $n the degree of the polynomial
  * @param double $u the value for which the polynomial will be computed
  * @return double[], the computed value
  */
  
public function getAllValues($n$u) {

    
$b[0] = 1.0;
    
$u1   1.0 $u;

    for (
$i 1$i <= $n$i++) {
      
$saved 0.0;
      for (
$k 0$k $i$k++) {
        
$temp  $b[$k];
        
$b[$k] = $saved $u1 $temp;
        
$saved $u $temp;
      }
      
$b[$i] = $saved;
    }
    
    return 
$b;

  }  
  

}

?>

The bernstein_test.php example below illustrates usage. In this example, we use the Bernstein class to find the coefficients to multiply the $X[$i] (i.e., $X = array(1, 2, 4, 3)) and $Y[$i] coordinates (i.e., $Y = array(1, 3, 3, 1)) of 4 control points. After we separately multiply the $X[$i] and $Y[$i] terms by the Bernstein cofficients we get a new point at $t that approximates the control points.

The setup value of $n is equal to 1 minus the number of control points. $T is an array of values in the range [0,1] that identify the left (0) to right (1) positions along the approximating and interpolating curve where we want to compute it's coordinates.

<?php
/**
* Testing Bernstein output against worked example in:
*
* Rogers, David F. An Introduction to NURBS With Historical Perspective.  
* 2001. Morgan Kaufmann, San Francisco. pp. 22-24
*/

include "Bernstein.php";

$berstein = new Bernstein;

$n 3;

$X = array(1243);
$Y = array(1331);

?>
<table border='1' cellspacing='1' cellpadding='2'>
  <tr>
    <th>t</th>
    <th>B<sub>3,0</sub></th>
    <th>B<sub>3,1</sub></th>
    <th>B<sub>3,2</sub></th>
    <th>B<sub>3,3</sub></th>
    <th>P(t)</th>   
  </tr>
  <?php
  $T 
= array(00.150.350.500.650.851);  
  foreach (
$T as $t) {
    
$B $berstein->getAllValues($n$t);
    
?>
    <tr>
      <td align='center'><?php echo $t ?></td>
      <?php
      $i
=0;
      
$x=0.0;
      
$y=0.0;     
      foreach (
$B as $b) {
        
?>
        <td align='center'><?php echo $b ?></td>
        <?php
        $x 
+= $b *$X[$i];
        
$y += $b *$Y[$i];       
        
$i++;       
      }
      
?>
      <td align='center'><?php echo "($x, $y)"?></td>
    </tr>
    <?php
  

  
?>
</table>

The test script above produces the following HTML table as output:

t B3,0 B3,1 B3,2 B3,3 P(t)
0 1 0 0 0 (1, 1)
0.15 0.614125 0.325125 0.057375 0.003375 (1.504, 1.765)
0.35 0.274625 0.443625 0.238875 0.042875 (2.246, 2.365)
0.5 0.125 0.375 0.375 0.125 (2.75, 2.5)
0.65 0.042875 0.238875 0.443625 0.274625 (3.119, 2.365)
0.85 0.003375 0.057375 0.325125 0.614125 (3.261, 1.765)
1 0 0 0 1 (3, 1)

The P(t) column above contains a set of points which we can graph that will approximate the set of control points a user has specified for the curve they want drawn. The x coordinate of each point was obtained by multiplying the 4 Berstein coefficients by the 4 x coordinates. Similarly, the y coordinate of each point was obtained by multiplying the 4 Bernstein coefficients by the 4 y coordinates.

Permalink 

Monte Carlo Integration [Calculus
Posted on January 9, 2008 @ 03:09:37 AM by Paul Meagher

Monte carlo integration is a way to find the area under a curve using random numbers. The program below provides a nice introductory example of a monte carlo method and how it can be used to solve a calculus problem.

<?php
/**
* Approximates area of the region under the graph of a
* monotonic function and above the x-axis beween a and b.
*
* Adaptation of Monte Carlo Method program by James M. Sconyers
* found in:
*
* Larson, Hostetler, Edwards (2006). Calculus I with Precalculus.
* 2nd Edition. Houghton Mifflin Co, Boston.  p. 490
*/
class MonteCarlo_Integration {

  
/**
  * @see http://ca.php.net/manual/en/function.mt-rand.php#75793
  */
  
function random_float($min$max) {
    return (
$min+lcg_value()*(abs($max-$min)));
  }
  
  
/**
  * Compute area under the curve using monte carlo method.
  * 
  * @param object $f function object
  * @param float $a lower limit
  * @param float $b upper limit
  * @param int $n number of samples
  * @return float area under the curve
  */
  
function integrate($f$a$b$n) {
  
    
$p 0;
  
    if (
$f->valueAt($a) > $f->valueAt($b))
      
$ymax $f->valueAt($a);
    else 
      
$ymax $f->valueAt($b);
    
    for (
$i=1$i $n$i++) {
      
$x $a + ($b $a) * $this->random_float(01);
      
$y $ymax $this->random_float(01);
      if (
$y $f->valueAt($x)) 
        
$p $p+1;
    }
     
    
$area = ($p/$n)*($b-$a)*$ymax;
    
    return 
$area;
  
  }

}   

?>    

Here is a program to test the MonteCarlo_Integration.php class:

<?php

include "MonteCarlo_Integration.php";

class 
{    
    function 
valueAt($x) {
    return 
pow($x2);
  }
}


$f = new F;
$a 0;
$b 2;
$n 10000;

$mc   = new MonteCarlo_Integration;
$area $mc->integrate($f$a$b$n);

echo 
"Approximate Area: $area";

?> 

The output of one execution of the program is:

Approximate Area: 2.6632

The exact answer is 8/3.

One programming project would involve histogramming and summarizing the area estimates produced by the above class where $n, the number of samples, is varied. Doing so would give you a feel for how variable the estimates are with different sample sizes and how the monte carlo estimates are distributed. A project like this might offer some insight into more advanced monte carlo applications.

Permalink 

Differentiation with PHP (Part 7): Polynomial Modelling [Calculus
Posted on October 15, 2007 @ 02:37:38 PM by Paul Meagher

Why should we care about polynomials? Polynomials are easy to manipulate, understand, and they are very expressive. The expressiveness property of polynomials means that you can use polynomials to represent or formally model many domains of experience.

In part 7 of this blog series, I will discuss the calculation of retirement savings using a polynomial model. I will be assuming that you have read my previous blog so I don't have to go over some basic polynomial notation and ideas.

Retirement Savings: Math

In Connally, Hughes-Hallett, Gleason, et al. (2004). "Functions modelling change", John Wiley & Sons, p. 384, they pose and provide an answer to the following question:

You make five separate deposits of $1000 each into a savings account, one deposit per year, beginning today. What annual interest rate gives a balance in the account of $6000 five years from today? (Assume the interest rate is constant over these 5 years).

In answering this question, Connally et al. assume a 5% annual interest rate r, compounded annually, and use an induction proof to make the argument for a polynomial model of retirement savings. The induction proof begins with the simplest case; computing your balance after 1 year on the savings plan:

Balance = (100% of initial deposit) + (5% of Initial deposit) plus Second deposit
Balance = 105% of Initial deposit + Second deposit
Balance = 1000 * 1.05 + 1000
Balance = 1050 + 1000
Balance = 2050

We can represent this relationship a bit more abstractly as:

B1 = 1000x + 1000

Where B1 denotes the balance after one year on the plan, r denotes the annual interest rate, and x = 1 + r (i.e., 1 + 0.05). The savings balance equation B1 is a first order polynomial formula.

After two years on the plan, our savings balance B2 would be:

B2 = (1000x+1000)x + 1000
B2 = 1000x2 + 1000x + 1000 (second order polynomial formula)
B2 = 1000 * 1.05 * 1.05 + 1000 * 1.05 + 1000
B2 = 1000 * 1.1025 + 1050 * 1.05 + 1000
B2 = 1102.5 + 1050 + 1000
B2 = 3152.5

Observing the polynomial pattern in computing B1 and B2, we take the induction step and assert that all future instances of this computation can be done using a polynomial series with the same polynomial order as the year you want a balance for. Continuing on with our example, after 3 years, our balance would be computed with the following third order polynomial:

B3 = (1000x2 + 1000x +1000) x + 1000
B3 = 1000x3 + 1000x2 +1000x + 1000 (third order polynomial formula)

An so on for the 4th, 5th, etc... year of the savings plan.

Retirement Savings: PHP

So now that we know the math behind computing retirement savings, we can save ourselves some labor evaluating the polynomial for each year of our savings plan by using the Math_Polynomial package to do so.

<?php

include 'Math/Polynomial.php';
include 
'Math/PolynomialOp.php';

$x 1.05;

// lazy way of generating polynomial equations corresponding to each year
$EQ[0] = "1000";
$EQ[1] = "1000x + 1000";
$EQ[2] = "1000x^2 + 1000x + 1000";
$EQ[3] = "1000x^3 + 1000x^2 + 1000x + 1000";
$EQ[4] = "1000x^4 + 1000x^3 + 1000x^2 + 1000x + 1000";
$EQ[5] = "1000x^5 + 1000x^4 + 1000x^3 + 1000x^2 + 1000x + 1000";

$num_years count($EQ);

for(
$year=0$year $num_years$year++) {
  
$p = new Math_Polynomial($EQ[$year]);
  echo 
"<p>";
  echo 
"B<sub>$year</sub> = "$p->toString() ."<br />";
  echo 
"B<sub>$year</sub> evaluated at x = $x is equal to "Math_PolynomialOp::evaluate($p$x);
  echo 
"</p>";
}

?>

Which generates as output the savings balance for each year in a 5 year plan:

B0 = 1000
B0 evaluated at x = 1.05 is equal to 1000

B1 = 1000x + 1000
B1 evaluated at x = 1.05 is equal to 2050

B2 = 1000x^2 + 1000x + 1000
B2 evaluated at x = 1.05 is equal to 3152.5

B3 = 1000x^3 + 1000x^2 + 1000x + 1000
B3 evaluated at x = 1.05 is equal to 4310.125

B4 = 1000x^4 + 1000x^3 + 1000x^2 + 1000x + 1000
B4 evaluated at x = 1.05 is equal to 5525.63125

B5 = 1000x^5 + 1000x^4 + 1000x^3 + 1000x^2 + 1000x + 1000
B5 evaluated at x = 1.05 is equal to 6801.9128125

Conclusion

The script above does not demonstrate many of the capabilities of the Math_Polynomial package. This is because I spent most of my allotted time wrestling with the retirement savings math until 1) I recognized the form of argument as an induction proof and 2) I substituted example values into the calculation to confirm the plausibility of the polynomial series solution. The aim of this blog was to show that polynomial models are "expressive". To drive this point home I would have to provide more examples of polynomial modelling. The exercises below will hopefully substitute for the time which I do not have to show that polynomial models can be usefully applied to a variety of domains.

Exercises

1. We are half way done answering the original retirement savings question: What annual interest rate gives a balance in the account of $6000 five years from today?

Develop some numerical and graphical php-based scripts to compute your own answer.

2. Can you model the effect of food and drugs on the body using a polynomial model? If you need to maintain a certain level of a drug in a patient's body how do you compute the required dosages to maintain that level?

3. A popular numerical technique based on polynomial theory is cubic spline interpolation. Download JpGraph and experiment with the functionality it provides for cubic spline interpolation.

4. Explore the use of polynomials to imitate standard probability distributions (e.g., gaussian, binomial, poisson, etc...). How does the polynomial order and cooefficient weights affect the shape of the curve? Could the order and weight parameters of your polynomial model have a deeper relationship to the domain you are modelling than a probability model based upon a more standard probability distribution and its parameters?

Permalink 

Differentiation with PHP (Part 6): Architecture of the Math_Polynomial package [Calculus
Posted on October 2, 2007 @ 02:12:06 PM by Paul Meagher

In the next couple of blogs I will be discussing polynomial functions and software for manipulating them. We will begin by examining the relationship between power functions, polynomial functions, rational functions, and algrebraic functions. After we cover these typical calculus topics, we will ready to examine the principal data structure and code layout of Keith Palmer's Math_Polynomial package.

Power Functions

A power function is a function of the form:

f(x) = axn where a and n are constants.

To understand polynomial functions you need to understand how its graph changes when the power n of the function changes. When the power is 1 we have linear function which graphs as a straight line. When the power is 2 we have a quadratic function which graphs as a parabola. When the power is 3 we have a cubic function. The graph of a cubic function has an inflection point where the concavity of the curve changes. We can also talk about quartic, quintic and higher order powers functions; however, there is not much practical need to study power functions of higher order than a cubic function.

The derivative of a cubic function is a quadratic function. The derivative of a quadratic function is a linear function. The general rule is that the derivative of a power function is one order less than the original power function.

Polynomial Functions

Sums and differences of power functions lead to the family of polynomial functions:

f(x) = axn + bx(n-1) + cx(n-2) + ... + yx + z

Graphically, when you add the graph of a cubic function to the graph of a quadratic function you add the height of each graph at each point. This results in a flexible curve that is linear combination of the two generating curves. It is this flexibility that makes polynomials so important as they can be used to interpolate or fit very complex curves using relatively low order polynomicals such a cubic polynomial.

A cubic power function is different than a cubic polynomial function in that the later includes the sums and differences of lower order power functions. A power function only includes one term; a polynomial function includes more than one power function term.

Rational and Algebraic Functions

Ratios of polynomials lead to the family of rational functions. The term rational is used because rational functions are like rational numbers which are the ratio of two integers. Just as the ratio of two integers defines a larger set of numbers than just the integers, so to the ratio of two polynomials defines a larger set of functions than just the polynomials.

Because polynomials can be easily added, subtracted, multiplied and divided they belong to a still larger class of functions called the Algebraic functions. The fact that polynomials are algebraic functions leads to the expectation that a polynomial software package should have excellent support for manipulating polynomials algrebraically (i.e., support for adding, subtracting, dividing and multiplying polynomials).

Architecture of the Math_Polynomial package

Now that we know what polynomials are we are in a better position to understand the data structure and layout of Keith Palmer's Math_Polynomial package.

Lets start with the data structure, Math_PolynomialTerm, which is used to represent individal terms of a polyomial expression.

<?php

/**
 * Class representing a term of a Polynomial
 * 
 * Represents a polynomial term of the form: 
 *     Cx^E where C is a constant (float), and E is an exponent (int)
 * 
 * @see Math_Polynomial
 * 
 * @package Math_Polynomial
 * @author Keith Palmer <keith@UglySlug.com>
 */
class Math_PolynomialTerm
{
    
/**
     * The coefficient of the term
     * 
     * @var float
     * @access protected
     */
    
var $_coef;
    
    
/**
     * The exponent of the term
     * 
     * @var integer
     * @access protected
     */
    
var $_exp;
    
    
/**
     * Construct a polynomial term object
     * 
     * @access public
     * 
     * @param float $flt_coef
     * @param integer $int_exp
     */
    
function Math_PolynomialTerm($flt_coef 0.0$int_exp 0)
    {
        
$this->_coef = (float) $flt_coef;
        
$this->_exp = (int) $int_exp;
    }
    
    
/**
     * Get the exponent from the Polynomial Term
     * 
     * @access public
     * 
     * @return integer
     */
    
function getExponent()
    {
        return 
$this->_exp;
    }
    
    
/**
     * Get the coefficient from the Polynomial Term
     * 
     * @access public
     * 
     * @return float
     */
    
function getCoefficient()
    {
        return 
$this->_coef;
    }
    
    
/**
     * Set the exponent of the term
     * 
     * @access public
     * 
     * @param integer $int_exp
     * @return void
     */
    
function setExponent($int_exp)
    {
        
$this->_exp = (int) $int_exp;
    }
    
    
/**
     * Set the coefficient of the term
     * 
     * @access public
     * 
     * @param float $flt_coef
     * @return void
     */
    
function setCoefficient($flt_coef)
    {
        
$this->_coef = (float) $flt_coef;
    }
    
    
/**
     * Get a string representation of just this term
     * 
     * @access public
     * 
     * @return string
     */
    
function toString()
    {
        if (
$this->_coef != and $this->_exp 1) {
            return 
$this->_coef 'x^' $this->_exp;
        } else if (
$this->_coef != and $this->_exp == 1) {
            return 
$this->_coef 'x';
        } else if (
$this->_coef != 0) {
            return 
$this->_coef;
        } else {
            return 
'0';
        }
    }
}

?>

Examination of this class leads to the expection that the other classes in this package will regularly invoke the Math_PolynomialTerm class to represent and manipulate the individual terms of a polynomial expression.

The task of representing and manipulating complete polynomial expressions is handled by the Math_Polynomial and Math_PolynomialOp classes respectively. The Math_Polynomial class requires the Math_PolynomialTerm class but not the Math_PolynomialOp class. The Math_Polynomial class contains the package constructor, handles expression parsing and output duties, and contains some introspective capabilities. Below is an abridged version of the Math_Polynomial class that includes the constructor code but removes 1) some documentation and 2) implementation code for other methods (to make it easier to obtain an overview of the class functionality):

<?php

/**
 * Requires PolynomialTerm class to represent individual terms of the Polynomial
 */
require_once 'Math/Polynomial/PolynomialTerm.php';

/**
 * Require PEAR for PEAR errors
 */
require_once 'PEAR.php';

/**
 * Methods having a Polynomial or mixed type as a parameter should be able
 * to take either a string representation of a polynomial or a Polynomial
 * object as a parameter. 
 *
 * String representation of the Polynomial object can be retrieved with the 
 * toString() method. 
 * 
 * @package Math_Polynomial
 * @author Keith Palmer <keith@AcademicKeys.com>
 */
class Math_Polynomial
{
    
/**
     * An array of PolynomialTerm objects
     * 
     * @var array
     * @access protected
     */
    
var $_terms
    
    
/**
     * Whether or not Polynomial may contain multiple terms of the same degree
     * 
     * @var bool
     * @access protected
     */
    
var $_needs_combining;
    
    
/**
     * Whether or nothe Polynomial terms list needs to be sorted 
     * 
     * @var bool
     * @access protected
     */
    
var $_needs_sorting;
    
    
/**
     * Polynomial Constructor
     * 
     * Constructs an empty ( equal to 0 )  Polynomial object OR constructs a 
     * Polynomial from another Polynomial object or from a string representation
     * of a polynomial. 
     * 
     * @param mixed $mixed_poly String or Polynomial to construct from 
     */
    
function Math_Polynomial($mixed_poly null)
    {
        
// Init vars
        
$this->_terms = array();
        
$this->_needs_combining false;
        
$this->_needs_sorting false;
        
        
// Load from another Polynomail
        
if (is_a($mixed_poly'Math_Polynomial')) { 
            
$this->_terms $mixed_poly->_terms;
            
$this->_needs_combining $mixed_poly->_needs_combining;
            
$this->_needs_sorting $mixed_poly->_needs_sorting;
        } else { 
// Parse from string/integer
            
$this->_parsePolynomial(trim((string) $mixed_poly));
        }
    }
    
    
/**
     * Tell the number of terms in the Polynomial object
     * 
     * @access public
     * 
     * @see Math_Polynomial::getTerm(), Math_Polynomial::addTerm(), Math_PolynomialTerm
     * 
     * @return integer
     */
    
function numTerms() {}
    
    
/**
     * Get the term in the nth position from the Polynomial
     * 
     * @access public
     * 
     * @see Math_Polynomial::numTerms(), Math_Polynomial::addTerm(), Math_PolynomialTerm
     * 
     * @param integer $n
     * @return object
     */
    
function getTerm($n) {}
    
    
/**
     * Add a term to the Polynomial
     * 
     * @access public
     * 
     * @see Math_Polynomial::getTerm(), Math_Polynomial::numTerms(), Math_PolynomialTerm
     * 
     * @param object $term
     * @return bool
     */
    
function addTerm($term) {}
    
    
/**
     * Parse an individual term of Math_Polynomial into a Math_PolynomialTerm object
     * 
     * @access protected
     * 
     * @param string $str_term
     * 
     * @return object A Math_PolynomialTerm object
     */
    
function _parseTerm($str_term) {}
    
    
/**
     * Parse a string
     * 
     * @uses Math_Polynomial::_parseTerm()
     * 
     * @access protected
     * 
     * @param string $str_poly A string representation of a Polynomial
     * 
     * @return void
     */
    
function _parsePolynomial($str_poly) {}
    
    
/**
     * Retrieve a string representation of the Polynomial
     * 
     * @access public
     * 
     * @param bool $spaces Whether or not the string should contain spaces
     *                         ( i.e.: 4x^2 - 2x vs. 4x^2-2x )
     * 
     * @return string String representation
     */
    
function toString($spaces true) {}
    
    
/**
     * Sort the terms list
     * 
     * Uses a bubble-sort to sort the list of terms in descending order by their 
     * exponent values. 
     * 
     * @access protected
     * 
     * @return void
     */
    
function _sortTerms() {}
    
    
/**
     * Combine like terms inside the Polynomial object
     * 
     * Examines each term of the Polynomial and, if any of the terms have 
     * equivalent exponents, adds the two terms together ( add coefficients, 
     * keep the same exponent ) This is used to simplify the Polynomial for 
     * output. 
     * 
     * @access protected
     * 
     * @return void
     */
    
function _combineLikeTerms() {}
    
    
/**
     * Retrieve the degree ( highest exponent value ) of the polynomial
     * 
     * @access public
     * 
     * @return integer The degree of the polynomial
     */
    
function degree() {}
    
    
/**
     * Retrieve a string naming the degree of the polynomial
     * 
     * Tries to retrieve a string to name the polynomial degree, for example:
     * <samp>
     * 3 - Contant 
     * 3x - Linear equation
     * 3x^2 - Quadtratic
     * 3x^3 - Cubic
     * etc.
     * </samp>
     * 
     * @access public 
     * 
     * @return string String naming the polynomial degree
     */
    
function degreeString() {}

}

?>

What's Next

We will focus on the Math_PolynomialOp class (among other things) in the next blog which contains a large set of methods for manipulating polynomials. In other words, we will focus more on usage of the Math_Polynomial package in the next blog.

It is a happy situation when your math package provides a clear and logical mapping between a math domain and software objects. Keith Palmer has provided a nice software factorization of the domain of polynomials.

I encourage you to follow up this blog by reading Sedgewick & Wayne's discussion of Polynomials and reviewing their Polynomial.java class. You might consider whether a single Polynomial.java class or a set of classes (i.e., Math_Polynomial package) is the best way to represent polynomails and the algebraic operations that can be performed on them.

Permalink 

Differentiation with PHP (Part 5): Root Finding with Newton Raphson [Calculus
Posted on September 14, 2007 @ 01:27:07 PM by Paul Meagher

Solving a single variable equation means finding the x value that generates a 0 result in your equation. The general way to depict this idea is:

f(x) = 0

Geometrically, the value of x that generates a 0 result corresponds to a location on the graph of the equation where the curve intersects the x-axis (i.e., y = 0). This value is often a quantity of interest in scientific, engineering and business problems.

Finding the x value that generates a 0 result is also called "finding the root" of an equation. In this blog we will mostly be concerned with finding roots to single variable equations using the Newton Raphson method.

Analytic solutions for finding roots

What are the roots of this equation:

0 = x2 - 3

One way to tackle this problem is via root-finding formulas for this particular type of formula. This formula can be classified as a quadratic equation:

0 = ax2 + bx + c

We can use the quadratic root formulas we learned in high school to find the positive and negative roots:

x+ = (-b + sqrt(b*b - 4*a*c)) / 2*a
x_ = (-b - sqrt(b*b - 4*a*c)) / 2*a

Where a=1, b=1, and c = -3.

Here is a simple script that uses these quadratic root formulas to compute the roots of the equation 0 = x2 - 3.

<?php

// Script to find roots of the quadratic equation:
// 0 = x^2 - 3

$a 1;
$b 0;
$c = -3;

$radical sqrt($b*$b 4*$a*$c);

$x_pos = (-$b $radical) / 2*$a;
$x_neg = (-$b $radical) / 2*$a;

echo 
"x pos = $x_pos <br />";
echo 
"x neg = $x_neg <br />";

// Answers:
// x pos = 1.7320508075689
// x neg = -1.7320508075689 

?>

A more general approach to root finding involves graphing your equation and then using an iterative numerical root finding method such as the Newton Raphson method to approximate the numerical value of the roots (often to a very high degree of accuracy).

Mathematics of the Newton Raphson procedure

The mathematics involved in the Newton Raphson procedure is relatively simple. We begin by making an educated guess at the x value that generates a 0 result (i.e., a root value) by looking at a graph of our equation and noting where the curve crosses the x axis. We then feed our initial guess into the following iterative formula:

xn+1 = xn - f(xn)/f'(xn)

Where f'(xn) denotes the derivative of f(xn). The fact that we have to compute a derivative value on each iteration makes this a calculus-type problem.

The Newton Raphson iterative formula itself is simple a rearrangement of the formula for estimating a derivative when one of your y values (rise value) is 0:

f'(xn) = rise/run = f(xn) - 0 / xn - xn+1

It is left as an exercise for you to derive the Newton Raphson iterative formula from this derivative formula.

Usability of the Newton Raphson procedure

There are different ways that we might approach implementing the Newton Raphson procedure. One common method involves manually supplying the derivative formula to be used when computing the derivative. This adds a bit of extra work in deriving the derivative formula and passing it to the procedure. The tack that I decided to take was to have a NewtonRaphson class derive the derivative formula using Etienne Kneuss' Math_Derivative and Math_Expression packages discussed in previous blogs in this series. This means that we can invoke the Newton Raphson procedure by simply specifying the formula we want to find the root for along with an initial guess as to where that root might be located.

It should be noted that the Newton Raphson procedure differs from other root finding methods in that you don't have to specify the interval where the root is to be found, just an x value (denoted x0) close to the root x value. Compared to other root finding methods this makes the Newton Raphson procedure one of the easist to use, especially when our NewtonRaphson class will take care of deriving the derivative formula for us.

The NewtonRaphson class

The NewtonRaphson class takes is shown below. It consists of a few getters and a constructor where all the action takes place. The constructor finds the root and stores it in a root instance variable when it is found.

<?php

require_once 'Math/Expression.php';
require_once 
'Math/Derivative.php'

class 
NewtonRaphson extends Math_Derivative {
    
  
// maximum number of iterations allowed
  
public $max  100;
  
  
// tolerance used to terminate search
  
public $tol  0.00001;    

  
// storage for x value where root is located
  
public $root;  
  
  
// storage for x approximations
  
public $x;

  
// storage for y approximations
  
public $y;

  
// storage for dy approximations
  
public $dy;
    
  function 
getTolerance() {
    return 
$this->tol;
  }

  function 
getMaxIterations() {
    return 
$this->max;
  }

  function 
getRoot() {
      return 
$this->root;
  }

  function 
NewtonRaphson($fx$x0) {      
  
    
// reset arrays
    
$x        = array();
    
$y        = array();
    
$dy       = array();        

    
// clear the root 
    
$this->root  "";

    
// initial educated guess
    
$this->x[0]  = $x0;

    
// get function value at x0
    
$f_literal   str_replace('x'$x0$fx);
    
$express     = new Math_Expression($f_literal);
    
$this->y[0]  = $express->evaluate()->getValue(); 
    
    
// get derivative value at x0    
    
$fxdx        $this->getDerivative($fx'x');  
    
$df_literal  str_replace('x'$x0$fxdx);
    
$express     = new Math_Expression($df_literal);    
    
$this->dy[0] = $express->evaluate()->getValue();    
        
    for(
$i 0$i $this->max$i++) {

      
// approximate solution                
      
$this->x[$i+1]  = $this->x[$i] - ($this->y[$i] / $this->dy[$i]);
      
      
// get function value at x[i+1]
      
$f_literal      str_replace('x'$this->x[$i+1], $fx);
      
$express        = new Math_Expression($f_literal);
      
$this->y[$i+1]  = $express->evaluate()->getValue();      

      
// solution is good enough
      
if(abs($this->y[$i+1]) < $this->tol
        break;

      
// get derivative value at x[i+1]
      
$fxdx           $this->getDerivative($fx'x');         
      
$express        = new Math_Expression($df_literal);
      
$this->dy[$i+1] = $express->evaluate()->getValue();             

    }

    
$this->root array_pop($this->x);          

  } 
  
}

?>

Note that all intermediate computations are stored in the $this->x, $this->y, and $this->dy arrays. This gives you the option of studying the route the class took to converge on the root. The Newton Raphson method converges quickly on the root. Convergence is at least quadratic "which intuitively means that the number of correct digits roughly at least doubles in every step" (source Wikipedia).

Demonstrating usage

We will demonstrate usage by finding the roots of this equation:

0 = x2 - 3

Using analytic methods we established the roots to be:

x pos = 1.7320508075689
x neg = -1.7320508075689 

We will now use the NewtonRaphson class to see if we can approximate these roots using a more general iterative method of solution.

The first step in demonstrating usage of the NewtonRaphson class involves obtaining an overview of the function behavior by graphing it. I used JpGraph to graph the function and its derivative function:

Based on this graph, I offerred a guess that the positive root was located at x=1.2 and the negative root was located at -1.2. It is obvious that these aren't the exact values for the root but the Newton Raphson procedure should find them.

Here is the script I used to find the positive and negative roots.

<?php

include "NewtonRaphson.php";

$fx "x^2 - 3";   

$x0 1.2;
$nr = new NewtonRaphson($fx$x0);  

echo 
"The positive root is ".$nr->root."<br />";

$x0 = -1.2;
$nr = new NewtonRaphson($fx$x0);  

echo 
"The negative root is ".$nr->root;
    
?>

If you point your browser at this newton_raphson.php example script you should see the following output:

The positive root is 1.7320522345858
The negative root is -1.7320522345858

This is very close to the answer we obtained using analytic methods. We could have achieved a closer match by modifying the tolerance setting ($tol = 0.00001) in the NewtonRaphson class to a tinier value.

Updates to the Math_Expression package

In the course of developing and testing the NewtonRaphson class I discovered a couple of bugs in the Math_Expression package. It was not properly representing reals with exponent values and it was not handling unary minus operations correctly. I addressed the former issue by suggesting a new regex patterns to use for matching reals and Etienne addressed the unary minus issue. The result is that the Math_Expression package is more robust and usable than before and is available as a new version to be downloaded. I'd like to thank Etienne Kneuss for his assistance in addressing these bugs with me.

It would be great if other members of the opensource PHP community offered their assistance in further developing the Math_Expression package. The Math_Expression package offers an excellent foundation upon which to build a full-featured expression parser and evaluator along the lines of JEP - Java Math Expression Parser. Working on a highly useful project like this would pay rich dividends to php community and provide you with an excellent context in which to learn, practice and demonstrate compiler design skills.

Permalink 

Differentiation with PHP (Part 4): Math_Expression Package [Calculus
Posted on September 8, 2007 @ 12:21:03 PM by Paul Meagher

After installing Etienne Kneuss's Math_Derivative and Math_Expression packages under the "php/Math" system folder, I experimented with using the Math_Expression package to evaluate results generated by the Math_Derivative package. I had the suspicion that the two packages were built to work together and I was not disappointed in that suspicion.

According to Etienne, the "The goal [of the Math_Expression package] is to provide an easy way to have internal representations of mathematical objects, so any specific math package could be focused on the operations, leaving the parsing/output process to Math_Expression". Usage of the Math_Expression package is illustrated in the script below where the Math_Expression package is used to evaluate the results generated by the Math_Derivative package.

<?php

require_once 'Math/Expression.php';
require_once 
'Math/Derivative.php';   

$fx "3*x^3 + 4*x - 9"
$x 5;

$df  = new Math_Derivative();  

$fxdx $df->getDerivative($fx'x');   

$literal str_replace('x'$x$fxdx);

$ex = new Math_Expression($literal);

$result $ex->evaluate();

echo 
"fx = $fx <br />";
echo 
"(fx)dx = ($fx)dx = $fxdx <br />";   
echo 
"x = $x <br />";
echo 
"result = $result";
    
?>

Point your web browser at this script and this is what you should see:

fx = 3*x^3 + 4*x - 9
(fx)dx = (3*x^3 + 4*x - 9)dx = (3*1*3*x^2)+(4)
x = 5
result = 229

As Borat would say, "Very niiiiiice!!".

Literal Expressions

There is one step in the script above that you should pay close attention to - the step where you turn the symbolic expression (i.e., with abstract terms) into a literal formula (i.e., with constant terms) so that the Math_Expression package can evaluate the expression.

$literal = str_replace('x', $x, $fxdx); 

One wonders whether the Math_Expression package should perform this literalization step for you. Perhaps you could pass in an associative array of substitutions $subs to use when evaluating an expression:

$subs = array('x'=>5);
$ex = new Math_Expression($fxdx, $subs);
$result = $ex->evaluate();

Very Promising So Far

The prospects for performing fairly sophisticated differentiation with PHP are actually quite good given what I have shown you so far. There are lots of interesting applications that can be built on top of the Math_Expression and Math_Derivative packages. The extensibility of the Math_Expression library to handle vector quantities is also quite interesting and warrants further investigation in the context of more complex differentation tasks (e.g., partial differentiation).

We began this series by talking about a package that uses the Finite Difference Method (aka FDM) to estimate the first, second, third, and forth derivatives of a function. It would seem that the FDM is rendered useless by the powers of these symbolic differentiation packages to provide exact solutions. While this is true for basic differentation tasks, it turns out that the Finite Difference Method is quite useful in more complex tasks like partial differentation. So we can't write off Finite Difference Methods until we explore the requirements of more complex differentation tasks.

Math blogging

Writing about and programming calculus is a good way to learn calculus. It is also a good way to purge yourself of the calculus demon that will block your entry into the higher realms of math and science.

Not sure yet what I will write on next. The two calculus-related topics that I have spent some recent time thinking about/working on are 1) the Newton Raphson method of root finding, 2) the Math_Polynomial package, and 3) calculus graphing with JPGraph. There are many interesting calculus-topics to explore further so "Differentiation with PHP" may be a long ride :-)

I'd like to thank Etienne Kneuss for his technically and conceptually elegant packages and for reinvigorating this journey into differentiation with PHP. Keep up the good work!

Permalink 

Differentiation with PHP (Part 3): Math_Derivative Package [Calculus
Posted on August 31, 2007 @ 12:57:43 AM by Paul Meagher

The Finite Difference Method Package, aka the FDM Package, currently contains one subpackage: the Math_FDM_Derivative Package which estimates the derivative (i.e., rate of change) of a function at a particular x value by evaluating the function at neighboring points along the curve, computing differences among the coordinates, and dividing those differences to provide estimates of the first, second, third or fourth derivative. This finite difference method of computing a derivative can be contrasted with the symbolic differentiation approach used by the Math_Derivative Package developed by Etienne Kneuss. The Math_Derivative package accepts as input a string representing a function and returns the formula for the derivative of that function. To find the derivative of a curve at a particular point you simply plug an x value into this symbolically derived derivative formula.

The Math_Derivative package consists of a single Derivative.php class. The heart of the class is the parse method which consists of two main parts: 1) a code black to parse the supplied formula into a data structure that can be symbolically manipulated, and 2) a code block that applies the releavant differentiation rules to the parsed formula to compute it's derivative.

<?php

/**
  * Parses the expression and recursively calculates its derivative
  *
  *
  * @param string $expression expressionyou want to parse
  *
  * @return string the derivative
  *
  * @access protected
  */     
 
function parse($expression
 {
     
     
$expression $this->cleanExpression($expression);

     
$initial_expression $expression;
     
     
// if the expression doesn't rely on dx -> 0
     
if (!$this->reliesOndx($expression)) {
         return 
'0';
     }

     
// checks if it already exists in the cache
     
if ($this->_useCache && isset($this->_cache[$initial_expression])) {
         return 
$this->_cache[$initial_expression];
     }
     
     
// begins parsing
     
$depth 0;

     
// array containing the positions of the main operators
     
$main_operators = array('precedence' => 0'positions' => array());
     
     for (
$i=0$i strlen($expression); $i++) {
         
$char $expression[$i];

         if (
$char === '(') {
             
$depth++;
             continue;
         }
         
         if (
$char === ')') {
             
$depth--;
             continue;
         } 

         if (
$depth) {
             continue;
         }
         
         
         if (isset(
$this->_operatorsPrecedences[$char])) {
             
// current char is an operator
             
if ($main_operators['precedence'] == $this->_operatorsPrecedences[$char]) {
                 
// same type that the main operators : add its position to the list.
                 
$main_operators['positions'][] = array($char$i);
             } else if (
$main_operators['precedence'] > $this->_operatorsPrecedences[$char] || !$main_operators['precedence']) {
                 
// lower precedence than the list : this operator becomes a main operator
                 
$main_operators = array('precedence' => $this->_operatorsPrecedences[$char], 
                                         
'positions' => array(array($char$i)));
             }

         }
     }
     unset(
$depth);
     
     
// splits the expression using operators' positions
     
$pos 0;
     
$expression_parts = array();
     foreach (
$main_operators['positions'] as $operator) {
         
$expression_parts[] = substr($expression$pos$operator[1]-$pos);
         
$expression_parts[] = $operator[0];
         
$pos $operator[1]+1;
     }
     
$expression_parts[] = substr($expression$pos);
     
     unset(
$pos);
     
     
// dispatchs to the rule corresponding to the main operator
     
switch ($main_operators['precedence']) {

         case 
// + -
             
$expression $this->ruleAddition($expression_parts);
             break;
             
         case 
// *
             
$expression $this->ruleMultiplication($expression_parts);
             break;
             
         case 
3// /
             
$expression $this->ruleDivision($expression_parts);
             break;

         case 
4// ^
             
$expression $this->rulePower($expression_parts);
             break;
         
         default: 
// term
             
$expression $this->ruleTerm($expression);
             

     }
     
     
// put in cache
     
if ($this->_useCache) {
         
$this->_cache[$initial_expression] = $expression;
     }

     return 
$expression;
 }
 
?>

Rules of Differentiation

To understand the scope of this class requires you to understand something about the rules of differentiation. In the table below we provide a mathematical definition of each symbolic rule that the Math_Derivative.php class provides.

Math_Derivative Method Mathematical Definition
$this->rulePower($expression_parts) d/dx xn = n xn-1
$this->ruleAddition($expression_parts) d/dx (u + v) = du/dx + dv/dx
$this->ruleMultiplication($expression_parts) d/dx (u * v) = u * dv/dx + v * du/dx
$this->ruleDivision($expression_parts) d/dx (u / v) = [ v * du/dx + u * dv/dx ] / v2

The Derivative.php class also offers symbolic differentiation of transcentental functions. The relevant transformations are defined in a $_registeredFunctions associative array:

<?php

/**
  * Contains derivative forms of functions
  * that could appear in expressions
  *
  * @var array
  * @access protected
  */
 
var $_registeredFunctions = array('sin'  => 'cos(arg)*d(arg)',
                                   
'cos'  => '-sin(arg)*d(arg)',
                                   
'tan'  => '1/cos(arg)^2*d(arg)',
                                   
'ln'   => 'd(arg)/(arg)',
                                   
'log'  => 'd(arg)/(arg)',
                                   
'e'    => 'd(arg)*e(arg)',
                                   
'sqrt' => 'd((arg)^(1/2))',
                                   
'acos' => '-1/((1-(arg)^2)^(1/2))',
                                   
'asin' => '1/((1-(arg)^2)^(1/2))',
                                   
'atan' => '1/(1+(arg)^2)',
                                   
'sinh' => 'cosh(arg)*d(arg)',
                                   
'cosh' => 'sinh(arg)*d(arg)',
                                   
'tanh' => '(sech(arg)^2)*d(arg)',
                                   
'coth' => '-(csch(arg)^2)*d(arg)',
                                   
'sech' => '-sech(arg)*tanh(arg)*d(arg)',
                                   
'csch' => '-csch(arg)*coth(arg)*d(arg)',
                                   );
 
?>

Until next time....

I've not had time to turn the Math_FDM_Derivative code base into a nice package. Perhaps as I further explore the topic of "Differentiation with PHP" I will have time to do so. There are many interesting additional topics to explore including perhaps a return to the Math_Derviative package when I have had more time to absorb it and use it.

Exercises

1. Use the Math_Derivative package to compute the second derivative of the following function: f(x) = 5x3 - 3x5.

2. Graph a formula and the corresponding derivative of the function.

Permalink 

Differentiation with PHP (Part 2): Estimating Derivatives [Calculus
Posted on August 24, 2007 @ 03:15:14 PM by Paul Meagher

In Part 2 of "Differentiation with PHP" we will demonstrate how to estimate the derivative of a function using the FDD Package which I introduced in Part 1 of this series.

The script below is called estimate_first_derivative.php. It illustrates how to:

  1. Specify the function object whose derivative you will be estimating.
  2. Specify the $derivative_order and $estimation_method to use when instantiating the Derivative class.
  3. Obtain the derivative estimate via the estimate method.

<?php

include_once "../Derivative.php";

class 
{
  function 
valueAt($x) {
    return -
0.1*pow($x,4) - 0.15*pow($x,3) - 0.5*pow($x,2) - 0.25*$x 1.2;
  }
}

// Instantiate function to be passed to the derivative estimation methods.
$f = new F;
// Specify where you want to estimate the deriviative at.
$x 0.5;  
// Specify the increment or step size you want to use.
$h 0.25;

// Estimate "first" order derivative using "forward" method.
$forward = new FDD_Derivative("first""forward");
echo 
"Forward Estimate: ".$forward->estimate($f$x$h)."<br />";

// Estimate "first" order derivative using "backward" method.
$backward = new FDD_Derivative("first""backward");
echo 
"Backward Estimate: ".$backward->estimate($f$x$h)."<br />";

// Estimate "first" order derivative using "centered" method.
$centered = new FDD_Derivative("first""centered");
echo 
"Centered Estimate: ".$centered->estimate($f$x$h);

?>

The output of this script is:

Forward Estimate: -0.859375
Backward Estimate: -0.878125
Centered Estimate: -0.9125

All of the classes for computing the first derivative (i.e., forward, backward, centered) look similiar. They only differ in the actual formula used to compute the first derivative. Here is an example of what the class for computing the centered first derivative looks like:

<?php

class FDD_Derivative_Centered_First extends FDD_Derivative {

  public function 
least_accurate($f$x$h) {
    return (
$f->valueAt($x+$h)-$f->valueAt($x-$h))/(2*$h);      
  }

  public function 
more_accurate($f$x$h) {       
    return (-
$f->valueAt($x+2*$h) + 8*$f->valueAt($x+$h) - 8*$f->valueAt($x-$h) + $f->valueAt($x-2*$h))/(12*$h);      
  }

}

?>

All the estimation formulas contain a least_accurate estimation method and a more_accurate estimation method. The more accurate formula includes an extra Taylor series term. I'm a bit ambivalent about using the terms "least accurate" and "more accurate" as you can probably find cases where the "least accurate" estimation method produces a more accurate result than the "most accurate" estimation method. Perhaps more neutral method names like taylor1 and taylor2 would be better.

Package Layout

The codebase for the FDD Package has the following structure:

FDD/Derivative.php

FDD/Backward/First.php
FDD/Backward/Second.php
FDD/Backward/Third.php
FDD/Backward/Fourth.php

FDD/Centered/First.php
FDD/Centered/Second.php
FDD/Centered/Third.php
FDD/Centered/Fourth.php

FDD/Forward/First.php
FDD/Forward/Second.php
FDD/Forward/Third.php
FDD/Forward/Fourth.php

FDD/examples/estimate_first_derivative.php
FDD/examples/estimate_second_derivative.php
FDD/examples/estimate_third_derivative.php
FDD/examples/estimate_forth_derivative.php

Next week

I will tar up this package next week sometime when I write another installment of this series. I want to discuss how the FDD Package relates to symbolic differentiation, polynomial curve fitting, and computing derivatives from tabulated data. You should understand these concepts if you want to properly explore the vast domain of Differentiation and where the FDD package, and other PEAR packages, fit in this landscape.

Permalink 

Differentiation with PHP (Part 1): Derivative Base Class [Calculus
Posted on August 20, 2007 @ 12:37:52 AM by Paul Meagher

I developed a Finite Divided Differences Package, aka the FDD Package, that contains 24 methods that can be used to estimate derivatives of a function at a supplied x value. The user invokes the particular derivative method they are looking for by specifying the derivative order, estimation method, and accuracy they want.

  • In the constructor method, the user can specify what order of derivative they want to estimate {$derivative_order = "first", "second", "third", "fourth"}.
  • In the constructor method, the user can specify what estimation method to use {$estimation_method = "forward", "centered", "backward"}.
  • In the estimate method, the user can specify how may Taylor series terms to use in the estimation formula {$accuracy = "least", "more"}.

The base class of the FDD Package package, the Derivative.php class, implements a factory design pattern that invokes the code for the requested estimation method:

<?php
/**
* The FDD Package uses a Finite-Divided-Difference approach to estimating
* a derivative. The package implements the 24 FDD formulas mentioned 
* in:
*
* Chapra, S.C., Canale, R.P. (1998). Numerical Methods For Engineers. 
* 3rd Edition. McGraw-Hill. Chapter 23. Numerical Differentiation. 
* pp. 630-631.
*
* @author Paul Meagher
* @version 0.2
* @modified June 23, 2007
*/
class FDD_Derivative {

  
/**
  * The order of the derivative to be estimated.  Options are 
  * first, second, third, or fourth.
  *
  * @var string $order
  */
  
private $order;

  
/**
  * Method used to estimate the derivative.  Options are forward, 
  * centered, or backward.
  *
  * @var string $method
  */
  
private $method;

  
/**  
  * Constructor sets the derivative order and estimation method
  * to be used by the factory method to instantiate the relevant 
  * derivative estimation class.  The estimation method argument 
  * defaults to the centered method.  
  * 
  * @param string $derivative_order 
  * @param string $estimation_method
  */
  
public function __construct($derivative_order$estimation_method="centered") {
      
$this->order  ucfirst(strtolower($derivative_order));
      
$this->method ucfirst(strtolower($estimation_method));      
  }
  
  
/**
  * Lightweight factory instantiates a derivative estimation class for: 
  *
  * 1. The order of the derivative you want to estimate: first, 
  *    second, third, or fourth.
  *
  * 2. The estimation method you want to use: forward, backward, 
  *    or centered.  
  *
  * return instantiated estimation class
  */ 
  
private function factory() {   
        
    
// FDD estimation classes are housed in three subfolders 
    // corresponding to the names of the estimation method.
    // Within these subfolders are classes for each order of the 
    // derivative.
    
include_once $this->method."/".$this->order.".php";
    
    
// Create string corresponding to the class name of the 
    // estimation method.
    
$class_name "FDD_Derivative_".$this->method."_".$this->order;

    
// Return instantiated class.
    
return new $class_name;
  
  } 

  
/**
  * Invokes the least accurate or a more accurate formula to approximate
  * the derivative at a supplied x value. Defaults to the more accurate
  * formula (i.e, formula with the most terms) although many calculus 
  * textbooks tend to use the least accurate estimation formula.
  *
  * @param object $f function object
  * @param float $x compute derivative at this point
  * @param float $h step size
  * @param integer $accuracy options include {least, more}
  * @return float estimate of derivative
  */   
  
public function estimate($f$x$h$accuracy="more") {        
    
$fdd =& FDD_Derivative::factory();
    if (
$accuracy=="more"
        return 
$fdd->more_accurate($f$x$h);      
      elseif (
$accuracy=="least")
        return 
$fdd->least_accurate($f$x$h);
      else
        die(
"Error: Unrecognized accuracy setting.");
  }  

  
/**
  * Interface method that all subclasses must implement.
  *
  * Estimates the derivative using the least accurate version
  * of the relevant formula (i.e., uses the minimum number of 
  * terms in the Taylor series expansion of the formula).
  *
  * @param object $f function object
  * @param float $x compute derivative at this point
  * @param float $h step size
  * @return float least accurate estimate of derivative  
  */   
  
public function least_accurate($f$x$h) {}

  
/**
  * Interface method that all subclasses must implement.
  *  
  * Estimates the derivative using a more accurate version
  * of the relevant formula (i.e., uses one more term in the 
  * Taylor series expansion of the formula).
  *
  * @param object $f function object
  * @param float $x compute derivative at this point
  * @param float $h step size
  * @return float more accurate estimate of derivative
  */   
  
public function more_accurate($f$x$h) {}

}
?>

In my next blog I will illustrate usage of the Derivative.php base class.

Permalink 

 Archive 
 


php/Math Project
© 2011. All rights reserved.