A
adrianv
Tue, May 18, 2021 2:38 PM
I see that the original poster said the input should be monotonic, but the
bisection method does not require a monotonic input. It requires that the
function be continuous and that the function values at the interval
endpoints have opposite sign. A properly written implementation would check
for this and fail immediately if this condition doesn't hold.
The secant slope is not a good approximation to the derivative where by
"good" I mean that we can bound its error. If you want the derivative at
zero and use the interval [0,1] and consider a function like ln(x+c) then
the secant approximation is ln(1+c)-ln(c) and the actual derivative is 1/c.
This approximation is arbitrarily bad, in the limit as c->0. If c = 0.1
then the derivative is 10 and the secant approximation is 2.4. If c = 0.01
the derivative is 100 and the secant approximation is 4.6.
What you have implemented is actually the Secant Method:
https://en.wikipedia.org/wiki/Secant_method and the referenced page notes
that this method (like Newton's method) is not guaranteed to converge. This
is in contrast to the bisection method, which is guaranteed to converge,
albeit with only linear convergence. The rate of convergence of the secant
method when it does converge (to a simple root of a C2 function) is
superlinear, but not quadratic like Newton's method.
Unlike the original method, the solution need not lie in the given interval.
You're really initializing this method with two guesses for the root, not an
interval that contains the root. So where is the function supposed to be
monotonic? Over all the real numbers? If you try to solve for ln(x)=-2.3
(monotonic) on the interval [.001,1] it goes into an infinite loop with nans
because the iteration takes it outside the valid bound for the log function.
If you ask for exp(x)=0.5 and start it on [0,1] it correctly finds -0.69 as
the solution---not in the starting interval. If you give it atan you can
get it to diverge by giving it an overly large interval, even though atan is
monotonic. The requirement of monotonicity doesn't suffice to force good
behavior. If you try to solve sin(x)=0 and give it the interval [89,90]
where sin is monotonic (but has no solutions) you get the solution -133920
(which is correct, being 372*360). I suppose you could keep track of the
interval and return failure if the final solution is not in the interval.
Is it possible that the algorithm finds a root outside the starting interval
even though one exists there? (I think the answer is yes for the
non-monotonic case.)
Colin49 wrote
Hi,
In the original email it was stated that it was used for a monotonic
function. The same restriction applies to this.
The derivative is (HighArg-LowArg)/(FF(HighArg)-FF(LowArg)) or the inverse
of this. It's an approximation for the Hi/Lo range and it's not terrible.
It should work fine for any non constant monotonic function where the
values of FF at High & Low are not infinite.
The binary chop/search will not work if the function is not monotonic
and always
returns HighArg for a constant function.
For a constant function the results are undefined. It should fault with a
divide by zero. It would be easy to add a check for this. I can see there
is a possible issue if the function has a region of constant values
(derivative = 0) and we hit that with our Hi&Lo points.
In the original algorithm the user had to choose a High/Low range that was
meant to contain the solution.
I also think that it works better starting with a large range but my code
doesn't require that. With functions like Sine that are not monotonic over
their entire range it's important that the hi and low cover a monotonic
region and that the target value is inside the region. Much the same as
the
original binary chop algorithm.
I can see there is a possible issue if the function has a region of
constant values (derivative = 0) Just adding a test for HiValue ==
LowValue and returning (HighArg+LowArg)/2 would solve this. But again the
assumption was monotonically increasing function.
On Tue, May 18, 2021 at 10:09 AM adrianv <
Note that this isn't really Newton's method, which requires use of the
actual derivative, or even a standard approximate Newton's method, where
the derivative would be approximated numerically, e.g. by
(f(x+h)-f(x-h))/2h. Anybody planning to use it should be aware that the
solution can be any real number and is by no means guaranteed to lie in
the
interval [LowArg,HighArg]. Is there a convergence proof for this
algorithm? Can you prove that the interval shrinks? I have this feeling
that it might be possible to construct some kind of infinite loop case
where it fails to converge because the derivative approximation is so
horrible. It will blow up with an error if the value of the function is
equal at both endpoints, and it seems likely to misbehave if the function
values are almost equal at the endpoints.
Given these various issues, one might ask, why not just use a regular
approximate Newton's method, which will be better behaved, and easier to
understand, and wouldn't require the mysterious use of two parameters
(the
high and low values) to kick off the iteration. If you invoke this
with
the "high" and "low" values different by only a small amount like 1e-5
then
I think might get an actual approximate Newton method. (I'm haven't
tried
to prove that the interval doesn't grow, which would be necessary to
ensure
that the derivative approximation stays good.)
Colin49 wrote
Accuracy isn't req'd but I forgot to remove it from the recursive call of
InvFF. It still worked and I missed the warning message. Sorry, it was
very
late here.
There is now no accuracy setting and it still does less iterations than
the
Binary Chop. Put and echo(MidArg) after the let() to watch it in action.
function InvFF(Value,LowArg,HighArg)
= let(MidArg=LowArg +
(Value-FF(LowArg))*(HighArg-LowArg)/(FF(HighArg)-FF(LowArg)))
MidArg==LowArg || MidArg==HighArg
? MidArg
: InvFF(Value,HighArg,MidArg);
On Tue, May 18, 2021 at 3:06 AM LenStruttmann <[hidden email]
<http:///user/SendEmail.jtp?type=email&email=LenStruttmann%40>>
wrote:
Ummmm... "...even simpler.." is broke.
What happened to Accuracy parameter?
Sent from the OpenSCAD mailing list archive
<http://forum.openscad.org/>
at Nabble.com.
OpenSCAD mailing list
To unsubscribe send an email to [hidden email]
<http:///user/SendEmail.jtp?type=email&email=discuss-leave%40.openscad>
OpenSCAD mailing list
To unsubscribe send an email to [hidden email]
<http:///user/SendEmail.jtp?type=email&email=discuss-leave%40.openscad>
Sent from the OpenSCAD mailing list archive
<http://forum.openscad.org/>
at Nabble.com.
OpenSCAD mailing list
To unsubscribe send an email to
OpenSCAD mailing list
To unsubscribe send an email to
I see that the original poster said the input should be monotonic, but the
bisection method does not require a monotonic input. It requires that the
function be continuous and that the function values at the interval
endpoints have opposite sign. A properly written implementation would check
for this and fail immediately if this condition doesn't hold.
The secant slope is not a good approximation to the derivative where by
"good" I mean that we can bound its error. If you want the derivative at
zero and use the interval [0,1] and consider a function like ln(x+c) then
the secant approximation is ln(1+c)-ln(c) and the actual derivative is 1/c.
This approximation is arbitrarily bad, in the limit as c->0. If c = 0.1
then the derivative is 10 and the secant approximation is 2.4. If c = 0.01
the derivative is 100 and the secant approximation is 4.6.
What you have implemented is actually the Secant Method:
https://en.wikipedia.org/wiki/Secant_method and the referenced page notes
that this method (like Newton's method) is not guaranteed to converge. This
is in contrast to the bisection method, which *is* guaranteed to converge,
albeit with only linear convergence. The rate of convergence of the secant
method when it does converge (to a simple root of a C2 function) is
superlinear, but not quadratic like Newton's method.
Unlike the original method, the solution need not lie in the given interval.
You're really initializing this method with two guesses for the root, not an
interval that contains the root. So where is the function supposed to be
monotonic? Over all the real numbers? If you try to solve for ln(x)=-2.3
(monotonic) on the interval [.001,1] it goes into an infinite loop with nans
because the iteration takes it outside the valid bound for the log function.
If you ask for exp(x)=0.5 and start it on [0,1] it correctly finds -0.69 as
the solution---not in the starting interval. If you give it atan you can
get it to diverge by giving it an overly large interval, even though atan is
monotonic. The requirement of monotonicity doesn't suffice to force good
behavior. If you try to solve sin(x)=0 and give it the interval [89,90]
where sin is monotonic (but has no solutions) you get the solution -133920
(which is correct, being 372*360). I suppose you could keep track of the
interval and return failure if the final solution is not in the interval.
Is it possible that the algorithm finds a root outside the starting interval
even though one exists there? (I think the answer is yes for the
non-monotonic case.)
Colin49 wrote
> Hi,
>
> In the original email it was stated that it was used for a monotonic
> function. The same restriction applies to this.
>
> The derivative is (HighArg-LowArg)/(FF(HighArg)-FF(LowArg)) or the inverse
> of this. It's an approximation for the Hi/Lo range and it's not terrible.
>
> It should work fine for any non constant monotonic function where the
> values of FF at High & Low are not infinite.
>
> The binary chop/search will not work if the function is not monotonic
> and always
> returns HighArg for a constant function.
>
> For a constant function the results are undefined. It should fault with a
> divide by zero. It would be easy to add a check for this. I can see there
> is a possible issue if the function has a region of constant values
> (derivative = 0) and we hit that with our Hi&Lo points.
>
> In the original algorithm the user had to choose a High/Low range that was
> meant to contain the solution.
>
> I also think that it works better starting with a large range but my code
> doesn't require that. With functions like Sine that are not monotonic over
> their entire range it's important that the hi and low cover a monotonic
> region and that the target value is inside the region. Much the same as
> the
> original binary chop algorithm.
>
> I can see there is a possible issue if the function has a region of
> constant values (derivative = 0) Just adding a test for HiValue ==
> LowValue and returning (HighArg+LowArg)/2 would solve this. But again the
> assumption was monotonically increasing function.
>
> On Tue, May 18, 2021 at 10:09 AM adrianv <
> avm4@
> > wrote:
>
>> Note that this isn't really Newton's method, which requires use of the
>> actual derivative, or even a standard approximate Newton's method, where
>> the derivative would be approximated numerically, e.g. by
>> (f(x+h)-f(x-h))/2h. Anybody planning to use it should be aware that the
>> solution can be any real number and is by no means guaranteed to lie in
>> the
>> interval [LowArg,HighArg]. Is there a convergence proof for this
>> algorithm? Can you prove that the interval shrinks? I have this feeling
>> that it might be possible to construct some kind of infinite loop case
>> where it fails to converge because the derivative approximation is so
>> horrible. It will blow up with an error if the value of the function is
>> equal at both endpoints, and it seems likely to misbehave if the function
>> values are almost equal at the endpoints.
>>
>> Given these various issues, one might ask, why not just use a regular
>> approximate Newton's method, which will be better behaved, and easier to
>> understand, and wouldn't require the mysterious use of two parameters
>> (the
>> high and low values) to kick off the iteration. If you invoke this
>> with
>> the "high" and "low" values different by only a small amount like 1e-5
>> then
>> I think might get an actual approximate Newton method. (I'm haven't
>> tried
>> to prove that the interval doesn't grow, which would be necessary to
>> ensure
>> that the derivative approximation stays good.)
>>
>> Colin49 wrote
>> Accuracy isn't req'd but I forgot to remove it from the recursive call of
>> InvFF. It still worked and I missed the warning message. Sorry, it was
>> very
>> late here.
>>
>> There is now no accuracy setting and it still does less iterations than
>> the
>> Binary Chop. Put and echo(MidArg) after the let() to watch it in action.
>>
>> function InvFF(Value,LowArg,HighArg)
>> = let(MidArg=LowArg +
>> (Value-FF(LowArg))*(HighArg-LowArg)/(FF(HighArg)-FF(LowArg)))
>> MidArg==LowArg || MidArg==HighArg
>> ? MidArg
>> : InvFF(Value,HighArg,MidArg);
>>
>>
>> On Tue, May 18, 2021 at 3:06 AM LenStruttmann <[hidden email]
>> <http:///user/SendEmail.jtp?type=email&email=LenStruttmann%40>>
>> wrote:
>>
>> > Ummmm... "...even simpler.." is broke.
>> >
>> > What happened to Accuracy parameter?
>> > ------------------------------
>> > Sent from the OpenSCAD mailing list archive
>> <http://forum.openscad.org/>
>> > at Nabble.com.
>> > _______________________________________________
>> > OpenSCAD mailing list
>> > To unsubscribe send an email to [hidden email]
>> <http:///user/SendEmail.jtp?type=email&email=discuss-leave%40.openscad>
>> >
>>
>> _______________________________________________
>> OpenSCAD mailing list
>> To unsubscribe send an email to [hidden email]
>> <http:///user/SendEmail.jtp?type=email&email=discuss-leave%40.openscad>
>>
>>
>> ------------------------------
>> Sent from the OpenSCAD mailing list archive
>> <http://forum.openscad.org/>
>> at Nabble.com.
>> _______________________________________________
>> OpenSCAD mailing list
>> To unsubscribe send an email to
> discuss-leave@.openscad
>>
>
> _______________________________________________
> OpenSCAD mailing list
> To unsubscribe send an email to
> discuss-leave@.openscad
--
Sent from: http://forum.openscad.org/
A
adrianv
Thu, May 20, 2021 9:21 PM
I got curious about the right way to solve the root finding problem and did a
little investigation. So here's a stab at it. This version is a
bisection/inverse quadratic interpolation approach, which seems to be the
standard method for root finding. That is, you use an inverse quadratic
step if it works and if it doesn't you use bisection. This gives a "best of
both worlds" approach, where the algorithm is robust, but also gives fast
convergence. I'm not sure I like the convergence criterion, though. At
least, it can give unexpectedly poor results for functions with a large
range unless you really set the tolerance small. Anyway, give it a shot
and let me know if you can break it.
// Function: rootfind()
// Usage:
// x = rootfind(f, x0, x1, [tol])
// Description:
// Find a root of the continuous function f where the sign of f(x0) is
different
// from the sign of f(x1). The tolerance (tol) specifies the accuracy of
the solution:
// abs(f(x)) < tol * yrange, where yrange is the range of observed
function values.
// This function can only find roots that cross the x axis: it cannot
find the
// the root of x^2.
function rootfind(f,x0,x1,tol=1e-15) =
let(
y0 = f(x0),
y1 = f(x1),
yrange = y0<y1 ? [y0,y1] : [y1,y0]
)
// Check endpoints
y0==0 || _rfcheck(y0,yrange,tol) ? x0 :
y1==0 || _rfcheck(y1,yrange,tol) ? x1 :
assert(y0*y1<0, "Sign of function must be different at the interval
endpoints")
_rootfind(f,[x0,x1],[y0,y1],yrange,tol);
function _rfcheck(y,range,tol) = abs(y) < tol*(range[1]-range[0]);
function _rootfind(f, xpts, ypts, yrange, tol, i=0) =
//echo(i, x=[xpts[0],xpts[1]], y=[ypts[0],ypts[1]]) // uncomment to
see steps
assert(i<100, "rootfind did not converge to a solution")
let(
xmid = (xpts[0]+xpts[1])/2,
ymid = f(xmid),
yrange = [min(ymid, yrange[0]), max(ymid, yrange[1])]
)
_rfcheck(ymid, yrange, tol) ? xmid :
let(
// Force root to be between x0 and midpoint
y = ymid * ypts[0] < 0 ? [ypts[0], ymid, ypts[1]]
: [ypts[1], ymid, ypts[0]],
x = ymid * ypts[0] < 0 ? [xpts[0], xmid, xpts[1]]
: [xpts[1], xmid, xpts[0]],
v = y[2](y[2]-y[0]) - 2y[1](y[1]-y[0])
)
v <= 0 ? _rootfind(f,x,y,yrange,tol,i+1) // Root is between first two
points, extra 3rd point doesn't hurt
:
let( // Do quadratic approximation
B = (x[1]-x[0]) / (y[1]-y[0]),
C = y[-1,2,-1] / (y[2]-y[1]) / (y[2]-y[0]),
newx = x[0] - B * y[0] (1-Cy[1]),
newy = f(newx),
new_yrange = [min(yrange[0],newy), max(yrange[1], newy)],
// select interval that contains the root by checking sign
yinterval = newyy[0] < 0 ? [y[0],newy] : [newy,y[1]],
xinterval = newyy[0] < 0 ? [x[0],newx] : [newx,x[1]]
)
_rfcheck(newy, new_yrange, tol)
? newx
: _rootfind(f, xinterval, yinterval, new_yrange, tol, i+1);
--
Sent from: http://forum.openscad.org/
I got curious about the right way to solve the root finding problem and did a
little investigation. So here's a stab at it. This version is a
bisection/inverse quadratic interpolation approach, which seems to be the
standard method for root finding. That is, you use an inverse quadratic
step if it works and if it doesn't you use bisection. This gives a "best of
both worlds" approach, where the algorithm is robust, but also gives fast
convergence. I'm not sure I like the convergence criterion, though. At
least, it can give unexpectedly poor results for functions with a large
range unless you really set the tolerance small. Anyway, give it a shot
and let me know if you can break it.
// Function: rootfind()
// Usage:
// x = rootfind(f, x0, x1, [tol])
// Description:
// Find a root of the continuous function f where the sign of f(x0) is
different
// from the sign of f(x1). The tolerance (tol) specifies the accuracy of
the solution:
// abs(f(x)) < tol * yrange, where yrange is the range of observed
function values.
// This function can only find roots that cross the x axis: it cannot
find the
// the root of x^2.
function rootfind(f,x0,x1,tol=1e-15) =
let(
y0 = f(x0),
y1 = f(x1),
yrange = y0<y1 ? [y0,y1] : [y1,y0]
)
// Check endpoints
y0==0 || _rfcheck(y0,yrange,tol) ? x0 :
y1==0 || _rfcheck(y1,yrange,tol) ? x1 :
assert(y0*y1<0, "Sign of function must be different at the interval
endpoints")
_rootfind(f,[x0,x1],[y0,y1],yrange,tol);
function _rfcheck(y,range,tol) = abs(y) < tol*(range[1]-range[0]);
function _rootfind(f, xpts, ypts, yrange, tol, i=0) =
//echo(i, x=[xpts[0],xpts[1]], y=[ypts[0],ypts[1]]) // uncomment to
see steps
assert(i<100, "rootfind did not converge to a solution")
let(
xmid = (xpts[0]+xpts[1])/2,
ymid = f(xmid),
yrange = [min(ymid, yrange[0]), max(ymid, yrange[1])]
)
_rfcheck(ymid, yrange, tol) ? xmid :
let(
// Force root to be between x0 and midpoint
y = ymid * ypts[0] < 0 ? [ypts[0], ymid, ypts[1]]
: [ypts[1], ymid, ypts[0]],
x = ymid * ypts[0] < 0 ? [xpts[0], xmid, xpts[1]]
: [xpts[1], xmid, xpts[0]],
v = y[2]*(y[2]-y[0]) - 2*y[1]*(y[1]-y[0])
)
v <= 0 ? _rootfind(f,x,y,yrange,tol,i+1) // Root is between first two
points, extra 3rd point doesn't hurt
:
let( // Do quadratic approximation
B = (x[1]-x[0]) / (y[1]-y[0]),
C = y*[-1,2,-1] / (y[2]-y[1]) / (y[2]-y[0]),
newx = x[0] - B * y[0] *(1-C*y[1]),
newy = f(newx),
new_yrange = [min(yrange[0],newy), max(yrange[1], newy)],
// select interval that contains the root by checking sign
yinterval = newy*y[0] < 0 ? [y[0],newy] : [newy,y[1]],
xinterval = newy*y[0] < 0 ? [x[0],newx] : [newx,x[1]]
)
_rfcheck(newy, new_yrange, tol)
? newx
: _rootfind(f, xinterval, yinterval, new_yrange, tol, i+1);
--
Sent from: http://forum.openscad.org/
