Roots of a function: Difference between revisions
(common lisp entry) |
(omit m4) |
||
| Line 711: | Line 711: | ||
In this case, the roots are exact; inexact results are marked by decimal points. |
In this case, the roots are exact; inexact results are marked by decimal points. |
||
{{omit from|M4}} |
|||
Revision as of 01:47, 5 September 2009
You are encouraged to solve this task according to the task description, using any language you may know.
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)=x3-3x2+2x.
Ada
<lang ada>with Ada.Text_Io; use Ada.Text_Io;
procedure Roots_Of_Function is
package Real_Io is new Ada.Text_Io.Float_Io(Long_Float);
use Real_Io;
function F(X : Long_Float) return Long_Float is
begin
return (X**3 - 3.0*X*X + 2.0*X);
end F;
Step : constant Long_Float := 1.0E-6;
Start : constant Long_Float := -1.0;
Stop : constant Long_Float := 3.0;
Value : Long_Float := F(Start);
Sign : Boolean := Value > 0.0;
X : Long_Float := Start + Step;
begin
if Value = 0.0 then
Put("Root found at ");
Put(Item => Start, Fore => 1, Aft => 6, Exp => 0);
New_Line;
end if;
while X <= Stop loop
Value := F(X);
if (Value > 0.0) /= Sign then
Put("Root found near ");
Put(Item => X, Fore => 1, Aft => 6, Exp => 0);
New_Line;
elsif Value = 0.0 then
Put("Root found at ");
Put(Item => X, Fore => 1, Aft => 6, Exp => 0);
New_Line;
end if;
Sign := Value > 0.0;
X := X + Step;
end loop;
end Roots_Of_Function;</lang>
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)
AutoHotkey
Poly(x) is a test function of one variable.
We search for its roots.
roots() searches for intervals within given limits,
shifted by a given “step”, where our function has different
signs at the endpoints. Having found such an interval, the root() function searches for
a value where our
function is 0, within a given tolerance.
It also sets ErrorLevel to info about the root found.
ahk
discussion
<lang autohotkey>MsgBox % roots("poly", -0.99, 2, 0.1, 1.0e-5)
MsgBox % roots("poly", -1, 3, 0.1, 1.0e-5)
roots(f,x1,x2,step,tol) { ; search for roots in intervals of length "step", within tolerance "tol"
x := x1, y := %f%(x), s := (y>0)-(y<0)
Loop % ceil((x2-x1)/step) {
x += step, y := %f%(x), t := (y>0)-(y<0)
If (s=0 || s!=t)
res .= root(f, x-step, x, tol) " [" ErrorLevel "]`n"
s := t
}
Sort res, UN ; remove duplicate endpoints
Return res
}
root(f,x1,x2,d) { ; find x in [x1,x2]: f(x)=0 within tolerance d, by bisection
If (!y1 := %f%(x1))
Return x1, ErrorLevel := "Exact"
If (!y2 := %f%(x2))
Return x2, ErrorLevel := "Exact"
If (y1*y2>0)
Return "", ErrorLevel := "Need different sign ends!"
Loop {
x := (x2+x1)/2, y := %f%(x)
If (y = 0 || x2-x1 < d)
Return x, ErrorLevel := y ? "Approximate" : "Exact"
If ((y>0) = (y1>0))
x1 := x, y1 := y
Else
x2 := x, y2 := y
}
}
poly(x) {
Return ((x-3)*x+2)*x
}</lang>
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>
Common Lisp
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>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
<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>
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)
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:
<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>
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
<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
<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>
Ruby
<lang ruby>def sign(x)
x > 0 ? 1 : x == 0 ? 0 : -1
end
def find_roots(f, start, stop, step)
prev_sign = sign(f.call(start))
start.step(stop, step) do |x|
value = f.call(x)
if value == 0
# hit a root
puts "Root found at #{x}"
elsif prev_sign == -sign(value)
# passed a root in this iteration
puts "Root found between #{x-step} and #{x}"
end
prev_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 between -0.000999999999999979 and 2.08166817117217e-17 Root found at 1.0 Root found at 2.0
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.