Roots of a function: Difference between revisions
Content added Content deleted
(Added Perl version) |
(added D code) |
||
Line 41: | Line 41: | ||
} |
} |
||
</pre> |
</pre> |
||
=={{header|D}}== |
|||
<pre>module findroot ; |
|||
import std.stdio ; |
|||
import std.math ; |
|||
void report(T)(T[] r, T function(T) f, T tolerance = cast(T) 1e-4L) { |
|||
if (r.length) { |
|||
writefln("Root found (tolerance = %1.4g) :", tolerance) ; |
|||
foreach(x ; r) { |
|||
T y = f(x) ; |
|||
if (nearZero(y)) |
|||
writefln("... EXACTLY at %+1.20f, f(x) = %+1.4g", x, y) ; |
|||
else if (nearZero(y, tolerance)) |
|||
writefln(".... MAY-BE at %+1.20f, f(x) = %+1.4g", x, y) ; |
|||
else |
|||
writefln("Verify needed, f(%1.4g) = %1.4g > tolerance in magnitude", x, y) ; |
|||
} |
|||
} else |
|||
writefln("No root found.") ; |
|||
} |
|||
T nearZero(T)(T a, T b = T.epsilon * 4) { return abs(a) <= b ; } |
|||
T[] findroot(T)(T function(T) f, T start, T end, T step = cast(T) 0.001L, |
|||
T tolerance = cast(T) 1e-4L) { |
|||
T[T] result ; |
|||
if (nearZero(step)) |
|||
writefln("WARNING: step size may be too small.") ; |
|||
T searchRoot(T a, T b) { // search root by simple bisection |
|||
T root ; |
|||
int limit = 49 ; |
|||
T gap = b - a ; |
|||
while (!nearZero(gap) && limit--) { |
|||
if (nearZero(f(a))) return a ; |
|||
if (nearZero(f(b))) return b ; |
|||
root = (b + a)/2.0L ; |
|||
if (nearZero(f(root))) return root ; |
|||
if (f(a) * f(root) < 0) |
|||
b = root ; |
|||
else |
|||
a = root ; |
|||
gap = b - a ; |
|||
} |
|||
return root ; |
|||
} |
|||
T dir = cast(T) (end > start ? 1.0 : -1.0) ; |
|||
step = (end > start) ? abs(step) : - abs(step) ; |
|||
for(T x = start ; x*dir <= end*dir ; x = x + step) |
|||
if (f(x)*f(x + step) <= 0) { |
|||
T r = searchRoot(x, x+ step) ; |
|||
result[r] = f(r) ; |
|||
} |
|||
return result.keys.sort ; // reduce duplacated root, if any |
|||
} |
|||
real f(real x){ |
|||
return x*x*x - 3*x*x + 2*x ; |
|||
} |
|||
void main(){ |
|||
findroot(&f, -1.0L, 3.0L, 0.001L).report(&f) ; |
|||
}</pre> |
|||
Output ( NB:smallest increment for real type in D is real.epsilon = 1.0842e-19 ): |
|||
<pre>Root found (tolerance = 0.0001) : |
|||
.... MAY-BE at -0.00000000000000000080, f(x) = -1.603e-18 |
|||
... EXACTLY at +1.00000000000000000020, f(x) = -2.168e-19 |
|||
.... MAY-BE at +1.99999999999999999950, f(x) = -8.674e-19</pre> |
|||
=={{header|Perl}}== |
=={{header|Perl}}== |
||
<pre>sub f |
<pre>sub f |