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)

 

=={{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|FreeBASIC}}==

Simple bisection method.

const iterations=20000000

 

next n

end if

{{out}}

<pre>in range -20 to 20

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

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

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