Roots of a function
Roots of a function
You are encouraged to solve this task according to the task description, using any language you may know.
You are encouraged to solve this task according to the task description, using any language you may know.
Create a program that finds an 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++
#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 ); } }
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) ; }
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
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 ); }