ASKSAGE: Sage Q&A Forum - RSS feedhttps://ask.sagemath.org/questions/Q&A Forum for SageenCopyright Sage, 2010. Some rights reserved under creative commons license.Sun, 12 Feb 2017 19:49:49 +0100find_root does not fulfill tolerancehttps://ask.sagemath.org/question/36505/find_root-does-not-fulfill-tolerance/In Sage 7.3 the result of:
find_root((4*x^4-x^2).function(x),-0.25,0.4,xtol=1e-13,rtol=1e-13)
is
1.6619679287440101e-10
One would expect that the error with the exact solution (x=0) would be lower than 1e-13. I imagine that the multiplicity of the root produces this effect. Does somebody know the numerical algorithm behind find_root?
Thanks in advance.
**UPDATE:** Thanks to @kcrisman, I see that find_root calls `brentq`. The documentation of such function says:
*Safer algorithms are brentq, brenth, ridder, and bisect, but they **all require** that the root first be bracketed in an **interval where the function changes sign**. The brentq algorithm is recommended for general use in one dimensional problems when such an interval has been found.*
Since 4x^4-x^2 has a root in 0 with multiplicity 2, this condition is not fulfilled. Is there any other command in Sage different from `find_root` to numerically calculate roots for general functions that can handle with roots with multiplicity bigger than 1?Thu, 09 Feb 2017 17:44:53 +0100https://ask.sagemath.org/question/36505/find_root-does-not-fulfill-tolerance/Answer by kcrisman for <p>In Sage 7.3 the result of:</p>
<pre><code>find_root((4*x^4-x^2).function(x),-0.25,0.4,xtol=1e-13,rtol=1e-13)
</code></pre>
<p>is</p>
<pre><code>1.6619679287440101e-10
</code></pre>
<p>One would expect that the error with the exact solution (x=0) would be lower than 1e-13. I imagine that the multiplicity of the root produces this effect. Does somebody know the numerical algorithm behind find_root?</p>
<p>Thanks in advance. </p>
<p><strong>UPDATE:</strong> Thanks to <a href="/users/41/kcrisman/">@kcrisman</a>, I see that find_root calls <code>brentq</code>. The documentation of such function says:</p>
<p><em>Safer algorithms are brentq, brenth, ridder, and bisect, but they <strong>all require</strong> that the root first be bracketed in an <strong>interval where the function changes sign</strong>. The brentq algorithm is recommended for general use in one dimensional problems when such an interval has been found.</em></p>
<p>Since 4x^4-x^2 has a root in 0 with multiplicity 2, this condition is not fulfilled. Is there any other command in Sage different from <code>find_root</code> to numerically calculate roots for general functions that can handle with roots with multiplicity bigger than 1?</p>
https://ask.sagemath.org/question/36505/find_root-does-not-fulfill-tolerance/?answer=36508#post-id-36508If you do `find_root??` you will see that the final line of code is a Scipy algorithm:
return scipy.optimize.brentq(f, a, b,
full_output=full_output, xtol=xtol, rtol=rtol, maxiter=maxiter)
The [Scipy doc](https://docs.scipy.org/doc/scipy-0.18.1/reference/generated/scipy.optimize.brentq.html) indicates it should indeed return within the tolerance expected, so I'm not sure what is going on, but at least this is where you should start your search.Thu, 09 Feb 2017 19:07:26 +0100https://ask.sagemath.org/question/36505/find_root-does-not-fulfill-tolerance/?answer=36508#post-id-36508Comment by mforets for <p>If you do <code>find_root??</code> you will see that the final line of code is a Scipy algorithm:</p>
<pre><code>return scipy.optimize.brentq(f, a, b,
full_output=full_output, xtol=xtol, rtol=rtol, maxiter=maxiter)
</code></pre>
<p>The <a href="https://docs.scipy.org/doc/scipy-0.18.1/reference/generated/scipy.optimize.brentq.html">Scipy doc</a> indicates it should indeed return within the tolerance expected, so I'm not sure what is going on, but at least this is where you should start your search.</p>
https://ask.sagemath.org/question/36505/find_root-does-not-fulfill-tolerance/?comment=36509#post-id-36509`brentq` requires a sign change of $f$ in the given interval, but this cannot be achieved in $[-0.25, 0.4]$ (note that $f$ is tangent to $y=0$ at $0$). Maybe that's the reason?Thu, 09 Feb 2017 20:08:04 +0100https://ask.sagemath.org/question/36505/find_root-does-not-fulfill-tolerance/?comment=36509#post-id-36509Comment by kcrisman for <p>If you do <code>find_root??</code> you will see that the final line of code is a Scipy algorithm:</p>
<pre><code>return scipy.optimize.brentq(f, a, b,
full_output=full_output, xtol=xtol, rtol=rtol, maxiter=maxiter)
</code></pre>
<p>The <a href="https://docs.scipy.org/doc/scipy-0.18.1/reference/generated/scipy.optimize.brentq.html">Scipy doc</a> indicates it should indeed return within the tolerance expected, so I'm not sure what is going on, but at least this is where you should start your search.</p>
https://ask.sagemath.org/question/36505/find_root-does-not-fulfill-tolerance/?comment=36516#post-id-36516Maybe - I see what you mean here. But some experimentation with some perturbations of your example don't give me confidence, since it doesn't satisfy its own documentation - http://sagecell.sagemath.org/?z=eJxtjTEKxCAQAPtA_rClHiom8Uq_kkMSA4LuitmA9_tLd01gimmGSaVSY8Cr1C-EE7COQwAPR8L904hYCPfqq9N9naU5Ltw4EYoulbZmfitrnOpM2U9RT4tqfy2hJ47NL9bKcagt4b2pJuS8ZTqjsArCzVMrfwpBLlQ=&lang=sageThu, 09 Feb 2017 21:36:28 +0100https://ask.sagemath.org/question/36505/find_root-does-not-fulfill-tolerance/?comment=36516#post-id-36516Comment by mforets for <p>If you do <code>find_root??</code> you will see that the final line of code is a Scipy algorithm:</p>
<pre><code>return scipy.optimize.brentq(f, a, b,
full_output=full_output, xtol=xtol, rtol=rtol, maxiter=maxiter)
</code></pre>
<p>The <a href="https://docs.scipy.org/doc/scipy-0.18.1/reference/generated/scipy.optimize.brentq.html">Scipy doc</a> indicates it should indeed return within the tolerance expected, so I'm not sure what is going on, but at least this is where you should start your search.</p>
https://ask.sagemath.org/question/36505/find_root-does-not-fulfill-tolerance/?comment=36522#post-id-36522..still i think this is a border case, because the root cannot be properly "enclosed". For completeness: it is also possible to do `R.<x> = RDF[]`, `f = 4*x^4 - x^2`, `f.roots()`.Thu, 09 Feb 2017 22:59:37 +0100https://ask.sagemath.org/question/36505/find_root-does-not-fulfill-tolerance/?comment=36522#post-id-36522Comment by kcrisman for <p>If you do <code>find_root??</code> you will see that the final line of code is a Scipy algorithm:</p>
<pre><code>return scipy.optimize.brentq(f, a, b,
full_output=full_output, xtol=xtol, rtol=rtol, maxiter=maxiter)
</code></pre>
<p>The <a href="https://docs.scipy.org/doc/scipy-0.18.1/reference/generated/scipy.optimize.brentq.html">Scipy doc</a> indicates it should indeed return within the tolerance expected, so I'm not sure what is going on, but at least this is where you should start your search.</p>
https://ask.sagemath.org/question/36505/find_root-does-not-fulfill-tolerance/?comment=36524#post-id-36524Hmm, but RDF doesn't even do more than a certain tolerance. Maybe it only supports tol up to a certain point?Thu, 09 Feb 2017 23:53:34 +0100https://ask.sagemath.org/question/36505/find_root-does-not-fulfill-tolerance/?comment=36524#post-id-36524Comment by mforets for <p>If you do <code>find_root??</code> you will see that the final line of code is a Scipy algorithm:</p>
<pre><code>return scipy.optimize.brentq(f, a, b,
full_output=full_output, xtol=xtol, rtol=rtol, maxiter=maxiter)
</code></pre>
<p>The <a href="https://docs.scipy.org/doc/scipy-0.18.1/reference/generated/scipy.optimize.brentq.html">Scipy doc</a> indicates it should indeed return within the tolerance expected, so I'm not sure what is going on, but at least this is where you should start your search.</p>
https://ask.sagemath.org/question/36505/find_root-does-not-fulfill-tolerance/?comment=36587#post-id-36587@kcrisman Yes, I think that's the case; some experiments with the function `4*x^4-x^2+0.01` around the exact root `-1/2*sqrt(-1/10*sqrt(21) + 1/2)` (this time in $[-0.25, 0.0]$) show it works ok until `rtol` is too small, giving `ValueError: rtol too small (1e-16 < 4.44089e-16)`. But alternatively for the OP's question, let us now that using [Newton's method](https://docs.scipy.org/doc/scipy/reference/generated/scipy.optimize.newton.html), as in `newton(4*x^4-x^2, -0.25, tol=1e-14, maxiter=3000)`, it will pass the tests with `1e-13` (and smaller).Sun, 12 Feb 2017 19:49:49 +0100https://ask.sagemath.org/question/36505/find_root-does-not-fulfill-tolerance/?comment=36587#post-id-36587