Roots of a function: Difference between revisions

{{trans|Python}}

 

R x^3 - 3 * x^2 + 2 * x

 

 

sgn = value > 0

 

{{out}}

 

=={{header|Ada}}==

procedure Roots_Of_Function is

X := X + Step;

end loop;

 

=={{header|ALGOL 68}}==

{{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)$;

 

)

OUT printf($"No third root found"l$); stop

Output:

<pre>1st root found at x = 9.1557112297752398099031e-1 (Approximately)

3rd root found at x = 0.0000000000000000000000e 0 (Exactly)

</pre>

 

=={{header|ATS}}==

#include

"share/atspre_staload.hats"

main0 () =

findRoots (~1.0, 3.0, 0.001, lam (x) => x*x*x - 3.0*x*x + 2.0*x, 0, 0.0)

 

=={{header|AutoHotkey}}==

 

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

 

poly(x) {

Return ((x-3)*x+2)*x

 

=={{header|Axiom}}==

Using a polynomial solver:

Output:

Using the secant method in the interpreter:

secant(eq: Equation Expression Float, binding: SegmentBinding(Float)):Float ==

eps := 1.0e-30

abs(fx2)<eps => return x2

(x1, fx1, x2) := (x2, fx2, x2 - fx2 * (x2 - x1) / (fx2 - fx1))

The example can now be called using:

 

=={{header|BBC BASIC}}==

rangemin = -1

rangemax = 3

oldsign% = sign%

NEXT x

Output:

<pre>Root found near x = 2.29204307E-9

=== Secant Method ===

 

#include <stdio.h>

 

}

return 0;

 

=== GNU Scientific Library ===

 

#include <stdio.h>

 

 

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>

{{trans|C++}}

 

 

class Program

}

}

 

{{trans|Java}}

 

class Program

PrintRoots(f, -1.0, 4, 0.002);

}

 

===Brent's Method===

 

{{trans|C++}}

 

class Program

// end brents_fun

}

 

=={{header|C++}}==

 

double f(double x)

sign = ( value > 0 );

}

 

===Brent's Method===

 

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>

} // end brents_fun

 

 

=={{header|Clojure}}==

 

{{trans|Haskell}}

 

(defn findRoots [f start stop step eps]

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

 

<pre>

=={{header|CoffeeScript}}==

{{trans|Python}}

print_roots = (f, begin, end, step) ->

# Print approximate roots of f between x=begin and x=end,

f4 = (x) -> x*x - 2

print_roots f4, -2, 2, step

 

output

 

> coffee roots.coffee

-----

Root found near -1.4140625

Root found near 1.41796875

 

=={{header|Common Lisp}}==

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

(format t "~&Root found near ~w." x)

(push (cons (- x step) x) roots)))

 

<pre>> (find-roots #'(lambda (x) (+ (* x x x) (* -3 x x) (* 2 x))) -1 3)

 

=={{header|D}}==

 

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

 

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

{{out}}

<pre>Root found (tolerance = 0.0001):

=={{header|Dart}}==

{{trans|Scala}}

 

findRoots(Function(double) f, double start, double stop, double step, double epsilon) sync* {

// Vector(-9.381755897326649E-14, 0.9999999999998124, 1.9999999999997022)

print(findRoots(fn, -1.0, 3.0, 0.0001, 0.000000001));

 

=={{header|Delphi}}==

=={{header|DWScript}}==

{{trans|C}}

 

function f(x : Float) : Float;

end;

x += fstep;

 

=={{header|EchoLisp}}==

We use the 'math' library, and define f(x) as the polynomial : x³ -3x² +2x

 

(lib 'math.lib)

Lib: math.lib loaded.

(root f -1000 (- 2 epsilon)) → 1.385559938161431e-7 ;; 0

(root f epsilon (- 2 epsilon)) → 1.0000000002190812 ;; 1

 

=={{header|Elixir}}==

{{trans|Ruby}}

def find_roots(f, range, step \\ 0.001) do

first .. last = range

 

f = fn x -> x*x*x - 3*x*x + 2*x end

 

{{out}}

 

=={{header|Erlang}}==

-module(roots).

-export([main/0]).

while(X+Step, Step, Start, Stop, Value > 0,F)

end.

{{out}}

<pre>Root found near 8.81239525796218e-16

 

=={{header|ERRE}}==

PROGRAM ROOTS_FUNCTION

 

END LOOP

END PROGRAM

Note: Outputs are calculated in single precision.

{{out}}

=={{header|Fortran}}==

{{works with|Fortran|90 and later}}

 

IMPLICIT NONE

END DO

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

x2 = x2 - d

i = i + 1

 

=={{header|Go}}==

Secant method. No error checking.

 

import (

}

return 0, ""

Output:

<pre>

 

=={{header|Haskell}}==

 

findRoots start stop step eps =

Executed in GHCi:

 

Or using package [http://hackage.haskell.org/package/hmatrix hmatrix] from HackageDB.

import Data.Complex

 

(-5.421010862427522e-20) :+ 0.0

2.000000000000001 :+ 0.0

No complex roots, so:

-5.421010862427522e-20

2.000000000000001

 

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)

findRoots f [] = empty

findRoots f [x] = if f x == 0 then pure (Exact x) else empty

 

It is possible to use these functions with different Alternative functors: IO, Maybe or List:

=={{header|HicEst}}==

HicEst's [http://www.HicEst.com/SOLVE.htm SOLVE] function employs the Levenberg-Marquardt method:

 

1 DLG(NameEdit=x0, DNum=3)

 

WRITE(FIle='test.txt', LENgth=6, Name) x0, x, solution, chi2, iterations

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=1.55; x=2; solution=exact; chi2=79E-32 iterations=55;

x0=1E10; x=2; solution=exact; chi2=18E-31 iterations=511;

 

=={{header|Icon}} and {{header|Unicon}}==

 

Works in both languages:

showRoots(f, -1.0, 4, 0.002)

end

procedure sign(x)

return (x<0, -1) | (x>0, 1) | 0

 

Output:

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:

 

 

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, <tt>p.</tt> 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.

 

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

0 1 2

0=blackbox 0 1 2

 

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

printRoots(poly, -1.0, 4, 0.002);

}

Produces this output:

<pre>~2.616794878713638E-18

{{works with|SpiderMonkey|22}}

{{works with|Firefox|22}}

// This function notation is sorta new, but useful here

// Part of the EcmaScript 6 Draft

 

printRoots(poly, -1.0, 4, 0.002);

 

=={{header|jq}}==

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;

 

end )

| .[3] ;

We present two examples, one where step is a power of 1/2, and one where it is not:

{{Out}}

 

[

"~1.0000000000000002",

"~1.9999999999999993"

 

=={{header|Julia}}==

Assuming that one has the Roots package installed:

 

 

 

{{out}}

 

Without the Roots package, Newton's method may be defined in this manner:

##f: the function of x

##fp: the derivative of f

end

end

 

Finding the roots of f(x) = x3 - 3x2 + 2x:

 

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)

{{out}}

 

=={{header|Kotlin}}==

{{trans|C}}

 

typealias DoubleToDouble = (Double) -> Double

x += step

}

 

{{out}}

 

=={{header|Lambdatalk}}==

1) defining the function:

{def func {lambda {:x} {+ {* 1 :x :x :x} {* -3 :x :x} {* 2 :x}}}}

- a root found at 540°

- a root found at 720°

 

=={{header|Liberty BASIC}}==

' within a range of x values.

 

function f( x)

f =x^3 -3 *x^2 +2 *x

 

=={{header|Lua}}==

function f (x) return x^3 - 3*x^2 + 2*x end

 

for _, r in pairs(root(f, -1, 3, 10^-6)) do

print(string.format("%0.12f", r.val), r.err)

{{out}}

<pre>Root (to 12DP) Max. Error

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)

{{out}}

<pre>Root (to 12DP) Max. Error

=={{header|Maple}}==

 

 

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.

 

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>

===NSolve===

This requires merely the polynomial and will perform numerical operations if needed:

Output

<pre> {{x->0.},{x->1.},{x->2.}}</pre>

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

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

===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 */

[x=1.16154139999725193608791768724717407484314725802151429063617b0*%i + 3.41163901914009663684741869855524128445594290948999288901864b-1,

x=3.41163901914009663684741869855524128445594290948999288901864b-1 - 1.16154139999725193608791768724717407484314725802151429063617b0*%i,

 

=={{header|Nim}}==

import strformat

 

sign = value > 0

 

{{out}}

=={{header|Objeck}}==

{{trans|C++}}

bundle Default {

class Roots {

}

}

 

=={{header|OCaml}}==

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

 

x fx (if fx = 0.0 then "exact" else "approx") in

let f x = ((x -. 3.0)*.x +. 2.0)*.x in

 

Output:

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

 

r = roots(a);

% let's print it

endif

printf(" exact)\n");

 

Otherwise we can program our (simple) method:

 

{{trans|Python}}

y = x.^3 -3.*x.^2 + 2.*x;

endfunction

se = sign(v);

x = x + step;

 

=={{header|Oforth}}==

 

| x y lasty |

a f perform dup ->y ->lasty

] ;

 

 

{{out}}

 

=={{header|ooRexx}}==

a*x**3+b*x**2+c*x+d =(x-x1)*(x-x2)*(x-x3)

*/

Else

Return xr

{{out}}

<pre>p=-3

=={{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)

 

===Factorization to linear factors===

my(f=factor(P),t);

for(i=1,#f[,1],

)

};

 

===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],

)

};

 

=={{header|Pascal}}==

{{trans|Fortran}}

 

var

end;

end.

Output:

<pre>

 

=={{header|Perl}}==

{

my $x = shift;

# Update our sign

$sign = ( $value > 0 );

 

=={{header|Phix}}==

{{trans|CoffeeScript}}

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

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

{{out}}

<pre>

=={{header|PicoLisp}}==

{{trans|Clojure}}

(filter

'((N) (> Eps (abs (F N))))

(findRoots

'((X) (+ (*/ X X X `(* 1.0 1.0)) (*/ -3 X X 1.0) (* 2 X)))

Output:

<pre>-> ("0.000" "1.000" "2.000")</pre>

 

=={{header|PL/I}}==

f: procedure (x) returns (float (18));

declare x float (18);

end;

end locate_root;

 

=={{header|PureBasic}}==

{{trans|C++}}

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

EndProcedure

EndProcedure

 

 

=={{header|Python}}==

{{trans|Perl}}

 

step = 0.001 # Smaller step values produce more accurate and precise results

sign = value > 0

 

 

=={{header|R}}==

{{trans|Octave}}

 

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

}

 

 

=={{header|Racket}}==

 

#lang racket

 

 

(define tolerance (make-parameter 5e-16))

 

Different root-finding methods

 

(define (secant f a b)

(let next ([x1 a] [y1 (f a)] [x2 b] [y2 (f b)] [n 50])

c

(or (divide a c) (divide c b)))))))

 

Examples:

-> (find-root (λ (x) (- 2. (* x x))) 1 2)

1.414213562373095

-> (find-roots f -3 4 #:divisions 50)

'(2.4932181969624796e-33 1.0 2.0)

 

In order to provide a comprehensive code the given solution does not optimize the number of function calls.

 

Simple memoization operator

(define (memoized f)

(define tbl (make-hash))

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

 

The profiling shows that memoization reduces the number of function calls

(formerly Perl 6)

Uses exact arithmetic.

 

my $start = -1;

NEXT $sign = $next;

}

{{out}}

<pre>Root found at 0

Both of these REXX versions use the &nbsp; '''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.*/

f: parse arg x; return x * (x * (x-3) +2) /*formula used ──► x^3 - 3x^2 + 2x */

/*with factoring ──► x{ x^2 -3x + 2 } */

{{out|output|text=&nbsp; when using the defaults for input:}}

<pre>

===function coded in-line===

This version is about &nbsp; '''40%''' &nbsp; faster than the 1^st 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.*/

else if !\==$ then if !\==0 then say 'passed a root at' x/1

!= $ /*use the new sign for the next compare*/

{{out|output|text=&nbsp; is the same as the 1^st REXX version.}} <br><br>

 

=={{header|Ring}}==

load "stdlib.ring"

function = "return pow(x,3)-3*pow(x,2)+2*x"

oldsign = num

next

Output:

<pre>

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,

f = function(x)

{

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

0

 

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

{{trans|Python}}

 

x <=> 0

end

 

f = lambda { |x| x**3 - 3*x**2 + 2*x }

 

{{out}}

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

 

def sign

self <=> 0

end

 

 

=={{header|Rust}}==

 

use roots::find_roots_cubic;

 

println!("Result : {:?}", roots);

{{out}}

<pre>

 

Another without external crates:

use num::Float;

 

}

 

{{out}}

<pre>

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

 

printRoots(poly, -1.0, 4, 0.002)

 

===Functional version (Recommended)===

def findRoots(fn: Double => Double, start: Double, stop: Double, step: Double, epsilon: Double) = {

for {

 

println(findRoots(fn, -1.0, 3.0, 0.0001, 0.000000001))

{{out}}

Vector(-9.381755897326649E-14, 0.9999999999998124, 1.9999999999997022)

 

=={{header|Sidef}}==

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

}

}

sign = value>0

{{out}}

<pre>Root found at 0

=={{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]]

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

 

Result and timing:

<pre>/Tcl $ time ./froots.tcl

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

puts [frootsNR \

{x {expr {$x**3 - 3*$x**2 + 2*$x}}} \

 

=={{header|TI-89 BASIC}}==

{{trans|Go}}

{{libheader|Wren-fmt}}

 

var secant = Fn.new { |f, x0, x1|

 

var example = Fn.new { |x| x*x*x - 3*x*x + 2*x }

 

{{out}}

<pre>

=={{header|zkl}}==

{{trans|Haskell}}

[start..stop,step].filter('wrap(x){ f(x).closeTo(0.0,eps) })

{{out}}

<pre>L(-9.38176e-14,1,2)</pre>

{{trans|C}}

reg e=1.0e-12;

 

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

{{out}}

<pre>L(-0.032,0.968,1.068,1.968) --> L(1.87115e-19,1,1,2)</pre>