Anonymous user
Roots of a function: Difference between revisions
added D code
(Added Perl version) |
(added D code) |
||
Line 41:
}
</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}}==
<pre>sub f
| |||