Roots of a function: Difference between revisions

Content added Content deleted

Inline

Revision as of 14:40, 30 March 2009

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³-3x²+2x.

ALGOL 68

Finding 3 roots using the secant method:

MODE DBL = LONG REAL;
FORMAT dbl = $g(-long real width, long real width-6, -2)$;

MODE XY = STRUCT(DBL x, y);
FORMAT xy root = $f(dbl)" ("b("Exactly", "Approximately")")"$;

MODE DBLOPT = UNION(DBL, VOID);
MODE XYRES = UNION(XY, VOID);

PROC find root = (PROC (DBL)DBL f, DBLOPT in x1, in x2, in x error, in y error)XYRES:(
  INT limit = ENTIER (long real width / log(2)); # worst case of a binary search) #
  DBL x1 := (in x1|(DBL x1):x1|-5.0), # if x1 is EMPTY then -5.0 #
      x2 := (in x2|(DBL x2):x2|+5.0), 
      x error := (in x error|(DBL x error):x error|small real), 
      y error := (in y error|(DBL y error):y error|small real);
  DBL y1 := f(x1), y2;
  DBL dx := x1 - x2, dy;

  IF y1 = 0 THEN
    XY(x1, y1) # we already have a solution! #
  ELSE
    FOR i WHILE
      y2 := f(x2);
      IF y2 = 0 THEN stop iteration FI;
      IF i = limit THEN value error FI;
      IF y1 = y2 THEN value error FI;
      dy := y1 - y2;
      dx := dx / dy * y2;
      x1 := x2; y1 := y2; # retain for next iteration #
      x2 -:= dx;
  # WHILE # ABS dx > x error AND ABS dy > y error DO 
      SKIP
    OD;
    stop iteration: 
      XY(x2, y2) EXIT
    value error:
      EMPTY
  FI
);

PROC f = (DBL x)DBL: x UP 3 - LONG 3.1 * x UP 2 + LONG 2.0 * x;

DBL first root, second root, third root;

XYRES first result = find root(f, LENG -1.0, LENG 3.0, EMPTY, EMPTY);
CASE first result IN
  (XY first result): (
    printf(($"1st root found at x = "f(xy root)l$, x OF first result, y OF first result=0));
    first root := x OF first result
  )
  OUT printf($"No first root found"l$); stop
ESAC;

XYRES second result = find root( (DBL x)DBL: f(x) / (x - first root), EMPTY, EMPTY, EMPTY, EMPTY);
CASE second result IN
  (XY second result): (
    printf(($"2nd root found at x = "f(xy root)l$, x OF second result, y OF second result=0));
    second root := x OF second result
  )
  OUT printf($"No second root found"l$); stop
ESAC;

XYRES third result = find root( (DBL x)DBL: f(x) / (x - first root) / ( x - second root ), EMPTY, EMPTY, EMPTY, EMPTY);
CASE third result IN
  (XY third result): (
    printf(($"3rd root found at x = "f(xy root)l$, x OF third result, y OF third result=0));
    third root := x OF third result
  )
  OUT printf($"No third root found"l$); stop
ESAC

Output:


1st root found at x =  9.1557112297752398099031e-1 (Approximately)
2nd root found at x =  2.1844288770224760190097e 0 (Approximately)
3rd root found at x =  0.0000000000000000000000e 0 (Exactly)

C++

<lang 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 ); } }</lang>

D

<lang 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) ;

}</lang>

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

Fortran

Works with: Fortran version 90 and later

<lang fortran> PROGRAM ROOTS_OF_A_FUNCTION

  IMPLICIT NONE

  INTEGER, PARAMETER :: dp = SELECTED_REAL_KIND(15)
  REAL(dp) :: f, e, x, step, value
  LOGICAL :: s 
   
  f(x) = x*x*x - 3.0_dp*x*x + 2.0_dp*x
   
  x = -1.0_dp ; step = 1.0e-6_dp ; e = 1.0e-9_dp
 
  s = (f(x) > 0)
  DO WHILE (x < 3.0)
    value = f(x)
    IF(ABS(value) < e) THEN
      WRITE(*,"(A,F12.9)") "Root found at x =", x
      s = .NOT. s
    ELSE IF ((value > 0) .NEQV. s) THEN
      WRITE(*,"(A,F12.9)") "Root found near x = ", x
      s = .NOT. s
    END IF
    x = x + step
  END DO
 
END PROGRAM ROOTS_OF_A_FUNCTION</lang>

The following approach uses the Secant Method[1] to numerically find one root. Which root is found will depend on the start values x1 and x2 and if these are far from a root this method may not converge. <lang fortran> INTEGER, PARAMETER :: dp = SELECTED_REAL_KIND(15)

 INTEGER :: i=1, limit=100
 REAL(dp) :: d, e, f, x, x1, x2
  
 f(x) = x*x*x - 3.0_dp*x*x + 2.0_dp*x

 x1 = -1.0_dp ; x2 = 3.0_dp ; e = 1.0e-15_dp
  
 DO 
   IF (i > limit) THEN
     WRITE(*,*) "Function not converging"
     EXIT
   END IF
   d = (x2 - x1) / (f(x2) - f(x1)) * f(x2)
   IF (ABS(d) < e) THEN
     WRITE(*,"(A,F18.15)") "Root found at x = ", x2
     EXIT    
   END IF
   x1 = x2
   x2 = x2 - d
   i = i + 1
 END DO</lang>

GNU Octave

If the equation is a polynomial, we can put the coefficients in a vector and use roots:

<lang matlab>a = [ 1, -3, 2, 0 ]; r = roots(a); % let's print it for i = 1:3

 n = polyval(a, r(i));
 printf("x%d = %f (%f", i, r(i), n);
 if (n != 0.0)
   printf(" not");
 endif
 printf(" exact)\n");

endfor</lang>

Otherwise we can program our (simple) method:

Translation of: Python

<lang matlab>function y = f(x)

 y = x.^3 -3.*x.^2 + 2.*x;

endfunction

step = 0.001; tol = 10 .* eps; start = -1; stop = 3; se = sign(f(start));

x = start; while (x <= stop)

 v = f(x);
 if ( (v < tol) && (v > -tol) )
   printf("root at %f\n", x);
 elseif ( sign(v) != se )
   printf("root near %f\n", x);
 endif
 se = sign(v);
 x = x + step;

endwhile</lang>

J

J has builtin a root-finding operator, p., whose input is the (reversed) coeffiecients of the polynomial. Hence:

   1{::p.  0 2 _3 1    
2 1 0

We can determine whether the roots are exact or approximate by evaluating the polynomial at the candidate roots, and testing for zero:

   (0=]p.1{::p.) 0 2 _3 1 
1 1 1

As you can see, p. is also the operator which evaluates polynomials. This is not a coincidence.

Java

<lang java>public class Roots{ private static final double epsilon= 1E-10; //error bound, change for more or less accuracy

public double f(double x){ return x * x * x - 3 * x * x + 2 * x; //any formula you want here }

public void roots(double lowerBound, double upperBound, double step){ for(double x= lowerBound;x <= upperBound;x+= step){ double val; if(Math.abs(val= f(x)) < epsilon){ System.out.println(val); } } } }</lang>

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

<lang 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 );

}</lang>

Python

From Perl: <lang 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</lang>