TLDR
Give calc's a R (calc-find-root)
[root] function
a range, instead of a single numerical guess. Make sure the function
evaluates to different signs on each end of the range. This gives calc
enough information to fall back to a different root finding algorithm that
will succeed.
For example, here I am attempting to solve utpn(x, 10, 1) = 0.5 for x, guessing that it might be x = 1.0.
` utpn(x,10,1)=0.5 RET
1.0 RET <---- the guess
a R RET
x RET
The above doesn't work, producing the error message "Newton's method failed to converge: 4.8639203627839e17".
Calc succeeds if the guess is the range [0..1000]:
` utpn(x, 10, 1) = 0.5 RET
[0..1000] RET <---- the better guess
a R RET
x RET
Success! Calc produces this:
2: utpn(x, 10, 1) = 0.5
1: [10., 0.]
Background
As pointed out in both NickD's
answer and dbx48's
answer, calc is having
difficulty due to the extremely flat portions of the utpn function. Let's graph utpn(x, 10, 1) for x over the range [0..20]:
[0..200]
10
/
' utpn(x, 10, 1)
g f
You can see that outside the range of about [7..13] the function is very flat.
The problem with Newton's algorithm
Clues into Calc's behavior are found in its Root
Finding
documentation.
The initial guess can be a real number, in which case Calc searches for a
real solution near that number, or a complex number, in which case Calc
searches the whole complex plane near that number for a solution, or it
can be an interval form which restricts the search to real numbers inside
that interval.
Recall that in my original question I was giving the a R function a
single guess:
` utpn(x, 10, 1) = 0.5 RET
1.0 RET <---- the guess
a R RET
x RET
On this the documentation says:
Calc tries to use a d to take the derivative of the equation. If this
succeeds, it uses Newton’s method. If the equation is not differentiable
Calc uses a bisection method.
So, let's see if utpn(x, 10, 1) - 0.5 is differentiable.
` utpn(x, 10, 1) - 0.5 RET
a d RET
x RET
Calc presents this derivative:
1: 0.398942280401 / exp((7.07106781188 - x / 1.41421356237)^2)
Let's plot that derivative over a range of [0 .. 20]:
[0..200]
10
/
' 0.398942280401 / exp((7.07106781188 - x / 1.41421356237)^2)
g f
The derivative is nearly zero outside of the range [7 .. 13].
Numerically, the derivative approaches zero asymptotically,
because utpn itself asymptotically approaches zero.
Because calc thinks it can compute the derivative it will use (Newton's
method)[https://en.wikipedia.org/wiki/Newton%27s_method) to find the root,
and here things go astray.
Why? There is nothing mathematically wrong with using Newton's method with utpn. My guess: calc's limited precision arithmetic are exceeded during the computation, which "confuses" the algorithm. One intuition: Newton's method computes the next guess for x as x - f(x)/f'(x). Above calc produced an error message "Newton's method failed to converge: 4.8639203627839e17". Let's see the values involved at that point:
utpn(4.8639203627839e17, 10, 1) - 0.5 evaluates to -0.50.398942280401 / exp((7.07106781188 - 4.8639203627839e17 / 1.41421356237)^2) yields an error: "Floating-point overflow occurred: exp(1.18288606478e35)"
I have no explanation for how Newton's algorithm has gone from the initial guess of 1.0 all the way out to 4.8639203627839e17 without realizing the root was at 10.0. Presumably calc's limited precision arithmetic is to blame for that, too.
Calc's fallback algorithm: Bisection
Calc's documentation claims it has code to
detect when Newton's method is failing:
(If Newton’s method appears to be going astray, Calc switches over to
bisection if it can, or otherwise gives up. In this case it may help to
try again with a slightly different initial guess.)
For some reason this is not working for the particular case discussed here.
Calc does have an alternate algorithm, which uses a bisection approach.
Read up on Bisection method on
Wikipedia. Because the
method does not use the derivative it is guaranteed to find a root if it
has an interval where f(x) is of different sign at each end. The calc
manual discusses this:
If the formula (or the difference between the sides of an equation) is
negative at one end of the interval you specify and positive at the other
end, the root finder is guaranteed to find a root. Otherwise, Calc
subdivides the interval into small parts looking for positive and
negative values to bracket the root. When your guess is an interval, Calc
will not look outside that interval for a root.
Helping Calc use the Bisection method
Let's try the original problem giving an interval guess instead of a
single number:
` utpn(x, 10, 1) = 0.5 RET
[0..1000] RET
a R RET
x RET
Success! Calc produces this:
2: utpn(x, 10, 1) = 0.5
1: [10., 0.]
Could wroot work too?
Calc has another root solving function, wroot, described in the manual as follows:
The H a R [wroot] command is similar to a R, except that if the initial
guess is an interval for which the function has the same sign at both
ends, then rather than subdividing the interval Calc attempts to widen it
to enclose a root. Use this mode if you are not sure if the function has
a root in your interval.
How might this fare with the original problematic case:
` utpn(x, 10, 1) = 0.5 RET
1.0 RET
H a R RET
x RET
No luck:
2: utpn(x, 10, 1) = 0.5
1: wroot(utpn(x, 10, 1) = 0.5, x, 1.)
With the message "Newton's method failed to converge: 4.8639203627839e17".
Presumably, what is going on here is that the extremely flat derivative
causes the root finding algorithm to quickly reach very large x values,
where the derivative becomes unstable.
Conclusion
The flat derivative of the utpn function is a worst case scenario for the
Newton's method
algorithm used by calc, given calc's limited precision arithmetic. Calc
can sometimes detect this and switch to a (Bisection
method)[https://en.wikipedia.org/wiki/Bisection_method], but in this case
that heuristic is failing. To work around this, supply an explicit range
to the "a R (calc-find-root)
[root]"
function, such that the sign of the result is different for each end of the
interval, and the Bisection method should kick in and work.
