Roots of a function: Difference between revisions
Rename Perl 6 -> Raku, alphabetize, minor clean-up
m (→function coded in-line: changed wording in the REXX section header.) |
Thundergnat (talk | contribs) (Rename Perl 6 -> Raku, alphabetize, minor clean-up) |
||
Line 226:
Return ((x-3)*x+2)*x
}</lang>
=={{header|Axiom}}==
Using a polynomial solver:
Line 372 ⟶ 373:
One can also use the GNU Scientific Library to find roots of functions. Compile with <pre>gcc roots.c -lgsl -lcblas -o roots</pre>
=={{header|C
{{trans|C++}}
Line 709 ⟶ 566:
// end brents_fun
}
</lang>
=={{header|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>
===Brent's Method===
Brent's Method uses a combination of the bisection method, inverse quadratic interpolation, and the secant method to find roots. It has a guaranteed run time equal to that of the bisection method (which always converges in a known number of steps (log2[(upper_bound-lower_bound)/tolerance] steps to be precise ) unlike the other methods), but the algorithm uses the much faster inverse quadratic interpolation and secant method whenever possible. The algorithm is robust and commonly used in libraries with a roots() function built in.
The algorithm is coded as a function that returns a double value for the root. The function takes an input that requires the function being evaluated, the lower and upper bounds, the tolerance one is looking for before converging (i recommend 0.0001) and the maximum number of iterations before giving up on finding the root (the root will always be found if the root is bracketed and a sufficient number of iterations is allowed).
The implementation is taken from the pseudo code on the wikipedia page for Brent's Method found here: https://en.wikipedia.org/wiki/Brent%27s_method.
<lang cpp>#include <iostream>
#include <cmath>
#include <algorithm>
#include <functional>
double brents_fun(std::function<double (double)> f, double lower, double upper, double tol, unsigned int max_iter)
{
double a = lower;
double b = upper;
double fa = f(a); // calculated now to save function calls
double fb = f(b); // calculated now to save function calls
double fs = 0; // initialize
if (!(fa * fb < 0))
{
std::cout << "Signs of f(lower_bound) and f(upper_bound) must be opposites" << std::endl; // throws exception if root isn't bracketed
return -11;
}
if (std::abs(fa) < std::abs(b)) // if magnitude of f(lower_bound) is less than magnitude of f(upper_bound)
{
std::swap(a,b);
std::swap(fa,fb);
}
double c = a; // c now equals the largest magnitude of the lower and upper bounds
double fc = fa; // precompute function evalutation for point c by assigning it the same value as fa
bool mflag = true; // boolean flag used to evaluate if statement later on
double s = 0; // Our Root that will be returned
double d = 0; // Only used if mflag is unset (mflag == false)
for (unsigned int iter = 1; iter < max_iter; ++iter)
{
// stop if converged on root or error is less than tolerance
if (std::abs(b-a) < tol)
{
std::cout << "After " << iter << " iterations the root is: " << s << std::endl;
return s;
} // end if
if (fa != fc && fb != fc)
{
// use inverse quadratic interopolation
s = ( a * fb * fc / ((fa - fb) * (fa - fc)) )
+ ( b * fa * fc / ((fb - fa) * (fb - fc)) )
+ ( c * fa * fb / ((fc - fa) * (fc - fb)) );
}
else
{
// secant method
s = b - fb * (b - a) / (fb - fa);
}
// checks to see whether we can use the faster converging quadratic && secant methods or if we need to use bisection
if ( ( (s < (3 * a + b) * 0.25) || (s > b) ) ||
( mflag && (std::abs(s-b) >= (std::abs(b-c) * 0.5)) ) ||
( !mflag && (std::abs(s-b) >= (std::abs(c-d) * 0.5)) ) ||
( mflag && (std::abs(b-c) < tol) ) ||
( !mflag && (std::abs(c-d) < tol)) )
{
// bisection method
s = (a+b)*0.5;
mflag = true;
}
else
{
mflag = false;
}
fs = f(s); // calculate fs
d = c; // first time d is being used (wasnt used on first iteration because mflag was set)
c = b; // set c equal to upper bound
fc = fb; // set f(c) = f(b)
if ( fa * fs < 0) // fa and fs have opposite signs
{
b = s;
fb = fs; // set f(b) = f(s)
}
else
{
a = s;
fa = fs; // set f(a) = f(s)
}
if (std::abs(fa) < std::abs(fb)) // if magnitude of fa is less than magnitude of fb
{
std::swap(a,b); // swap a and b
std::swap(fa,fb); // make sure f(a) and f(b) are correct after swap
}
} // end for
std::cout<< "The solution does not converge or iterations are not sufficient" << std::endl;
} // end brents_fun
</lang>
Line 1,939 ⟶ 1,940:
Note these roots are exact solutions with floating-point calculation.
=={{header|Octave}}==
Line 2,008 ⟶ 1,981:
x = x + step;
endwhile</lang>
=={{header|Oforth}}==
<lang Oforth>: findRoots(f, a, b, st)
| x y lasty |
a f perform dup ->y ->lasty
a b st step: x [
x f perform -> y
y ==0 ifTrue: [ System.Out "Root found at " << x << cr ]
else: [ y lasty * sgn -1 == ifTrue: [ System.Out "Root near " << x << cr ] ]
y ->lasty
] ;
: f(x) x 3 pow x sq 3 * - x 2 * + ; </lang>
{{out}}
<pre>
findRoots(#f, -1, 3, 0.0001)
Root found at 0
Root found at 1
Root found at 2
findRoots(#f, -1.000001, 3, 0.0001)
Root near 9.90000000000713e-005
Root near 1.000099
Root near 2.000099
</pre>
=={{header|ooRexx}}==
Line 2,275 ⟶ 2,276:
$sign = ( $value > 0 );
}</lang>
=={{header|Phix}}==
Line 2,410 ⟶ 2,386:
end locate_root;
</lang>
=={{header|PureBasic}}==
{{trans|C++}}
Line 2,588 ⟶ 2,565:
The profiling shows that memoization reduces the number of function calls
in this example from 184 to 67 (50 calls for primary interval division and about 6 calls for each point refinement).
=={{header|Raku}}==
(formerly Perl 6)
Uses exact arithmetic.
<lang perl6>sub f(\x) { x³ - 3*x² + 2*x }
my $start = -1;
my $stop = 3;
my $step = 0.001;
for $start, * + $step ... $stop -> $x {
state $sign = 0;
given f($x) {
my $next = .sign;
when 0.0 {
say "Root found at $x";
}
when $sign and $next != $sign {
say "Root found near $x";
}
NEXT $sign = $next;
}
}</lang>
{{out}}
<pre>Root found at 0
Root found at 1
Root found at 2</pre>
=={{header|REXX}}==
Line 2,736 ⟶ 2,740:
find_roots(-1..3) { |x| x**3 - 3*x**2 + 2*x }</lang>
=={{header|Scala}}==
| |||