Roots of a function: Difference between revisions
Added Easylang
Alextretyak (talk | contribs) (Added 11l) |
(Added Easylang) |
||
(17 intermediate revisions by 14 users not shown) | |||
Line 13:
{{trans|Python}}
<
R x
-V step = 0.001
Line 32:
sgn = value > 0
x += step</
{{out}}
Line 42:
=={{header|Ada}}==
<
procedure Roots_Of_Function is
Line 80:
X := X + Step;
end loop;
end Roots_Of_Function;</
=={{header|ALGOL 68}}==
Line 88:
{{wont work with|ELLA ALGOL 68|Any (with appropriate job cards) - tested with release [http://sourceforge.net/projects/algol68/files/algol68toc/algol68toc-1.8.8d/algol68toc-1.8-8d.fc9.i386.rpm/download 1.8-8d] - due to extensive use of FORMATted transput}}
Finding 3 roots using the secant method:
<
FORMAT dbl = $g(-long real width, long real width-6, -2)$;
Line 157:
)
OUT printf($"No third root found"l$); stop
ESAC</
Output:
<pre>1st root found at x = 9.1557112297752398099031e-1 (Approximately)
Line 163:
3rd root found at x = 0.0000000000000000000000e 0 (Exactly)
</pre>
=={{header|Arturo}}==
<syntaxhighlight lang="arturo">f: function [n]->
((n^3) - 3*n^2) + 2*n
step: 0.01
start: neg 1.0
stop: 3.0
sign: positive? f start
x: start
while [x =< stop][
value: f x
if? value = 0 ->
print ["root found at" to :string .format:".5f" x]
else ->
if sign <> value > 0 -> print ["root found near" to :string .format:".5f" x]
sign: value > 0
'x + step
]</syntaxhighlight>
{{out}}
<pre>root found near 0.00000
root found near 1.00000
root found near 2.00000</pre>
=={{header|ATS}}==
<syntaxhighlight lang="ats">
#include
"share/atspre_staload.hats"
Line 212 ⟶ 241:
main0 () =
findRoots (~1.0, 3.0, 0.001, lam (x) => x*x*x - 3.0*x*x + 2.0*x, 0, 0.0)
</syntaxhighlight>
=={{header|AutoHotkey}}==
Line 221 ⟶ 250:
[http://www.autohotkey.com/forum/viewtopic.php?t=44657&postdays=0&postorder=asc&start=139 discussion]
<
MsgBox % roots("poly", -1, 3, 0.1, 1.0e-5)
Line 256 ⟶ 285:
poly(x) {
Return ((x-3)*x+2)*x
}</
=={{header|Axiom}}==
Using a polynomial solver:
<
solve(expr,x)</
Output:
<
Type: List(Equation(Fraction(Polynomial(Integer))))</
Using the secant method in the interpreter:
<
secant(eq: Equation Expression Float, binding: SegmentBinding(Float)):Float ==
eps := 1.0e-30
Line 280 ⟶ 309:
abs(fx2)<eps => return x2
(x1, fx1, x2) := (x2, fx2, x2 - fx2 * (x2 - x1) / (fx2 - fx1))
error "Function not converging."</
The example can now be called using:
<
=={{header|BBC BASIC}}==
<
rangemin = -1
rangemax = 3
Line 311 ⟶ 340:
oldsign% = sign%
NEXT x
ENDPROC</
Output:
<pre>Root found near x = 2.29204307E-9
Line 321 ⟶ 350:
=== Secant Method ===
<
#include <stdio.h>
Line 380 ⟶ 409:
}
return 0;
}</
=== GNU Scientific Library ===
<
#include <stdio.h>
Line 400 ⟶ 429:
return 0;
}</
One can also use the GNU Scientific Library to find roots of functions. Compile with <pre>gcc roots.c -lgsl -lcblas -o roots</pre>
Line 408 ⟶ 437:
{{trans|C++}}
<
class Program
Line 441 ⟶ 470:
}
}
}</
{{trans|Java}}
<
class Program
Line 486 ⟶ 515:
PrintRoots(f, -1.0, 4, 0.002);
}
}</
===Brent's Method===
{{trans|C++}}
<
class Program
Line 597 ⟶ 626:
// end brents_fun
}
</syntaxhighlight>
=={{header|C++}}==
<
double f(double x)
Line 635 ⟶ 664:
sign = ( value > 0 );
}
}</
===Brent's Method===
Line 643 ⟶ 672:
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.
<
#include <cmath>
#include <algorithm>
Line 741 ⟶ 770:
} // end brents_fun
</syntaxhighlight>
=={{header|Clojure}}==
{{trans|Haskell}}
<syntaxhighlight lang="clojure">
(defn findRoots [f start stop step eps]
(filter #(-> (f %) Math/abs (< eps)) (range start stop step)))
</syntaxhighlight>
<pre>
Line 759 ⟶ 788:
=={{header|CoffeeScript}}==
{{trans|Python}}
<
print_roots = (f, begin, end, step) ->
# Print approximate roots of f between x=begin and x=end,
Line 795 ⟶ 824:
f4 = (x) -> x*x - 2
print_roots f4, -2, 2, step
</syntaxhighlight>
output
<syntaxhighlight lang="text">
> coffee roots.coffee
-----
Line 813 ⟶ 842:
Root found near -1.4140625
Root found near 1.41796875
</syntaxhighlight>
=={{header|Common Lisp}}==
Line 821 ⟶ 850:
<code>find-roots</code> prints roots (and values near roots) and returns a list of root designators, each of which is either a number <code><var>n</var></code>, in which case <code>(zerop (funcall function <var>n</var>))</code> is true, or a <code>cons</code> whose <code>car</code> and <code>cdr</code> are such that the sign of function at car and cdr changes.
<
(let* ((roots '())
(value (funcall function start))
Line 837 ⟶ 866:
(format t "~&Root found near ~w." x)
(push (cons (- x step) x) roots)))
(setf plusp (plusp value)))))</
<pre>> (find-roots #'(lambda (x) (+ (* x x x) (* -3 x x) (* 2 x))) -1 3)
Line 848 ⟶ 877:
=={{header|D}}==
<
bool nearZero(T)(in T a, in T b = T.epsilon * 4) pure nothrow {
Line 919 ⟶ 948:
findRoot(&f, -1.0L, 3.0L, 0.001L).report(&f);
}</
{{out}}
<pre>Root found (tolerance = 0.0001):
Line 929 ⟶ 958:
=={{header|Dart}}==
{{trans|Scala}}
<
findRoots(Function(double) f, double start, double stop, double step, double epsilon) sync* {
Line 940 ⟶ 969:
// Vector(-9.381755897326649E-14, 0.9999999999998124, 1.9999999999997022)
print(findRoots(fn, -1.0, 3.0, 0.0001, 0.000000001));
}</
=={{header|Delphi}}==
See [https://rosettacode.org/wiki/Roots_of_a_function#Pascal Pascal].
=={{header|DWScript}}==
{{trans|C}}
<
function f(x : Float) : Float;
Line 994 ⟶ 1,026:
end;
x += fstep;
end;</
=={{header|EasyLang}}==
<syntaxhighlight>
func f x .
return x * x * x - 3 * x * x + 2 * x
.
numfmt 6 0
proc findroot start stop step . .
x = start
while x <= stop
val = f x
if val = 0
print x & " (exact)"
elif sign val <> sign0 and sign0 <> 0
print x & " (err = " & step & ")"
.
sign0 = sign val
x += step
.
.
proc drawfunc start stop . .
linewidth 0.3
drawgrid
x = start
while x <= stop
line x * 10 + 50 f x * 10 + 50
x += 0.1
.
.
drawfunc -1 3
findroot -1 3 pow 2 -20
print ""
findroot -1 3 1e-6
</syntaxhighlight>
=={{header|EchoLisp}}==
We use the 'math' library, and define f(x) as the polynomial : x<sup>3</sup> -3x<sup>2</sup> +2x
<
(lib 'math.lib)
Lib: math.lib loaded.
Line 1,011 ⟶ 1,077:
(root f -1000 (- 2 epsilon)) → 1.385559938161431e-7 ;; 0
(root f epsilon (- 2 epsilon)) → 1.0000000002190812 ;; 1
</syntaxhighlight>
=={{header|Elixir}}==
{{trans|Ruby}}
<
def find_roots(f, range, step \\ 0.001) do
first .. last = range
Line 1,041 ⟶ 1,107:
f = fn x -> x*x*x - 3*x*x + 2*x end
RC.find_roots(f, -1..3)</
{{out}}
Line 1,051 ⟶ 1,117:
=={{header|Erlang}}==
<
-module(roots).
-export([main/0]).
Line 1,079 ⟶ 1,145:
while(X+Step, Step, Start, Stop, Value > 0,F)
end.
</syntaxhighlight>
{{out}}
<pre>Root found near 8.81239525796218e-16
Line 1,087 ⟶ 1,153:
=={{header|ERRE}}==
<syntaxhighlight lang="erre">
PROGRAM ROOTS_FUNCTION
Line 1,144 ⟶ 1,210:
END LOOP
END PROGRAM
</syntaxhighlight>
Note: Outputs are calculated in single precision.
{{out}}
Line 1,159 ⟶ 1,225:
=={{header|Fortran}}==
{{works with|Fortran|90 and later}}
<
IMPLICIT NONE
Line 1,184 ⟶ 1,250:
END DO
END PROGRAM ROOTS_OF_A_FUNCTION</
The following approach uses the [[wp:Secant_method|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.
<
INTEGER :: i=1, limit=100
REAL(dp) :: d, e, f, x, x1, x2
Line 1,207 ⟶ 1,273:
x2 = x2 - d
i = i + 1
END DO</
=={{header|FreeBASIC}}==
Simple bisection method.
<syntaxhighlight lang="freebasic">#Include "crt.bi"
const iterations=20000000
sub bisect( f1 as function(as double) as double,min as double,max as double,byref O as double,a() as double)
dim as double last,st=(max-min)/iterations,v
for n as double=min to max step st
v=f1(n)
if sgn(v)<>sgn(last) then
redim preserve a(1 to ubound(a)+1)
a(ubound(a))=n
O=n+st:exit sub
end if
last=v
next
end sub
function roots(f1 as function(as double) as double,min as double,max as double, a() as double) as long
redim a(0)
dim as double last,O,st=(max-min)/iterations,v
for n as double=min to max step st
v=f1(n)
if sgn(v)<>sgn(last) and n>min then bisect(f1,n-st,n,O,a()):n=O
last=v
next
return ubound(a)
end function
Function CRound(Byval x As Double,Byval precision As Integer=30) As String
If precision>30 Then precision=30
Dim As zstring * 40 z:Var s="%." &str(Abs(precision)) &"f"
sprintf(z,s,x)
If Val(z) Then Return Rtrim(Rtrim(z,"0"),".")Else Return "0"
End Function
function defn(x as double) as double
return x^3-3*x^2+2*x
end function
redim as double r()
print
if roots(@defn,-20,20,r()) then
print "in range -20 to 20"
print "All roots approximate"
print "number","root to 6 dec places","function value at root"
for n as long=1 to ubound(r)
print n,CRound(r(n),6),,defn(r(n))
next n
end if
sleep</syntaxhighlight>
{{out}}
<pre>in range -20 to 20
All roots approximate
number root to 6 dec places function value at root
1 0 -2.929925652002424e-009
2 1 1.477781779325033e-009
3 2 -2.897852187377925e-009</pre>
=={{header|Go}}==
Secant method. No error checking.
<
import (
Line 1,247 ⟶ 1,373:
}
return 0, ""
}</
Output:
<pre>
Line 1,256 ⟶ 1,382:
=={{header|Haskell}}==
<
findRoots start stop step eps =
[x | x <- [start, start+step .. stop], abs (f x) < eps]</
Executed in GHCi:
<
[-9.381755897326649e-14,0.9999999999998124,1.9999999999997022]</
Or using package [http://hackage.haskell.org/package/hmatrix hmatrix] from HackageDB.
<
import Data.Complex
Line 1,271 ⟶ 1,397:
(-5.421010862427522e-20) :+ 0.0
2.000000000000001 :+ 0.0
0.9999999999999996 :+ 0.0</
No complex roots, so:
<
-5.421010862427522e-20
2.000000000000001
0.9999999999999996</
It is possible to solve the problem directly and elegantly using robust bisection method and Alternative type class.
<
data Root a = Exact a | Approximate a deriving (Show, Eq)
Line 1,299 ⟶ 1,425:
findRoots f [] = empty
findRoots f [x] = if f x == 0 then pure (Exact x) else empty
findRoots f (a:b:xs) = bisection f a b <|> findRoots f (b:xs)</
It is possible to use these functions with different Alternative functors: IO, Maybe or List:
Line 1,321 ⟶ 1,447:
=={{header|HicEst}}==
HicEst's [http://www.HicEst.com/SOLVE.htm SOLVE] function employs the Levenberg-Marquardt method:
<
1 DLG(NameEdit=x0, DNum=3)
Line 1,330 ⟶ 1,456:
WRITE(FIle='test.txt', LENgth=6, Name) x0, x, solution, chi2, iterations
GOTO 1</
<
x0=0.4; x=2E-162 solution=exact; chi2=0; iterations=1E4;
x0=0.45; x=1; solution=exact; chi2=79E-32 iterations=67;
Line 1,339 ⟶ 1,465:
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;</
=={{header|Icon}} and {{header|Unicon}}==
Line 1,345 ⟶ 1,471:
Works in both languages:
<
showRoots(f, -1.0, 4, 0.002)
end
Line 1,372 ⟶ 1,498:
procedure sign(x)
return (x<0, -1) | (x>0, 1) | 0
end</
Output:
Line 1,387 ⟶ 1,513:
J has builtin a root-finding operator, '''<tt>p.</tt>''', 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:
<
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:
<
1 1 1</
As you can see, <tt>p.</tt> is also the operator which evaluates polynomials. This is not a coincidence.
Line 1,399 ⟶ 1,525:
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.
<
(#~ (=<./)@:|@blackbox) i.&.(1e6&*)&.(1&+) 3
0 1 2
0=blackbox 0 1 2
1 1 1</
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).
=={{header|Java}}==
<
public interface Function {
public double f(double x);
Line 1,445 ⟶ 1,571:
printRoots(poly, -1.0, 4, 0.002);
}
}</
Produces this output:
<pre>~2.616794878713638E-18
Line 1,455 ⟶ 1,581:
{{works with|SpiderMonkey|22}}
{{works with|Firefox|22}}
<
// This function notation is sorta new, but useful here
// Part of the EcmaScript 6 Draft
Line 1,486 ⟶ 1,612:
printRoots(poly, -1.0, 4, 0.002);
</syntaxhighlight>
=={{header|jq}}==
Line 1,499 ⟶ 1,625:
printRoots/3 emits an array of results, each of which is either a
number (representing an exact root within the limits of machine arithmetic) or a string consisting of "~" followed by an approximation to the root.
<
if . < 0 then -1 elif . > 0 then 1 else 0 end;
Line 1,521 ⟶ 1,647:
end )
| .[3] ;
</syntaxhighlight>
We present two examples, one where step is a power of 1/2, and one where it is not:
{{Out}}
<
[
Line 1,537 ⟶ 1,663:
"~1.0000000000000002",
"~1.9999999999999993"
]</
=={{header|Julia}}==
Line 1,543 ⟶ 1,669:
Assuming that one has the Roots package installed:
<
println(find_zero(x -> x^3 - 3x^2 + 2x, (-100, 100)))</
{{out}}
Line 1,553 ⟶ 1,679:
Without the Roots package, Newton's method may be defined in this manner:
<
##f: the function of x
##fp: the derivative of f
Line 1,574 ⟶ 1,700:
end
end
</syntaxhighlight>
Finding the roots of f(x) = x3 - 3x2 + 2x:
<syntaxhighlight lang="julia">
f(x) = x^3 - 3*x^2 + 2*x
fp(x) = 3*x^2-6*x+2
x_s, count = newton(f,fp,1.00)
</syntaxhighlight>
{{out}}
Line 1,590 ⟶ 1,716:
=={{header|Kotlin}}==
{{trans|C}}
<
typealias DoubleToDouble = (Double) -> Double
Line 1,639 ⟶ 1,765:
x += step
}
}</
{{out}}
Line 1,649 ⟶ 1,775:
=={{header|Lambdatalk}}==
<
1) defining the function:
{def func {lambda {:x} {+ {* 1 :x :x :x} {* -3 :x :x} {* 2 :x}}}}
Line 1,680 ⟶ 1,806:
- a root found at 540°
- a root found at 720°
</syntaxhighlight>
=={{header|Liberty BASIC}}==
<
' within a range of x values.
Line 1,728 ⟶ 1,854:
function f( x)
f =x^3 -3 *x^2 +2 *x
end function</
=={{header|Lua}}==
<
function f (x) return x^3 - 3*x^2 + 2*x end
Line 1,760 ⟶ 1,886:
for _, r in pairs(root(f, -1, 3, 10^-6)) do
print(string.format("%0.12f", r.val), r.err)
end</
{{out}}
<pre>Root (to 12DP) Max. Error
Line 1,768 ⟶ 1,894:
2.000000999934 1e-06</pre>
Note that the roots found are all near misses because fractional numbers that seem nice and 'round' in decimal (such as 10^-6) often have some rounding error when represented in binary. To increase the chances of finding exact integer roots, try using an integer start value with a step value that is a power of two.
<
print("Root (to 12DP)\tMax. Error\n")
for _, r in pairs(root(f, -1, 3, 2^-10)) do
print(string.format("%0.12f", r.val), r.err)
end</
{{out}}
<pre>Root (to 12DP) Max. Error
Line 1,782 ⟶ 1,908:
=={{header|Maple}}==
<
roots(f,x);</
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).
Line 1,793 ⟶ 1,919:
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.
=={{header|Mathematica}}/{{header|Wolfram Language}}==
There are multiple obvious ways to do this in Mathematica.
===Solve===
This requires a full equation and will perform symbolic operations only:
<
Output
<pre> {{x->0},{x->1},{x->2}}</pre>
Line 1,805 ⟶ 1,929:
===NSolve===
This requires merely the polynomial and will perform numerical operations if needed:
<
Output
<pre> {{x->0.},{x->1.},{x->2.}}</pre>
Line 1,812 ⟶ 1,936:
===FindRoot===
This will numerically try to find one(!) local root from a given starting point:
<
Output
<pre> {x->0.}</pre>
From a different start point:
<
Output
<pre>{x->1.}</pre>
Line 1,823 ⟶ 1,947:
===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:
<
Output
<pre>{{x->0}}</pre>
Line 1,829 ⟶ 1,953:
===Reduce===
This will (symbolically) reduce a given expression to the simplest possible form, solving equations and performing substitutions in the process:
<
<pre> x==0||x==1||x==2</pre>
(note that this doesn't yield a "solution" but a different expression that expresses the same thing as the original)
=={{header|Maxima}}==
<
/* Number of roots in a real interval, using Sturm sequences */
Line 1,887 ⟶ 2,011:
[x=1.16154139999725193608791768724717407484314725802151429063617b0*%i + 3.41163901914009663684741869855524128445594290948999288901864b-1,
x=3.41163901914009663684741869855524128445594290948999288901864b-1 - 1.16154139999725193608791768724717407484314725802151429063617b0*%i,
x=-6.82327803828019327369483739711048256891188581897998577803729b-1]</
=={{header|Nim}}==
<
import strformat
Line 1,911 ⟶ 2,035:
sign = value > 0
x += step</
{{out}}
Line 1,922 ⟶ 2,046:
=={{header|Objeck}}==
{{trans|C++}}
<
bundle Default {
class Roots {
Line 1,957 ⟶ 2,081:
}
}
</syntaxhighlight>
=={{header|OCaml}}==
Line 1,963 ⟶ 2,087:
A general root finder using the False Position (Regula Falsi) method, which will find all simple roots given a small step size.
<
((u > 0.0) && (v < 0.0)) || ((u < 0.0) && (v > 0.0));;
Line 1,995 ⟶ 2,119:
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);;</
Output:
Line 2,010 ⟶ 2,134:
If the equation is a polynomial, we can put the coefficients in a vector and use ''roots'':
<
r = roots(a);
% let's print it
Line 2,020 ⟶ 2,144:
endif
printf(" exact)\n");
endfor</
Otherwise we can program our (simple) method:
{{trans|Python}}
<
y = x.^3 -3.*x.^2 + 2.*x;
endfunction
Line 2,045 ⟶ 2,169:
se = sign(v);
x = x + step;
endwhile</
=={{header|Oforth}}==
<
| x y lasty |
a f perform dup ->y ->lasty
Line 2,060 ⟶ 2,184:
] ;
: f(x) x 3 pow x sq 3 * - x 2 * + ; </
{{out}}
Line 2,076 ⟶ 2,200:
=={{header|ooRexx}}==
<
a*x**3+b*x**2+c*x+d =(x-x1)*(x-x2)*(x-x3)
*/
Line 2,130 ⟶ 2,254:
Else
Return xr
::requires 'rxmath' LIBRARY</
{{out}}
<pre>p=-3
Line 2,145 ⟶ 2,269:
=={{header|PARI/GP}}==
===Gourdon–Schönhage algorithm===<!-- X. Gourdon, "Algorithmique du théorème fondamental de l'algèbre" (1993). -->
<
===Newton's method===
This uses a modified version of the Newton–Raphson method.
<
===Brent's method===
<
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)</
===Factorization to linear factors===
<
my(f=factor(P),t);
for(i=1,#f[,1],
Line 2,175 ⟶ 2,299:
)
};
findRoots(x^3-3*x^2+2*x)</
===Factorization to quadratic factors===
Of course this process could be continued to degrees 3 and 4 with sufficient additional work.
<
my(f=factor(P),t);
for(i=1,#f[,1],
Line 2,224 ⟶ 2,348:
)
};
findRoots(x^3-3*x^2+2*x)</
=={{header|Pascal}}==
{{trans|Fortran}}
<
var
Line 2,292 ⟶ 2,416:
end;
end.
</syntaxhighlight>
Output:
<pre>
Line 2,304 ⟶ 2,428:
=={{header|Perl}}==
<
{
my $x = shift;
Line 2,340 ⟶ 2,464:
# Update our sign
$sign = ( $value > 0 );
}</
=={{header|Phix}}==
{{trans|CoffeeScript}}
<!--<syntaxhighlight lang="phix">(phixonline)-->
<span style="color: #008080;">procedure</span> <span style="color: #000000;">print_roots</span><span style="color: #0000FF;">(</span><span style="color: #004080;">integer</span> <span style="color: #000000;">f</span><span style="color: #0000FF;">,</span> <span style="color: #004080;">atom</span> <span style="color: #000000;">start</span><span style="color: #0000FF;">,</span> <span style="color: #000000;">stop</span><span style="color: #0000FF;">,</span> <span style="color: #000000;">step</span><span style="color: #0000FF;">)</span>
<span style="color: #000080;font-style:italic;">--
-- Print approximate roots of f between x=start and x=stop, using
-- sign changes as an indicator that a root has been encountered.
<span style="color: #004080;">atom</span> <span style="color: #000000;">x</span> <span style="color: #0000FF;">=</span> <span style="color: #000000;">start</span><span style="color: #0000FF;">,</span> <span style="color: #000000;">y</span> <span style="color: #0000FF;">=</span> <span style="color: #000000;">0</span>
<span style="color: #7060A8;">puts</span><span style="color: #0000FF;">(</span><span style="color: #000000;">1</span><span style="color: #0000FF;">,</span><span style="color: #008000;">"-----\n"</span><span style="color: #0000FF;">)</span>
<span style="color: #008080;">while</span> <span style="color: #000000;">x</span><span style="color: #0000FF;"><=</span><span style="color: #000000;">stop</span> <span style="color: #008080;">do</span>
<span style="color: #004080;">atom</span> <span style="color: #000000;">last_y</span> <span style="color: #0000FF;">=</span> <span style="color: #000000;">y</span>
<span style="color: #000000;">y</span> <span style="color: #0000FF;">=</span> <span style="color: #000000;">f</span><span style="color: #0000FF;">(</span><span style="color: #000000;">x</span><span style="color: #0000FF;">)</span>
<span style="color: #008080;">if</span> <span style="color: #000000;">y</span><span style="color: #0000FF;">=</span><span style="color: #000000;">0</span>
<span style="color: #008080;">or</span> <span style="color: #0000FF;">(</span><span style="color: #000000;">last_y</span><span style="color: #0000FF;"><</span><span style="color: #000000;">0</span> <span style="color: #008080;">and</span> <span style="color: #000000;">y</span><span style="color: #0000FF;">></span><span style="color: #000000;">0</span><span style="color: #0000FF;">)</span>
<span style="color: #008080;">or</span> <span style="color: #0000FF;">(</span><span style="color: #000000;">last_y</span><span style="color: #0000FF;">></span><span style="color: #000000;">0</span> <span style="color: #008080;">and</span> <span style="color: #000000;">y</span><span style="color: #0000FF;"><</span><span style="color: #000000;">0</span><span style="color: #0000FF;">)</span> <span style="color: #008080;">then</span>
<span style="color: #7060A8;">printf</span><span style="color: #0000FF;">(</span><span style="color: #000000;">1</span><span style="color: #0000FF;">,</span><span style="color: #008000;">"Root found %s %.10g\n"</span><span style="color: #0000FF;">,</span> <span style="color: #0000FF;">{</span><span style="color: #008080;">iff</span><span style="color: #0000FF;">(</span><span style="color: #000000;">y</span><span style="color: #0000FF;">=</span><span style="color: #000000;">0</span><span style="color: #0000FF;">?</span><span style="color: #008000;">"at"</span><span style="color: #0000FF;">:</span><span style="color: #008000;">"near"</span><span style="color: #0000FF;">),</span><span style="color: #000000;">x</span><span style="color: #0000FF;">})</span>
<span style="color: #008080;">end</span> <span style="color: #008080;">if</span>
<span style="color: #000000;">x</span> <span style="color: #0000FF;">+=</span> <span style="color: #000000;">step</span>
<span style="color: #008080;">end</span> <span style="color: #008080;">while</span>
<span style="color: #008080;">end</span> <span style="color: #008080;">procedure</span>
<span style="color: #000080;font-style:italic;">-- 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.</span>
<span style="color: #008080;">constant</span> <span style="color: #000000;">step</span> <span style="color: #0000FF;">=</span> <span style="color: #000000;">1</span><span style="color: #0000FF;">/</span><span style="color: #000000;">256</span>
<span style="color: #008080;">function</span> <span style="color: #000000;">f1</span><span style="color: #0000FF;">(</span><span style="color: #004080;">atom</span> <span style="color: #000000;">x</span><span style="color: #0000FF;">)</span> <span style="color: #008080;">return</span> <span style="color: #000000;">x</span><span style="color: #0000FF;">*</span><span style="color: #000000;">x</span><span style="color: #0000FF;">*</span><span style="color: #000000;">x</span><span style="color: #0000FF;">-</span><span style="color: #000000;">3</span><span style="color: #0000FF;">*</span><span style="color: #000000;">x</span><span style="color: #0000FF;">*</span><span style="color: #000000;">x</span><span style="color: #0000FF;">+</span><span style="color: #000000;">2</span><span style="color: #0000FF;">*</span><span style="color: #000000;">x</span> <span style="color: #008080;">end</span> <span style="color: #008080;">function</span>
<span style="color: #008080;">function</span> <span style="color: #000000;">f2</span><span style="color: #0000FF;">(</span><span style="color: #004080;">atom</span> <span style="color: #000000;">x</span><span style="color: #0000FF;">)</span> <span style="color: #008080;">return</span> <span style="color: #000000;">x</span><span style="color: #0000FF;">*</span><span style="color: #000000;">x</span><span style="color: #0000FF;">-</span><span style="color: #000000;">4</span><span style="color: #0000FF;">*</span><span style="color: #000000;">x</span><span style="color: #0000FF;">+</span><span style="color: #000000;">3</span> <span style="color: #008080;">end</span> <span style="color: #008080;">function</span>
<span style="color: #008080;">function</span> <span style="color: #000000;">f3</span><span style="color: #0000FF;">(</span><span style="color: #004080;">atom</span> <span style="color: #000000;">x</span><span style="color: #0000FF;">)</span> <span style="color: #008080;">return</span> <span style="color: #000000;">x</span><span style="color: #0000FF;">-</span><span style="color: #000000;">1.5</span> <span style="color: #008080;">end</span> <span style="color: #008080;">function</span>
<span style="color: #008080;">function</span> <span style="color: #000000;">f4</span><span style="color: #0000FF;">(</span><span style="color: #004080;">atom</span> <span style="color: #000000;">x</span><span style="color: #0000FF;">)</span> <span style="color: #008080;">return</span> <span style="color: #000000;">x</span><span style="color: #0000FF;">*</span><span style="color: #000000;">x</span><span style="color: #0000FF;">-</span><span style="color: #000000;">2</span> <span style="color: #008080;">end</span> <span style="color: #008080;">function</span>
<span style="color: #000000;">print_roots</span><span style="color: #0000FF;">(</span><span style="color: #000000;">f1</span><span style="color: #0000FF;">,</span> <span style="color: #0000FF;">-</span><span style="color: #000000;">1</span><span style="color: #0000FF;">,</span> <span style="color: #000000;">5</span><span style="color: #0000FF;">,</span> <span style="color: #000000;">step</span><span style="color: #0000FF;">)</span>
<span style="color: #000000;">print_roots</span><span style="color: #0000FF;">(</span><span style="color: #000000;">f2</span><span style="color: #0000FF;">,</span> <span style="color: #0000FF;">-</span><span style="color: #000000;">1</span><span style="color: #0000FF;">,</span> <span style="color: #000000;">5</span><span style="color: #0000FF;">,</span> <span style="color: #000000;">step</span><span style="color: #0000FF;">)</span>
<span style="color: #000000;">print_roots</span><span style="color: #0000FF;">(</span><span style="color: #000000;">f3</span><span style="color: #0000FF;">,</span> <span style="color: #000000;">0</span><span style="color: #0000FF;">,</span> <span style="color: #000000;">4</span><span style="color: #0000FF;">,</span> <span style="color: #000000;">step</span><span style="color: #0000FF;">)</span>
<span style="color: #000000;">print_roots</span><span style="color: #0000FF;">(</span><span style="color: #000000;">f4</span><span style="color: #0000FF;">,</span> <span style="color: #0000FF;">-</span><span style="color: #000000;">2</span><span style="color: #0000FF;">,</span> <span style="color: #000000;">2</span><span style="color: #0000FF;">,</span> <span style="color: #000000;">step</span><span style="color: #0000FF;">)</span>
<!--</syntaxhighlight>-->
{{out}}
<pre>
-----
Line 2,393 ⟶ 2,521:
=={{header|PicoLisp}}==
{{trans|Clojure}}
<
(filter
'((N) (> Eps (abs (F N))))
Line 2,403 ⟶ 2,531:
(findRoots
'((X) (+ (*/ X X X `(* 1.0 1.0)) (*/ -3 X X 1.0) (* 2 X)))
-1.0 3.0 0.0001 0.00000001 ) )</
Output:
<pre>-> ("0.000" "1.000" "2.000")</pre>
=={{header|PL/I}}==
<syntaxhighlight lang="pl/i">
f: procedure (x) returns (float (18));
declare x float (18);
Line 2,450 ⟶ 2,578:
end;
end locate_root;
</syntaxhighlight>
=={{header|PureBasic}}==
{{trans|C++}}
<
ProcedureReturn x*x*x-3*x*x+2*x
EndProcedure
Line 2,481 ⟶ 2,609:
EndProcedure
main()</
=={{header|Python}}==
{{trans|Perl}}
<
step = 0.001 # Smaller step values produce more accurate and precise results
Line 2,507 ⟶ 2,635:
sign = value > 0
x += step</
=={{header|R}}==
{{trans|Octave}}
<
findroots <- function(f, begin, end, tol = 1e-20, step = 0.001) {
Line 2,528 ⟶ 2,656:
}
findroots(f, -1, 3)</
=={{header|Racket}}==
<
#lang racket
Line 2,572 ⟶ 2,700:
(define tolerance (make-parameter 5e-16))
</syntaxhighlight>
Different root-finding methods
<
(define (secant f a b)
(let next ([x1 a] [y1 (f a)] [x2 b] [y2 (f b)] [n 50])
Line 2,594 ⟶ 2,722:
c
(or (divide a c) (divide c b)))))))
</syntaxhighlight>
Examples:
<
-> (find-root (λ (x) (- 2. (* x x))) 1 2)
1.414213562373095
Line 2,606 ⟶ 2,734:
-> (find-roots f -3 4 #:divisions 50)
'(2.4932181969624796e-33 1.0 2.0)
</syntaxhighlight>
In order to provide a comprehensive code the given solution does not optimize the number of function calls.
Line 2,612 ⟶ 2,740:
Simple memoization operator
<
(define (memoized f)
(define tbl (make-hash))
Line 2,620 ⟶ 2,748:
(hash-set! tbl x res)
res])))
</
To use memoization just call
<
-> (find-roots (memoized f) -3 4 #:divisions 50)
'(2.4932181969624796e-33 1.0 2.0)
</syntaxhighlight>
The profiling shows that memoization reduces the number of function calls
Line 2,634 ⟶ 2,762:
(formerly Perl 6)
Uses exact arithmetic.
<syntaxhighlight lang="raku"
my $start = -1;
Line 2,652 ⟶ 2,780:
NEXT $sign = $next;
}
}</
{{out}}
<pre>Root found at 0
Line 2,661 ⟶ 2,789:
Both of these REXX versions use the '''bisection method'''.
===function coded as a REXX function===
<
parse arg bot top inc . /*obtain optional arguments from the CL*/
if bot=='' | bot=="," then bot= -5 /*Not specified? Then use the default.*/
Line 2,678 ⟶ 2,806:
f: parse arg x; return x * (x * (x-3) +2) /*formula used ──► x^3 - 3x^2 + 2x */
/*with factoring ──► x{ x^2 -3x + 2 } */
/*more " ──► x{ x( x-3 ) + 2 } */</
{{out|output|text= when using the defaults for input:}}
<pre>
Line 2,688 ⟶ 2,816:
===function coded in-line===
This version is about '''40%''' faster than the 1<sup>st</sup> REXX version.
<
parse arg bot top inc . /*obtain optional arguments from the CL*/
if bot=='' | bot=="," then bot= -5 /*Not specified? Then use the default.*/
Line 2,701 ⟶ 2,829:
else if !\==$ then if !\==0 then say 'passed a root at' x/1
!= $ /*use the new sign for the next compare*/
end /*x*/ /*dividing by unity normalizes X [↑] */</
{{out|output|text= is the same as the 1<sup>st</sup> REXX version.}} <br><br>
=={{header|Ring}}==
<
load "stdlib.ring"
function = "return pow(x,3)-3*pow(x,2)+2*x"
Line 2,727 ⟶ 2,855:
oldsign = num
next
</syntaxhighlight>
Output:
<pre>
Line 2,741 ⟶ 2,869:
The script that finds a root of a scalar function <math>f(x) = x^3-3\,x^2 + 2\,x</math> of a scalar variable ''x''
using the bisection method on the interval -5 to 5 is,
<syntaxhighlight lang="rlab">
f = function(x)
{
Line 2,750 ⟶ 2,878:
>> findroot(f, , [-5,5])
0
</syntaxhighlight>
For a detailed description of the solver and its parameters interested reader is directed to the ''rlabplus'' manual.
Line 2,757 ⟶ 2,885:
{{trans|Python}}
<
x <=> 0
end
Line 2,775 ⟶ 2,903:
f = lambda { |x| x**3 - 3*x**2 + 2*x }
find_roots(f, -1..3)</
{{out}}
Line 2,786 ⟶ 2,914:
Or we could use Enumerable#inject, monkey patching and block:
<
def sign
self <=> 0
Line 2,804 ⟶ 2,932:
end
find_roots(-1..3) { |x| x**3 - 3*x**2 + 2*x }</
=={{header|Rust}}==
<syntaxhighlight lang="rust">// 202100315 Rust programming solution
use roots::find_roots_cubic;
fn main() {
let roots = find_roots_cubic(1f32, -3f32, 2f32, 0f32);
println!("Result : {:?}", roots);
}</syntaxhighlight>
{{out}}
<pre>
Result : Three([0.000000059604645, 0.99999994, 2.0])
</pre>
Another without external crates:
<syntaxhighlight lang="rust">
use num::Float;
/// Note: We cannot use `range_step` here because Floats don't implement
/// the `CheckedAdd` trait.
fn find_roots<T, F>(f: F, start: T, stop: T, step: T, epsilon: T) -> Vec<T>
where
T: Copy + PartialOrd + Float,
F: Fn(T) -> T,
{
let mut ret = vec![];
let mut current = start;
while current < stop {
if f(current).abs() < epsilon {
ret.push(current);
}
current = current + step;
}
ret
}
fn main() {
let roots = find_roots(
|x: f64| x * x * x - 3.0 * x * x + 2.0 * x,
-1.0,
3.0,
0.0001,
0.00000001,
);
println!("roots of f(x) = x^3 - 3x^2 + 2x are: {:?}", roots);
}
</syntaxhighlight>
{{out}}
<pre>
roots of f(x) = x^3 - 3x^2 + 2x are: [-0.00000000000009381755897326649, 0.9999999999998124, 1.9999999999997022]
</pre>
=={{header|Scala}}==
Line 2,810 ⟶ 2,994:
{{trans|Java}}
{{Out}}Best seen running in your browser either by [https://scalafiddle.io/sf/T63KUsH/0 (ES aka JavaScript, non JVM)] or [https://scastie.scala-lang.org/bh8von94Q1y0tInvEZ3cBQ Scastie (remote JVM)].
<
val poly = (x: Double) => x * x * x - 3 * x * x + 2 * x
Line 2,834 ⟶ 3,018:
printRoots(poly, -1.0, 4, 0.002)
}</
===Functional version (Recommended)===
<
def findRoots(fn: Double => Double, start: Double, stop: Double, step: Double, epsilon: Double) = {
for {
Line 2,847 ⟶ 3,031:
println(findRoots(fn, -1.0, 3.0, 0.0001, 0.000000001))
}</
{{out}}
Vector(-9.381755897326649E-14, 0.9999999999998124, 1.9999999999997022)
=={{header|Scheme}}==
For R7RS Scheme.
<syntaxhighlight lang="scheme">
;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;
;;;
;;; Finding (real) roots of a function.
;;;
;;; I follow the model that breaks the task into two distinct ones:
;;; isolating real roots, and then finding the isolated roots. The
;;; former task I will call "isolating roots", and the latter I will
;;; call "rootfinding".
;;;
;;; Isolating real roots of a polynomial can be done exactly, and the
;;; methods can handle infinite domains. Scheme (because it has exact
;;; rationals) is a relatively easy language in which to write such
;;; code.
;;;
;;; I have also isolated the roots of low-degree polynomials on finite
;;; intervals by the following method, in floating-point arithmetic:
;;; rewrite the polynomial in the Bernstein polynomials basis, then
;;; take derivatives to get critical points, working your way back up
;;; in degree from a straight line. This method goes back and forth
;;; between the two subtasks. (You could use the quadratic formula
;;; once you got down to degree two, but I wouldn’t bother. The cubic
;;; and quartic formulas are numerically very poor and should be
;;; avoided.)
;;;
;;; However, these methods require that the function be a
;;; polynomial. Here I will simply use a step size and the
;;; intermediate value theorem.
;;;
;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;
(cond-expand
(r7rs)
(chicken (import r7rs)))
;;;
;;; Step 1. Isolation of roots.
;;;
;;; I will simply step over the domain and look for intervals that
;;; contain at least one root. There is a risk of getting the error
;;; message "the root is not bracketed" when you run the
;;; rootfinder. (One could have the rootfinder raise a recoverable
;;; exception or return a special value, instead.)
;;;
(define-library (isolate-roots)
(export isolate-roots)
(import (scheme base))
(begin
(define (isolate-roots f x-min x-max x-step)
(define (ith-x i) (+ x-min (* i x-step)))
(let ((x0 x-min)
(y0 (f x-min)))
(let loop ((i1 1)
(x0 x0)
(y0 y0)
(accum (if (zero? y0)
`((,x0 . ,x0))
'())))
(let ((x1 (ith-x i1)))
(if (< x-max x1)
(reverse accum)
(let ((y1 (f x1)))
(cond
((zero? y1)
(loop (+ i1 1) x1 y1 `((,x1 . ,x1) . ,accum)))
((negative? (* y0 y1))
(loop (+ i1 1) x1 y1 `((,x0 . ,x1) . ,accum)))
(else
(loop (+ i1 1) x1 y1 accum)))))))))
)) ;; end library roots-isolator
;;;
;;; Step 2. Rootfinding.
;;;
;;; I will use the ITP method. I wrote this implementation shortly
;;; after the algorithm was published. See
;;; https://sourceforge.net/p/chemoelectric/itp-root-finder
;;;
;;; Reference:
;;;
;;; I.F.D. Oliveira and R.H.C. Takahashi. 2020. An Enhancement of
;;; the Bisection Method Average Performance Preserving Minmax
;;; Optimality. ACM Trans. Math. Softw. 47, 1, Article 5 (December
;;; 2020), 24 pages. https://doi.org/10.1145/3423597
;;;
(define-library (itp-root-finder)
(export itp-root-finder-epsilon
itp-root-finder-extra-steps
itp-root-finder-kappa1
itp-root-finder-kappa2
;; itp-root-bracket-finder returns two values that form a
;; bracket no wider than 2ϵ.
itp-root-bracket-finder
;; itp-root-finder returns the point midway between the ends
;; of the final bracket.
itp-root-finder)
(import (scheme base))
(import (scheme inexact))
(import (scheme case-lambda))
(import (srfi 143)) ; Fixnums.
(import (srfi 144)) ; Flonums.
(begin
(define ϕ
;; The Golden Ratio, (1 + √5)/2, rounded down by about
;; 0.00003398875.
1.618)
(define 1+ϕ (+ 1.0 ϕ))
(define itp-root-finder-epsilon
(make-parameter
(* 1000.0 fl-epsilon)
(lambda (ϵ)
(if (positive? ϵ)
ϵ
(error 'itp-root-finder-epsilon
"a positive value was expected"
ϵ)))))
(define itp-root-finder-extra-steps
;; Increase extra-steps above zero, if you wish to try to speed
;; up convergence, at the expense of that many more steps in the
;; worst case.
(make-parameter
0
(lambda (n₀)
(if (or (negative? n₀) (not (integer? n₀)))
(error 'itp-root-finder-extra-steps
"a non-negative integer was expected"
n₀)
n₀))))
(define itp-root-finder-kappa1
(make-parameter
0.1
(lambda (κ₁)
(if (positive? κ₁)
κ₁
(error 'itp-root-finder-kappa1
"a positive value was expected"
κ₁)))))
(define itp-root-finder-kappa2
(make-parameter
2.0
(lambda (κ₂)
(if (or (< κ₂ 1) (< 1+ϕ κ₂))
;; We allow <= 1+ϕ (instead of ‘< 1+ϕ’) because we
;; already rounded ϕ down.
(error 'itp-root-finder-kappa2
(string-append "a value 1 <= kappa2 <= "
(number->string 1+ϕ)
" was expected")
κ₂)
κ₂))))
(define (sign x)
(cond ((negative? x) -1)
((positive? x) 1)
(else 0)))
(define (apply-sign σ x)
(cond ((fxnegative? σ) (- x))
((fxpositive? σ) x)
(else 0)))
(define (itp-root-bracket-finder%% f a b ϵ n₀ κ₁ κ₂)
(let* ((2ϵ (inexact (* 2 ϵ)))
(n½ (exact (ceiling (log (/ (inexact (- b a)) 2ϵ) 2))))
(n_max (+ n½ n₀))
(ya (f a))
(yb (f b))
(σ_ya (sign ya))
(σ_yb (sign yb)))
(cond
((fxzero? σ_ya) (values a a))
((fxzero? σ_yb) (values b b))
(else
(when (fxpositive? (* σ_ya σ_yb))
(error 'itp-root-bracket-finder
"the root is not bracketed"
a b))
(let loop ((pow2 (expt 2 n_max))
(a (inexact a))
(b (inexact b))
(ya (inexact ya))
(yb (inexact yb)))
(if (or (= pow2 1) (fl<=? (fl- b a) 2ϵ))
(values a b)
(let* ( ;; x½ – the bisection.
(x½ (fl* 0.5 (fl+ a b)))
;; xf – interpolation by regula falsi.
(xf (fl/ (fl- (fl* yb a) (fl* ya b))
(fl- yb ya)))
(b-a (fl- b a))
(δ (fl* κ₁ (flabs (expt b-a κ₂))))
(x½-xf (fl- x½ xf))
(σ (sign x½-xf))
;; xt – the ‘truncation’ of xf.
(xt (if (fl<=? δ (flabs x½-xf))
(fl+ xf (apply-sign σ δ))
x½))
(r (- (* pow2 ϵ) (fl* 0.5 b-a)))
;; xp – the projection of xt onto [x½-r,x½+r].
(xp (if (fl<=? (flabs (fl- xt x½)) r)
xt
(fl- x½ (apply-sign σ r))))
(yp (inexact (f xp))))
(let ((pow2/2 (truncate-quotient pow2 2))
(σ_yp (sign yp)))
(cond ((fx=? σ_ya σ_yp)
;; yp has the same sign as ya. Make it the
;; new ya.
(loop pow2/2 xp b yp yb))
((fx=? σ_yb σ_yp)
;; yp has the same sign as yb. Make it the
;; new yb.
(loop pow2/2 a xp ya yp))
(else
;; yp is zero.
(values xp xp)))))))))))
(define (itp-root-bracket-finder% f a b ϵ n₀ κ₁ κ₂)
(cond
((< b a) (itp-root-bracket-finder% b a f ϵ n₀ κ₁ κ₂))
(else
(let* ((ϵ (or ϵ (itp-root-finder-epsilon)))
(n₀ (or n₀ (itp-root-finder-extra-steps)))
(κ₁ (or κ₁ (itp-root-finder-kappa1)))
(κ₂ (or κ₂ (itp-root-finder-kappa2))))
(when (negative? ϵ)
(error 'itp-root-bracket-finder
"a positive value was expected" ϵ))
(when (negative? κ₁)
(error 'itp-root-bracket-finder
"a positive value was expected" κ₁))
(when (or (< κ₂ 1) (< 1+ϕ κ₂))
;; We allow <= 1+ϕ (instead of ‘< 1+ϕ’) because we already
;; rounded ϕ down.
(error 'itp-root-bracket-finder
(string-append "a value 1 <= kappa2 <= "
(number->string 1+ϕ)
" was expected")
κ₂))
(when (or (negative? n₀) (not (integer? n₀)))
(error 'itp-root-bracket-finder
"a non-negative integer was expected" n₀))
(itp-root-bracket-finder%% f a b ϵ n₀ κ₁ κ₂)))))
(define itp-root-bracket-finder
(case-lambda
((f a b)
(itp-root-bracket-finder% f a b #f #f #f #f))
((f a b ϵ)
(itp-root-bracket-finder% f a b ϵ #f #f #f))
((f a b ϵ n₀)
(itp-root-bracket-finder% f a b ϵ n₀ #f #f))
((f a b ϵ n₀ κ₁)
(itp-root-bracket-finder% f a b ϵ n₀ κ₁ #f))
((f a b ϵ n₀ κ₁ κ₂)
(itp-root-bracket-finder% f a b ϵ n₀ κ₁ κ₂))))
(define (itp-root-finder% f a b ϵ n₀ κ₁ κ₂)
(call-with-values
(lambda ()
(itp-root-bracket-finder f a b ϵ n₀ κ₁ κ₂))
(lambda (a b)
(if (= a b)
a
(* 0.5 (+ a b))))))
(define itp-root-finder
(case-lambda
((f a b)
(itp-root-finder% f a b #f #f #f #f))
((f a b ϵ)
(itp-root-finder% f a b ϵ #f #f #f))
((f a b ϵ n₀)
(itp-root-finder% f a b ϵ n₀ #f #f))
((f a b ϵ n₀ κ₁)
(itp-root-finder% f a b ϵ n₀ κ₁ #f))
((f a b ϵ n₀ κ₁ κ₂)
(itp-root-finder% f a b ϵ n₀ κ₁ κ₂))))
)) ;; end library itp-root-finder
(import (scheme base))
(import (scheme write))
(import (isolate-roots))
(import (itp-root-finder))
(define (f x)
;; x³ - 3x² + 2x, written in Horner form.
(* x (+ 2 (* x (+ -3 x)))))
(define (find-root f interval)
(define (display-exactness root)
(display (if (and (exact? root)
(exact? (f root))
(zero? (f root)))
" (exact) "
" (inexact) ")))
(let ((x0 (car interval))
(x1 (cdr interval)))
(if (= x0 x1)
(begin
(let ((root (if (exact? x0) x0 x1)))
(display-exactness root)
(display "(rootfinder not used) ")
(display root)
(newline)))
(begin
;;
;; I am not careful here to avoid accidentally excluding the
;; root from the bracketing interval [x0,x1]. Floating point
;; is very tricky to work with.
;;
(let ((root (itp-root-finder f x0 x1)))
(display-exactness root)
(display "(rootfinder used) ")
(display root)
(newline))))))
;;; The following two demonstrations find all three roots exactly, as
;;; exact rationals, without the need for a rootfinding step.
(newline)
(display "Stepping by 1/1000 from 0 to 2:")
(newline)
(do ((p (isolate-roots f 0 2 1/1000) (cdr p)))
((not (pair? p)))
(find-root f (car p)))
(newline)
(display "Stepping by 1/1000 from -10 to 10:")
(newline)
(do ((p (isolate-roots f -10 10 1/1000) (cdr p)))
((not (pair? p)))
(find-root f (car p)))
;;; The following demonstration gives inexact results, because the
;;; step size is an inexact number.
(newline)
(display "Stepping by 0.001 from -10.0 to 10.0:")
(newline)
(do ((p (isolate-roots f -10.0 10.0 0.001) (cdr p)))
((not (pair? p)))
(find-root f (car p)))
;;; The following demonstration gives inexact results, because the
;;; rootfinder is needed.
(newline)
(display "Stepping by 13/3333 from -2111/1011 to 33/13:")
(newline)
(do ((p (isolate-roots f -2111/1011 33/13 13/3333) (cdr p)))
((not (pair? p)))
(find-root f (car p)))
(newline)
</syntaxhighlight>
{{out}}
<pre>$ gosh roots-of-a-function.scm
Stepping by 1/1000 from 0 to 2:
(exact) (rootfinder not used) 0
(exact) (rootfinder not used) 1
(exact) (rootfinder not used) 2
Stepping by 1/1000 from -10 to 10:
(exact) (rootfinder not used) 0
(exact) (rootfinder not used) 1
(exact) (rootfinder not used) 2
Stepping by 0.001 from -10.0 to 10.0:
(inexact) (rootfinder not used) 0.0
(inexact) (rootfinder not used) 1.0
(inexact) (rootfinder not used) 2.0
Stepping by 13/3333 from -2111/1011 to 33/13:
(inexact) (rootfinder used) -1.3380129580295458e-15
(inexact) (rootfinder used) 0.9999999999998657
(inexact) (rootfinder used) 1.9999999999999998
</pre>
=={{header|Sidef}}==
<
x*x*x - 3*x*x + 2*x
}
Line 2,871 ⟶ 3,463:
}
sign = value>0
}</
{{out}}
<pre>Root found at 0
Line 2,879 ⟶ 3,471:
=={{header|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).
<
set res {}
set lastsign [sgn [apply $lambda $start]]
Line 2,893 ⟶ 3,485:
proc sgn x {expr {($x>0) - ($x<0)}}
puts [froots {x {expr {$x**3 - 3*$x**2 + 2*$x}}}]</
Result and timing:
<pre>/Tcl $ time ./froots.tcl
Line 2,902 ⟶ 3,494:
sys 0m0.030s</pre>
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 [[wp:Newton's method|Newton-Raphson]], but that requires the differential of the function with respect to the search variable.
<
set res {}
set lastsign [sgn [apply $f $start]]
Line 2,930 ⟶ 3,522:
puts [frootsNR \
{x {expr {$x**3 - 3*$x**2 + 2*$x}}} \
{x {expr {3*$x**2 - 6*$x + 2}}}]</
=={{header|TI-89 BASIC}}==
Line 2,941 ⟶ 3,533:
{{trans|Go}}
{{libheader|Wren-fmt}}
<
var secant = Fn.new { |f, x0, x1|
Line 2,975 ⟶ 3,567:
var example = Fn.new { |x| x*x*x - 3*x*x + 2*x }
findRoots.call(example, -0.5, 2.6, 1)</
{{out}}
<pre>
0.000 approximate
1.000 exact
2.000 approximate
</pre>
=={{header|XPL0}}==
{{trans|Wren}}
<syntaxhighlight lang "XPL0">include xpllib; \for Print
func real F(X);
real X;
return X*X*X - 3.*X*X + 2.*X;
char Status;
func real Secant(X0, X1);
real X0, X1, F0, F1, T;
int I;
[F1:= F(X0);
for I:= 0 to 100-1 do
[F0:= F1;
F1:= F(X1);
if F1 = 0. then [Status:= "exact"; return X1];
if abs(X1-X0) < 1e-6 then [Status:= "approximate"; return X1];
T:= X0;
X0:= X1;
X1:= X1 - F1*(X1-T)/(F1-F0);
];
Status:= 0; return 0.;
];
func FindRoots(Lower, Upper, Step);
real Lower, Upper, Step;
real X0, X1, R;
[X0:= Lower;
X1:= Lower + Step;
while X0 < Upper do
[X1:= if X1 < Upper then X1 else Upper;
R:= Secant(X0, X1);
if Status # 0 and R >= X0 and R < X1 then
Print(" %2.3f %s\n", R, Status);
X0:= X1;
X1:= X1 + Step;
];
];
FindRoots(-0.5, 2.6, 1.)</syntaxhighlight>
{{out}}
<pre>
Line 2,986 ⟶ 3,627:
=={{header|zkl}}==
{{trans|Haskell}}
<
[start..stop,step].filter('wrap(x){ f(x).closeTo(0.0,eps) })
}</
<
findRoots(f, -1.0, 3.0, 0.0001, 0.00000001).println();</
{{out}}
<pre>L(-9.38176e-14,1,2)</pre>
{{trans|C}}
<
reg e=1.0e-12;
Line 3,007 ⟶ 3,648:
if(f(xB).closeTo(0.0,e)) xB
else "Function is not converging near (%7.4f,%7.4f).".fmt(xA,xB);
}</
<
xs:=findRoots(f, -1.032, 3.0, step, 0.1);
xs.println(" --> ",xs.apply('wrap(x){ secant(f,x-step,x+step) }));</
{{out}}
<pre>L(-0.032,0.968,1.068,1.968) --> L(1.87115e-19,1,1,2)</pre>
|