Roots of a function

Create a program that finds and outputs the roots of a given function, range and (if applicable) step width. The program should identify whether the root is exact or approximate.

For this example, use f(x)=x^3-3x^2+2x.

C++

<cpp>#include <iostream>

double f(double x) { return (x*x*x - 3*x*x + 2*x); }

int main() { double step = 0.001; // Smaller step values produce more accurate and precise results double start = -1; double stop = 3; double value = f(start); double sign = (value > 0);

// Check for root at start if ( 0 == value ) std::cout << "Root found at " << start << std::endl;

for( double x = start + step; x <= stop; x += step ) { value = f(x);

if ( ( value > 0 ) != sign ) // We passed a root std::cout << "Root found near " << x << std::endl; else if ( 0 == value ) // We hit a root std::cout << "Root found at " << x << std::endl;

// Update our sign sign = ( value > 0 ); } }</cpp>

D

<d>module findroot ; import std.stdio ; import std.math ;

void report(T)(T[] r, T function(T) f, T tolerance = cast(T) 1e-4L) {

 if (r.length) {
   writefln("Root found (tolerance = %1.4g) :", tolerance) ;
   foreach(x ; r) {
     T y = f(x) ;    
     if (nearZero(y)) 
       writefln("... EXACTLY at %+1.20f, f(x) = %+1.4g", x, y) ;
     else if (nearZero(y, tolerance))
       writefln(".... MAY-BE at %+1.20f, f(x) = %+1.4g", x, y) ;
     else
       writefln("Verify needed, f(%1.4g) = %1.4g > tolerance in magnitude", x, y) ;
   }
 } else 
   writefln("No root found.") ;

}

bool nearZero(T)(T a, T b = T.epsilon * 4) { return abs(a) <= b ; }

T[] findroot(T)(T function(T) f, T start, T end, T step = cast(T) 0.001L,

                                T tolerance = cast(T) 1e-4L) {
 T[T] result ;

 if (nearZero(step))
   writefln("WARNING: step size may be too small.") ;

 T searchRoot(T a, T b) { // search root by simple bisection
   T root ;
   int limit = 49 ;
   T gap = b - a ;
   while (!nearZero(gap) && limit--) {
     if (nearZero(f(a))) return a ;
     if (nearZero(f(b))) return b ;
     root = (b + a)/2.0L ;
     if (nearZero(f(root))) return root ;
     if (f(a) * f(root) < 0)
       b = root ;
     else
       a = root ;
     gap = b - a ;           
   }
   return root ;
 }
 
 T dir = cast(T) (end > start ? 1.0 : -1.0) ;
 step = (end > start) ? abs(step) : - abs(step) ;
 for(T x = start ; x*dir <= end*dir ; x = x + step)
   if (f(x)*f(x + step) <= 0) {
     T r = searchRoot(x, x+ step) ;        
     result[r] = f(r) ;
   } 
 return result.keys.sort ; // reduce duplacated root, if any

}

real f(real x){

 return x*x*x - 3*x*x + 2*x ;

}

void main(){

   findroot(&f, -1.0L, 3.0L, 0.001L).report(&f) ;

}</d>

Output ( NB:smallest increment for real type in D is real.epsilon = 1.0842e-19 ):

Root found (tolerance = 0.0001) :
.... MAY-BE at -0.00000000000000000080, f(x) = -1.603e-18
... EXACTLY at +1.00000000000000000020, f(x) = -2.168e-19
.... MAY-BE at +1.99999999999999999950, f(x) = -8.674e-19

Maple

 f := x^3-3*x^2+2*x;
 roots(f,x);

outputs:

 [[0, 1], [1, 1], [2, 1]]

which means there are three roots. Each root is named as a pair where the first element is the value (0, 1, and 2), the second one the multiplicity (=1 for each means none of the three are degenerate).

By itself (i.e. unless specifically asked to do so), Maple will only perform exact (symbolic) operations and not attempt to do any kind of numerical approximation.

Mathematica

There are multiple obvious ways to do this in Mathematica.

Solve

This requires a full equation and will perform symbolic operations only:

In[1]:= Solve[x^3-3*x^2+2*x==0,x]
Out[1]= {{x->0},{x->1},{x->2}}

NSolve

This requires merely the polynomial and will perform numerical operations if needed:

In[2]:= NSolve[x^3 - 3*x^2 + 2*x , x]
Out[2]= {{x->0.},{x->1.},{x->2.}}

(note that the results here are floats)

FindRoot

This will numerically try to find one(!) local root from a given starting point:

In[3]:= FindRoot[x^3 - 3*x^2 + 2*x , {x, 1.5}]
Out[3]= {x->0.}
In[4]:= FindRoot[x^3 - 3*x^2 + 2*x , {x, 1.1}]
Out[4]= {x->1.}

(note that there is no guarantee which one is found).

FindInstance

This finds a value (optionally out of a given domain) for the given variable (or a set of values for a set of given variables) that satisfy a given equality or inequality:

In[5]:= FindInstance[x^3 - 3*x^2 + 2*x == 0, x]
Out[5]= {{x->0}}

Reduce

This will (symbolically) reduce a given expression to the simplest possible form, solving equations and performing substitutions in the process:

In[6]:= Reduce[x^3 - 3*x^2 + 2*x == 0, x]
Out[6]= x==0||x==1||x==2

(note that this doesn't yield a "solution" but a different expression that expresses the same thing as the original)

Perl

<perl>sub f {

       my $x = shift;

       return ($x * $x * $x - 3*$x*$x + 2*$x);

}

my $step = 0.001; # Smaller step values produce more accurate and precise results my $start = -1; my $stop = 3; my $value = &f($start); my $sign = $value > 0;

Check for root at start

print "Root found at $start\n" if ( 0 == $value );

for( my $x = $start + $step;

       $x <= $stop;
       $x += $step )

{

       $value = &f($x);

       if ( 0 == $value )
       {
               # We hit a root
               print "Root found at $x\n";
       }
       elsif ( ( $value > 0 ) != $sign )
       {
               # We passed a root
               print "Root found near $x\n";
       }

       # Update our sign
       $sign = ( $value > 0 );

}</perl>

Python

From Perl: <python>f = lambda x: x * x * x - 3 * x * x + 2 * x

step = 0.001 # Smaller step values produce more accurate and precise results start = -1 stop = 3

sign = f(start) > 0

x = start while x <= stop:

   value = f(x)

   if value == 0:
       # We hit a root
       print "Root found at", x
   elif (value > 0) != sign:
       # We passed a root
       print "Root found near", x

   # Update our sign
   sign = value > 0

   x += step</python>

Ada

C++

D