Roots of a function: Difference between revisions

Content added Content deleted

Inline

Revision as of 23:06, 25 March 2012

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

Works with: ALGOL 68 version Revision 1 - no extensions to language used

Works with: ALGOL 68G version Any - tested with release 1.18.0-9h.tiny

Finding 3 roots using the secant method: <lang algol68>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</lang> 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)

Axiom

Using a polynomial solver: <lang Axiom>expr := x^3-3*x^2+2*x solve(expr,x)</lang> Output: <lang Axiom> (1) [x= 2,x= 1,x= 0]

                         Type: List(Equation(Fraction(Polynomial(Integer))))</lang>

Using the secant method in the interpreter: <lang Axiom>digits(30) secant(eq: Equation Expression Float, binding: SegmentBinding(Float)):Float ==

 eps := 1.0e-30
 expr := lhs eq - rhs eq
 x := variable binding
 seg := segment binding
 x1 := lo seg
 x2 := hi seg
 fx1 := eval(expr, x=x1)::Float
 abs(fx1)<eps => return x1
 for i in 1..100 repeat
   fx2 := eval(expr, x=x2)::Float
   abs(fx2)<eps => return x2
   (x1, fx1, x2) := (x2, fx2, x2 - fx2 * (x2 - x1) / (fx2 - fx1))
 error "Function not converging."</lang>

The example can now be called using: <lang Axiom>secant(expr=0,x=-0.5..0.5)</lang>

BBC BASIC

<lang bbcbasic> function$ = "x^3-3*x^2+2*x"

     rangemin = -1
     rangemax = 3
     stepsize = 0.001
     accuracy = 1E-8
     PROCroots(function$, rangemin, rangemax, stepsize, accuracy)
     END
     
     DEF PROCroots(func$, min, max, inc, eps)
     LOCAL x, sign%, oldsign%
     oldsign% = 0
     FOR x = min TO max STEP inc
       sign% = SGN(EVAL(func$))
       IF sign% = 0 THEN
         PRINT "Root found at x = "; x
         sign% = -oldsign%
       ELSE IF sign% <> oldsign% AND oldsign% <> 0 THEN
           IF inc < eps THEN
             PRINT "Root found near x = "; x
           ELSE
             PROCroots(func$, x-inc, x+inc/8, inc/8, eps)
           ENDIF
         ENDIF
       ENDIF
       oldsign% = sign%
     NEXT x
     ENDPROC</lang>

Output:

Root found near x = 2.29204307E-9
Root found near x = 1
Root found at x = 2

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>

Clojure

Translation of: Haskell

(defn findRoots [f start stop step eps]

     (filter #(-> (f %) Math/abs (< eps)) (range start stop step)))

</lang>

> (findRoots #(+ (* % % %) (* -3 % %) (* 2 %)) -1.0 3.0 0.0001 0.00000001)
(-9.381755897326649E-14 0.9999999999998124 1.9999999999997022)

CoffeeScript

Translation of: Python

<lang coffeescript> print_roots = (f, begin, end, step) ->

 # Print approximate roots of f between x=begin and x=end,
 # using sign changes as an indicator that a root has been
 # encountered. 
 x = begin
 y = f(x)
 last_y = y
 
 cross_x_axis = ->
   (last_y < 0 and y > 0) or (last_y > 0 and y < 0)
 
 console.log '-----'
 while x <= end
   y = f(x)
   if y == 0
     console.log "Root found at", x
   else if cross_x_axis()
     console.log "Root found near", x
   x += step
   last_y = y

do ->

 # Smaller steps produce more accurate/precise results in general,
 # but for many functions we'll never get exact roots, either due
 # to imperfect binary representation or irrational roots.
 step = 1 / 256

 f1 = (x) -> x*x*x - 3*x*x + 2*x
 print_roots f1, -1, 5, step
 f2 = (x) -> x*x - 4*x + 3
 print_roots f2, -1, 5, step
 f3 = (x) -> x - 1.5
 print_roots f3, 0, 4, step
 f4 = (x) -> x*x - 2
 print_roots f4, -2, 2, step

</lang>

output

<lang> > coffee roots.coffee

Root found at 0 Root found at 1 Root found at 2

Root found at 1 Root found at 3

Root found at 1.5

Root found near -1.4140625 Root found near 1.41796875 </lang>

find-roots prints roots (and values near roots) and returns a list of root designators, each of which is either a number n, in which case (zerop (funcall function n)) is true, or a cons whose car and cdr are such that the sign of function at car and cdr changes.

<lang lisp>(defun find-roots (function start end &optional (step 0.0001))

 (let* ((roots '())
        (value (funcall function start))
        (plusp (plusp value)))
   (when (zerop value)
     (format t "~&Root found at ~W." start))
   (do ((x (+ start step) (+ x step)))
       ((> x end) (nreverse roots))
     (setf value (funcall function x))
     (cond
      ((zerop value)
       (format t "~&Root found at ~w." x)
       (push x roots))
      ((not (eql plusp (plusp value)))
       (format t "~&Root found near ~w." x)
       (push (cons (- x step) x) roots)))
     (setf plusp (plusp value)))))</lang>

> (find-roots #'(lambda (x) (+ (* x x x) (* -3 x x) (* 2 x))) -1 3)
Root found near 5.3588345E-5.
Root found near 1.0000072.
Root found near 2.000073.
((-4.6411653E-5 . 5.3588345E-5)
 (0.99990714 . 1.0000072)
 (1.9999729 . 2.000073))

D

<lang d>import std.stdio, std.math, std.algorithm;

bool nearZero(T)(in T a, in T b = T.epsilon * 4) pure nothrow {

   return abs(a) <= b;

}

T[] findRoot(T)(immutable T function(in T) pure nothrow f,

               in T start, in T end, in T step=cast(T)0.001L,
               T tolerance = cast(T)1e-4L) {
   if (nearZero(step))
       writefln("WARNING: step size may be too small.");
   immutable fi = f;

   /// Search root by simple bisection.
   T searchRoot(T a, T b) pure nothrow {
       T root;
       int limit = 49;
       T gap = b - a;

       while (!nearZero(gap) && limit--) {
           if (nearZero(fi(a))) return a;
           if (nearZero(fi(b))) return b;
           root = (b + a) / 2.0L;
           if (nearZero(fi(root))) return root;
           if (fi(a) * fi(root) < 0)
               b = root;
           else
               a = root;
           gap = b - a;
       }

       return root;
   }

   immutable dir = cast(T)(end > start ? 1.0 : -1.0);
   immutable step2 = (end > start) ? abs(step) : -abs(step);
   T[T] result;
   for (T x = start; (x * dir) <= (end * dir); x += step2)
       if (f(x) * f(x + step2) <= 0) {
           immutable T r = searchRoot(x, x + step2);
           result[r] = f(r);
       }

   // reduce duplacated root, if any
   return result.keys.sort().release();

}

void report(T)(in T[] r, immutable T function(in T) pure f,

              in T tolerance = cast(T)1e-4L) {
   if (r.length) {
       writefln("Root found (tolerance = %1.4g):", tolerance);
       foreach (x; r) {
           immutable 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.");

}

real f(in real x) pure nothrow {

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

}

void main() {

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

}</lang>

Output:

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

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

DWScript

Translation of: C

<lang delphi>type TFunc = function (x : Float) : Float;

function f(x : Float) : Float; begin

  Result := x*x*x-3.0*x*x +2.0*x;

end;

const e = 1.0e-12;

function Secant(xA, xB : Float; f : TFunc) : Float; const

  limit = 50;

var

  fA, fB : Float;
  d : Float;
  i : Integer;

begin

  fA := f(xA);
  for i := 0 to limit do begin
     fB := f(xB);
     d := (xB-xA)/(fB-fA)*fB;
     if Abs(d) < e then
        Exit(xB);
     xA := xB;
     fA := fB;
     xB -= d;
  end;
  PrintLn(Format('Function is not converging near (%7.4f,%7.4f).', [xA, xB]));
  Result := -99.0;

end;

const fstep = 1.0e-2;

var x := -1.032; // just so we use secant method var xx, value : Float; var s := f(x)>0.0;

while (x < 3.0) do begin

  value := f(x);
  if Abs(value)<e then begin
     PrintLn(Format("Root found at x= %12.9f", [x]));
     s := (f(x+0.0001)>0.0);
  end else if (value>0.0) <> s then begin
     xx := Secant(x-fstep, x, f);
     if xx <> -99.0 then   // -99 meaning secand method failed
        PrintLn(Format('Root found at x = %12.9f', [xx]))
     else PrintLn(Format('Root found near x= %7.4f', [xx]));
     s := (f(x+0.0001)>0.0);
  end;
  x += fstep;

end;</lang>

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 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>

Go

Secant method. No error checking. <lang go>package main

import (

   "fmt"
   "math"

)

func f(x float64) float64 {

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

}

var (

   lower float64 = -.5
   upper float64 = 2.6
   step float64 = 1

)

func main() {

   for x0, x1 := lower, lower+step; x0 < upper; x0, x1 = x1, x1+step {
       if x1 > upper {
           x1 = upper
       }
       r, status := secant(x0, x1)
       if status != "" && r >= x0 && r < x1 {
           fmt.Printf("  %6.3f %s\n", r, status)
       }
   }

}

func secant(x0, x1 float64) (float64, string) {

   var f0 float64
   f1 := f(x0)
   for i := 0; i < 100; i++ {
       f0, f1 = f1, f(x1)
       switch {
       case f1 == 0:
           return x1, "exact"
       case math.Abs(x1-x0) < 1e-6:
           return x1, "approximate"
       }
       x0, x1 = x1, x1-f1*(x1-x0)/(f1-f0)
   }
   return 0, ""

}</lang> Output:

   0.000 approximate
   1.000 exact
   2.000 approximate

Haskell

<lang haskell>f x = x^3-3*x^2+2*x

findRoots start stop step eps =

 [x | x <- [start, start+step .. stop], abs (f x) < eps]</lang>

Executed in GHCi: <lang haskell>*Main> findRoots (-1.0) 3.0 0.0001 0.000000001 [-9.381755897326649e-14,0.9999999999998124,1.9999999999997022]</lang>

Or using package hmatrix from HackageDB. <lang haskell>import Numeric.GSL.Polynomials import Data.Complex

Main> mapM_ print $ polySolve [0,2,-3,1]

(-5.421010862427522e-20) :+ 0.0 2.000000000000001 :+ 0.0 0.9999999999999996 :+ 0.0</lang> No complex roots, so: <lang haskell>*Main> mapM_ (print.realPart) $ polySolve [0,2,-3,1] -5.421010862427522e-20 2.000000000000001 0.9999999999999996</lang>

HicEst

HicEst's SOLVE function employs the Levenberg-Marquardt method: <lang HicEst>OPEN(FIle='test.txt')

1 DLG(NameEdit=x0, DNum=3)

x = x0 chi2 = SOLVE(NUL=x^3 - 3*x^2 + 2*x, Unknown=x, I=iterations, NumDiff=1E-15) EDIT(Text='approximate exact ', Word=(chi2 == 0), Parse=solution)

WRITE(FIle='test.txt', LENgth=6, Name) x0, x, solution, chi2, iterations GOTO 1</lang> <lang HicEst>x0=0.5; x=1; solution=exact; chi2=79E-32 iterations=65; x0=0.4; x=2E-162 solution=exact; chi2=0; iterations=1E4; x0=0.45; x=1; solution=exact; chi2=79E-32 iterations=67; x0=0.42; x=2E-162 solution=exact; chi2=0; iterations=1E4; x0=1.5; x=1.5; solution=approximate; chi2=0.1406; iterations=14: x0=1.54; x=1; solution=exact; chi2=44E-32 iterations=63; x0=1.55; x=2; solution=exact; chi2=79E-32 iterations=55; x0=1E10; x=2; solution=exact; chi2=18E-31 iterations=511; x0=-1E10; x=0; solution=exact; chi2=0; iterations=1E4;</lang>

J

J has builtin a root-finding operator, p., whose input is the coeffiecients of the polynomial (where the exponent of the indeterminate variable matches the index of the coefficient: 0 1 2 would be 0 + x + (2 times x squared)). Hence:

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

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

That said, we could also implement the technique used by most others here. Specifically: we can implement the function as a black box and check every 1 millionth of a unit between minus one and three, and we can test that result for exactness.

<lang J> blackbox=: 0 2 _3 1&p.

  (#~ (=<./)@:|@blackbox) i.&.(1e6&*)&.(1&+) 3

0 1 2

  0=blackbox 0 1 2

1 1 1</lang>

Here, we see that each of the results (0, 1 and 2) are as accurate as we expect our computer arithmetic to be. (The = returns 1 where paired values are equal and 0 where they are not equal).

Java

<lang java>public class Roots {

   public interface Function {

public double f(double x);

   private static int sign(double x) {

return (x < 0.0) ? -1 : (x > 0.0) ? 1 : 0;

   public static void printRoots(Function f, double lowerBound,

double upperBound, double step) { double x = lowerBound, ox = x; double y = f.f(x), oy = y; int s = sign(y), os = s;

for (; x <= upperBound ; x += step) { s = sign(y = f.f(x)); if (s == 0) { System.out.println(x); } else if (s != os) { double dx = x - ox; double dy = y - oy; double cx = x - dx * (y / dy); System.out.println("~" + cx); } ox = x; oy = y; os = s; }

   public static void main(String[] args) {

Function poly = new Function () { public double f(double x) { return x*x*x - 3*x*x + 2*x; } }; printRoots(poly, -1.0, 4, 0.002);

}</lang> Produces this output:

~2.616794878713638E-18
~1.0000000000000002
~2.000000000000001

Liberty BASIC

<lang lb>' Finds and output the roots of a given function f(x), ' within a range of x values.

' [RC]Roots of an function

   mainwin 80 12

   xMin  =-1
   xMax  = 3
   y     =f( xMin) '   Since Liberty BASIC has an 'eval(' function the fn
   '                       and limits would be better entered via 'input'.
   LastY =y

   eps  =1E-12 '   closeness acceptable

   bigH=0.01

   print
   print " Checking for roots of x^3 -3 *x^2 +2 *x =0 over range -1 to +3"
   print

   x=xMin: dx = bigH
   do
       x=x+dx
       y = f(x)
       'print x, dx, y
       if y*LastY <0 then 'there is a root, should drill deeper
           if dx < eps then    'we are close enough
               print " Just crossed axis, solution f( x) ="; y; " at x ="; using( "#.#####", x)
               LastY = y
               dx = bigH   'after closing on root, continue with big step
           else
               x=x-dx  'step back
               dx = dx/10  'repeat with smaller step
           end if
       end if
   loop while x<xMax

   print
   print " Finished checking in range specified."

end

   function f( x)
       f =x^3 -3 *x^2 +2 *x
   end function</lang>

Maple

<lang maple>f := x^3-3*x^2+2*x; roots(f,x);</lang>

outputs:

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: <lang mathematica>Solve[x^3-3*x^2+2*x==0,x]</lang> Output

 {{x->0},{x->1},{x->2}}

NSolve

This requires merely the polynomial and will perform numerical operations if needed: <lang mathematica> NSolve[x^3 - 3*x^2 + 2*x , x]</lang> Output

 {{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: <lang mathematica>FindRoot[x^3 - 3*x^2 + 2*x , {x, 1.5}]</lang> Output

 {x->0.}

From a different start point: <lang mathematica>FindRoot[x^3 - 3*x^2 + 2*x , {x, 1.1}]</lang> Output

{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: <lang mathematica> FindInstance[x^3 - 3*x^2 + 2*x == 0, x]</lang> Output

{{x->0}}

Reduce

This will (symbolically) reduce a given expression to the simplest possible form, solving equations and performing substitutions in the process: <lang mathematica>Reduce[x^3 - 3*x^2 + 2*x == 0, x]</lang>

 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)

Objeck

Translation of: C++

<lang objeck> bundle Default {

 class Roots {
   function : f(x : Float) ~ Float
   {
     return (x*x*x - 3.0*x*x + 2.0*x);
   }
    
   function : Main(args : String[]) ~ Nil 
   {
     step := 0.001;
     start := -1.0;
     stop := 3.0;
     value := f(start);
     sign := (value > 0);
    
     if(0.0 = value) {
       start->PrintLine();
     };
     
     for(x := start + step; x <= stop;  x += step;)  {
       value := f(x);
       
       if((value > 0) <> sign) {
         IO.Console->Instance()->Print("~")->PrintLine(x);
       }
       else if(0 = value) {
         IO.Console->Instance()->Print("~")->PrintLine(x);
       };
       
       sign := (value > 0);
     };
   }
 }

} </lang>

OCaml

A general root finder using the False Position (Regula Falsi) method, which will find all simple roots given a small step size.

<lang ocaml>let bracket u v =

 ((u > 0.0) && (v < 0.0)) || ((u < 0.0) && (v > 0.0));;

let xtol a b = (a = b);; (* or use |a-b| < epsilon *)

let rec regula_falsi a b fa fb f =

 if xtol a b then (a, fa) else
 let c = (fb*.a -. fa*.b) /. (fb -. fa) in
 let fc = f c in
 if fc = 0.0 then (c, fc) else
 if bracket fa fc then
   regula_falsi a c fa fc f
 else
   regula_falsi c b fc fb f;;

let search lo hi step f =

 let rec next x fx =
   if x > hi then [] else
     let y = x +. step in
     let fy = f y in
     if fx = 0.0 then
       (x,fx) :: next y fy
     else if bracket fx fy then
       (regula_falsi x y fx fy f) :: next y fy
     else
       next y fy in
 next lo (f lo);;

let showroot (x,fx) =

 Printf.printf "f(%.17f) = %.17f [%s]\n" 
   x fx (if fx = 0.0 then "exact" else "approx") in

let f x = ((x -. 3.0)*.x +. 2.0)*.x in List.iter showroot (search (-5.0) 5.0 0.1 f);;</lang>

Output:

f(0.00000000000000000) = 0.00000000000000000 [exact]
f(1.00000000000000022) = 0.00000000000000000 [exact]
f(1.99999999999999978) = 0.00000000000000000 [exact]

Note these roots are exact solutions with floating-point calculation.

Octave

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

<lang octave>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 octave>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>

PARI/GP

Gourdon–Schönhage algorithm

<lang parigp>polroots(x^3-3*x^2+2*x)</lang>

Newton's method

This uses a modified version of the Newton–Raphson method. <lang parigp>polroots(x^3-3*x^2+2*x,1)</lang>

Brent's method

<lang parigp>solve(x=-.5,.5,x^3-3*x^2+2*x) solve(x=.5,1.5,x^3-3*x^2+2*x) solve(x=1.5,2.5,x^3-3*x^2+2*x)</lang>

Factorization to linear factors

<lang parigp>findRoots(P)={

 my(f=factor(P),t);
 for(i=1,#f[,1],
   if(poldegree(f[i,1]) == 1,
     for(j=1,f[i,2],
       print(-polcoeff(f[i,1], 0), " (exact)")
     )
   );
   if(poldegree(f[i,1]) > 1,
     t=polroots(f[i,1]);
     for(j=1,#t,
       for(k=1,f[i,2],
         print(if(imag(t[j]) == 0.,real(t[j]),t[j]), " (approximate)")
       )
     )
   )
 )

}; findRoots(x^3-3*x^2+2*x)</lang>

Factorization to quadratic factors

Of course this process could be continued to degrees 3 and 4 with sufficient additional work. <lang parigp>findRoots(P)={

 my(f=factor(P),t);
 for(i=1,#f[,1],
   if(poldegree(f[i,1]) == 1,
     for(j=1,f[i,2],
       print(-polcoeff(f[i,1], 0), " (exact)")
     )
   );
   if(poldegree(f[i,1]) == 2,
     t=solveQuadratic(polcoeff(f[i,1],2),polcoeff(f[i,1],1),polcoeff(f[i,1],0));
     for(j=1,f[i,2],
       print(t[1]" (exact)\n"t[2]" (exact)")
     )
   );
   if(poldegree(f[i,1]) > 2,
     t=polroots(f[i,1]);
     for(j=1,#t,
       for(k=1,f[i,2],
         print(if(imag(t[j]) == 0.,real(t[j]),t[j]), " (approximate)")
       )
     )
   )
 )

}; solveQuadratic(a,b,c)={

 my(t=-b/2/a,s=b^2/4/a^2-c/a,inner=core(numerator(s))/core(denominator(s)),outer=sqrtint(s/inner));
 if(inner < 0,
   outer *= I;
   inner *= -1
 );
 s=if(inner == 1,
   outer
 ,
   if(outer == 1,
     Str("sqrt(", inner, ")")
   ,
     Str(outer, " * sqrt(", inner, ")")
   )
 );
 if (t,
   [Str(t, " + ", s), Str(t, " - ", s)]
 ,
   [s, Str("-", s)]
 )

}; findRoots(x^3-3*x^2+2*x)</lang>

Pascal

Translation of: Fortran

<lang pascal>Program RootsFunction;

var

 e, x, step, value: double;
 s: boolean;
 i, limit: integer;
 x1, x2, d: double;

function f(const x: double): double;

 begin
   f := x*x*x - 3*x*x + 2*x;
 end;

begin

 x    := -1;
 step := 1.0e-6;
 e    := 1.0e-9;
 s    := (f(x) > 0);

 writeln('Version 1: simply stepping x:');
 while x < 3.0 do
 begin
   value := f(x);
   if abs(value) < e then
   begin
     writeln ('root found at x = ', x);
     s := not s;
   end
   else if ((value > 0) <> s) then
   begin
     writeln ('root found at x = ', x);
     s := not s;
   end;
   x := x + step;
 end;
 
 writeln('Version 2: secant method:');
 x1 := -1.0;
 x2 :=  3.0;
 e  :=  1.0e-15;
 i  :=  1;
 limit := 300;
 while true do
 begin
   if i > limit then
   begin
     writeln('Error: function not converging');
     exit;
   end;
   d := (x2 - x1) / (f(x2) - f(x1)) * f(x2);
   if abs(d) < e then
   begin
     if d = 0 then
       write('Exact ')
     else
       write('Approximate ');
     writeln('root found at x = ', x2);
     exit;
   end;
   x1 := x2;
   x2 := x2 - d;
   i  := i + 1;
 end;

end. </lang> Output:

Version 1: simply stepping x:
root found at x =  7.91830063542152E-012
root found at x =  1.00000000001584E+000
root found at x =  1.99999999993357E+000
Version 2: secant method:
Exact root found at x =  1.00000000000000E+000

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>

PicoLisp

Translation of: Clojure

<lang PicoLisp>(de findRoots (F Start Stop Step Eps)

  (filter
     '((N) (> Eps (abs (F N))))
     (range Start Stop Step) ) )

(scl 12)

(mapcar round

  (findRoots
     '((X) (+ (*/ X X X `(* 1.0 1.0)) (*/ -3 X X 1.0) (* 2 X)))
     -1.0 3.0 0.0001 0.00000001 ) )</lang>

Output:

-> ("0.000" "1.000" "2.000")

PL/I

<lang PL/I> f: procedure (x) returns (float (18));

  declare x float (18);
  return (x**3 - 3*x**2 + 2*x );

end f;

declare eps float, (x, y) float (18); declare dx fixed decimal (15,13);

eps = 1e-12;

do dx = -5.03 to 5 by 0.1;

  x = dx;
  if sign(f(x)) ^= sign(f(dx+0.1)) then
     call locate_root;

end;

locate_root: procedure;

  declare (left, mid, right) float (18);

  put skip list ('Looking for root in [' || x, x+0.1 || ']' );
  left = x; right = dx+0.1;
  PUT SKIP LIST (F(LEFT), F(RIGHT) );
  if abs(f(left) ) < eps then
     do; put skip list ('Found a root at x=', left); return; end;
  else if abs(f(right) ) < eps then
     do; put skip list ('Found a root at x=', right); return; end;
  do forever;
     mid = (left+right)/2;
     if sign(f(mid)) = 0 then
        do; put skip list ('Root found at x=', mid); return; end;
     else if sign(f(left)) ^= sign(f(mid)) then
        right = mid;
     else
        left = mid;
     /* put skip list (left || right); */
     if abs(right-left) < eps then
        do; put skip list ('There is a root near ' ||
           (left+right)/2); return;
        end;
  end;

end locate_root; </lang>

PureBasic

Translation of: C++

<lang PureBasic>Procedure.d f(x.d)

 ProcedureReturn x*x*x-3*x*x+2*x

EndProcedure

Procedure main()

 OpenConsole()
 Define.d StepSize= 0.001
 Define.d Start=-1, stop=3
 Define.d value=f(start), x=start
 Define.i oldsign=Sign(value)
 
 If value=0
   PrintN("Root found at "+StrF(start))
 EndIf
 
 While x<=stop
   value=f(x)
   If Sign(value) <> oldsign
     PrintN("Root found near "+StrF(x))
   ElseIf value = 0
     PrintN("Root found at "+StrF(x))
   EndIf
   oldsign=Sign(value)
   x+StepSize
 Wend

EndProcedure

main()</lang>

Python

Translation of: 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>

R

Translation of: Octave

<lang R>f <- function(x) x^3 -3*x^2 + 2*x

findroots <- function(f, begin, end, tol = 1e-20, step = 0.001) {

 se <- ifelse(sign(f(begin))==0, 1, sign(f(begin)))
 x <- begin
 while ( x <= end ) {
   v <- f(x)
   if ( abs(v) < tol ) {
     print(sprintf("root at %f", x))
   } else if ( ifelse(sign(v)==0, 1, sign(v)) != se )  {
     print(sprintf("root near %f", x))
   }
   se <- ifelse( sign(v) == 0 , 1, sign(v))
   x <- x + step
 }

}

findroots(f, -1, 3)</lang>

RLaB

RLaB implements a number of solvers from the GSL and the netlib that find the roots of a real or vector function of a real or vector variable. The solvers are grouped with respect whether the variable is a scalar, findroot, or a vector, findroots. Furthermore, for each group there are two types of solvers, one that does not require the derivative of the objective function (which root(s) are being sought), and one that does.

The script that finds a root of a scalar function $f(x)=x^{3}-3\,x^{2}+2\,x$ of a scalar variable x using the bisection method on the interval -5 to 5 is, <lang RLaB> f = function(x) {

 rval = x .^ 3 - 3 * x .^ 2 + 2 * x;
 return rval;

};

>> findroot(f, , [-5,5])

</lang>

For a detailed description of the solver and its parameters interested reader is directed to the rlabplus manual.

Ruby

Translation of: Python

<lang ruby>def sign(x)

 x <=> 0

end

def find_roots(f, range, step)

 sign = sign(f[range.begin])
 
 range.step(step) do |x|
   value = f[x]

   if value == 0
     puts "Root found at #{x}"
   elsif sign(value) == -sign
     puts "Root found between #{x-step} and #{x}"
   end

   sign = sign(value)
 end

end

f = lambda { |x| x**3 - 3*x**2 + 2*x } find_roots(f, -1..3, 0.001)</lang>

Output:

Root found at 0.0
Root found at 1.0
Root found at 2.0

Or we could use Enumerable#inject, monkey patching and block:

<lang ruby>class Numeric

 def sign
   zero? ? 0 : self / abs
 end

end

def find_roots(range, step = 1e-3)

 range.step( step ).inject( yield(range.begin).sign ) { |sign, x|
   value = yield(x)
   if value == 0
     puts "Root found at #{x}"
   elsif value.sign == -sign
     puts "Root found between #{x-step} and #{x}"
   end
   value.sign
 }

end

find_roots(-1..3) { |x| x**3 - 3*x**2 + 2*x }</lang>

Tcl

This simple brute force iteration marks all results, with a leading "~", as approximate. This version always reports its results as approximate because of the general limits of computation using fixed-width floating-point numbers (i.e., IEEE double-precision floats). <lang Tcl>proc froots {lambda {start -3} {end 3} {step 0.0001}} {

   set res {}
   set lastsign [sgn [apply $lambda $start]]
   for {set x $start} {$x <= $end} {set x [expr {$x + $step}]} {
       set sign [sgn [apply $lambda $x]]
       if {$sign != $lastsign} {
           lappend res [format ~%.11f $x]
       }
       set lastsign $sign
   }
   return $res

} proc sgn x {expr {($x>0) - ($x<0)}}

puts [froots {x {expr {$x**3 - 3*$x**2 + 2*$x}}}]</lang> Result and timing:

/Tcl $ time ./froots.tcl
~0.00000000000 ~1.00000000000 ~2.00000000000

real    0m0.368s
user    0m0.062s
sys     0m0.030s

A more elegant solution (and faster, because you can usually make the initial search coarser) is to use brute-force iteration and then refine with Newton-Raphson, but that requires the differential of the function with respect to the search variable. <lang Tcl>proc frootsNR {f df {start -3} {end 3} {step 0.001}} {

   set res {}
   set lastsign [sgn [apply $f $start]]
   for {set x $start} {$x <= $end} {set x [expr {$x + $step}]} {
       set sign [sgn [apply $f $x]]
       if {$sign != $lastsign} {
           lappend res [format ~%.15f [nr $x $f $df]]
       }
       set lastsign $sign
   }
   return $res

} proc sgn x {expr {($x>0) - ($x<0)}} proc nr {x1 f df} {

   # Newton's method converges very rapidly indeed
   for {set iters 0} {$iters < 10} {incr iters} {
       set x1 [expr {
           [set x0 $x1] - [apply $f $x0]/[apply $df $x0]
       }]
       if {$x0 == $x1} {
           break
       }
   }
   return $x1

}

puts [frootsNR \

   {x {expr {$x**3 - 3*$x**2 + 2*$x}}} \
   {x {expr {3*$x**2 - 6*$x + 2}}}]</lang>

TI-89 BASIC

Finding roots is a built-in function: zeros(x^3-3x^2+2x, x) returns {0,1,2}.

In this case, the roots are exact; inexact results are marked by decimal points.