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.Fri, 09 Aug 2013 19:21:48 +0200Finding root with difficultieshttps://ask.sagemath.org/question/10429/finding-root-with-difficulties/I'm working with this example:
sage: var('x a t', domain='real')
sage: assume(a > 0, a < 1)
sage: xt(t) = solve(integrate(-4*sqrt(17)/(9*x^2 + 6*x - 16),
... (x, 0, a)).subs(a = x) == t, x)[2].rhs()
sage: find_root(xt - 1, 2, 3)
0.0
sage: plot(xt - 1, (t, 2, 3))
verbose 0 (2424: plot.py, generate_plot_points) WARNING: When plotting, failed to evaluate function at 137 points.
verbose 0 (2424: plot.py, generate_plot_points) Last error message: ''
The plot reports errors, but it *looks* reasonable enough. The existance of a root appears very clear, but sage is unable to find it. Furthermore, sage does not even properly report its inability to find a root, but instead returns a value as if there were no problems. The fact that this value lies outside the interval is an indication that there might be something fishy happening here.
Is there any way to get at my root numerically, or perhaps even symbolically? I know that the function `xt` is an ugly beast, but some form of robust bisection should be possible, I think.Fri, 09 Aug 2013 03:46:22 +0200https://ask.sagemath.org/question/10429/finding-root-with-difficulties/Answer by Eviatar Bach for <p>I'm working with this example:</p>
<pre><code>sage: var('x a t', domain='real')
sage: assume(a > 0, a < 1)
sage: xt(t) = solve(integrate(-4*sqrt(17)/(9*x^2 + 6*x - 16),
... (x, 0, a)).subs(a = x) == t, x)[2].rhs()
sage: find_root(xt - 1, 2, 3)
0.0
sage: plot(xt - 1, (t, 2, 3))
verbose 0 (2424: plot.py, generate_plot_points) WARNING: When plotting, failed to evaluate function at 137 points.
verbose 0 (2424: plot.py, generate_plot_points) Last error message: ''
</code></pre>
<p>The plot reports errors, but it <em>looks</em> reasonable enough. The existance of a root appears very clear, but sage is unable to find it. Furthermore, sage does not even properly report its inability to find a root, but instead returns a value as if there were no problems. The fact that this value lies outside the interval is an indication that there might be something fishy happening here.</p>
<p>Is there any way to get at my root numerically, or perhaps even symbolically? I know that the function <code>xt</code> is an ugly beast, but some form of robust bisection should be possible, I think.</p>
https://ask.sagemath.org/question/10429/finding-root-with-difficulties/?answer=15331#post-id-15331Unfortunately, the current root-finding only works with real numbers (see [#14991](http://trac.sagemath.org/ticket/14991)), and `find_root` is encountering complex numbers.
With mpmath's root-finding you get the following answer, although the syntax is extremely clumsy at the moment:
from sage.libs.mpmath.utils import call
import mpmath
def mpmath_xt(t, prec):
t = ComplexField(prec)(t.real, t.imag)
ft = xt(t).n(prec)
return mpmath.mpc(ft.real(), ft.imag()) - 1
Then,
sage: call(mpmath.findroot, lambda t: mpmath_xt(t, 100), [2, 3], parent=ComplexField(100))
2.4629947809517837884000550617 + 2.0631855558523985234600213842e-88*I
Change 100 in the above to whatever precision you want.
Hope this helps!Fri, 09 Aug 2013 05:16:08 +0200https://ask.sagemath.org/question/10429/finding-root-with-difficulties/?answer=15331#post-id-15331Comment by Eviatar Bach for <p>Unfortunately, the current root-finding only works with real numbers (see <a href="http://trac.sagemath.org/ticket/14991">#14991</a>), and <code>find_root</code> is encountering complex numbers.</p>
<p>With mpmath's root-finding you get the following answer, although the syntax is extremely clumsy at the moment:</p>
<pre><code>from sage.libs.mpmath.utils import call
import mpmath
def mpmath_xt(t, prec):
t = ComplexField(prec)(t.real, t.imag)
ft = xt(t).n(prec)
return mpmath.mpc(ft.real(), ft.imag()) - 1
</code></pre>
<p>Then,</p>
<pre><code>sage: call(mpmath.findroot, lambda t: mpmath_xt(t, 100), [2, 3], parent=ComplexField(100))
2.4629947809517837884000550617 + 2.0631855558523985234600213842e-88*I
</code></pre>
<p>Change 100 in the above to whatever precision you want.</p>
<p>Hope this helps!</p>
https://ask.sagemath.org/question/10429/finding-root-with-difficulties/?comment=17169#post-id-17169Okay, I looked into it, and I think you're right; it's not the wrong way to do it. The problem is precisely that it is converting the expression to a `fast_float`, but then chokes when it encounters complex numbers.
Look what happens:
sage: (xt - 1)._fast_float_(t)(1)
nanFri, 09 Aug 2013 19:17:05 +0200https://ask.sagemath.org/question/10429/finding-root-with-difficulties/?comment=17169#post-id-17169Comment by Eviatar Bach for <p>Unfortunately, the current root-finding only works with real numbers (see <a href="http://trac.sagemath.org/ticket/14991">#14991</a>), and <code>find_root</code> is encountering complex numbers.</p>
<p>With mpmath's root-finding you get the following answer, although the syntax is extremely clumsy at the moment:</p>
<pre><code>from sage.libs.mpmath.utils import call
import mpmath
def mpmath_xt(t, prec):
t = ComplexField(prec)(t.real, t.imag)
ft = xt(t).n(prec)
return mpmath.mpc(ft.real(), ft.imag()) - 1
</code></pre>
<p>Then,</p>
<pre><code>sage: call(mpmath.findroot, lambda t: mpmath_xt(t, 100), [2, 3], parent=ComplexField(100))
2.4629947809517837884000550617 + 2.0631855558523985234600213842e-88*I
</code></pre>
<p>Change 100 in the above to whatever precision you want.</p>
<p>Hope this helps!</p>
https://ask.sagemath.org/question/10429/finding-root-with-difficulties/?comment=17168#post-id-17168I should also mention that you need to be careful with the precision; your expression is so big that there's probably roundoff error, and the answer is not correct to the precision you specify. So you should work at a higher precision than the one you want and then discard some of the trailing digits.Fri, 09 Aug 2013 19:21:48 +0200https://ask.sagemath.org/question/10429/finding-root-with-difficulties/?comment=17168#post-id-17168Comment by MvG for <p>Unfortunately, the current root-finding only works with real numbers (see <a href="http://trac.sagemath.org/ticket/14991">#14991</a>), and <code>find_root</code> is encountering complex numbers.</p>
<p>With mpmath's root-finding you get the following answer, although the syntax is extremely clumsy at the moment:</p>
<pre><code>from sage.libs.mpmath.utils import call
import mpmath
def mpmath_xt(t, prec):
t = ComplexField(prec)(t.real, t.imag)
ft = xt(t).n(prec)
return mpmath.mpc(ft.real(), ft.imag()) - 1
</code></pre>
<p>Then,</p>
<pre><code>sage: call(mpmath.findroot, lambda t: mpmath_xt(t, 100), [2, 3], parent=ComplexField(100))
2.4629947809517837884000550617 + 2.0631855558523985234600213842e-88*I
</code></pre>
<p>Change 100 in the above to whatever precision you want.</p>
<p>Hope this helps!</p>
https://ask.sagemath.org/question/10429/finding-root-with-difficulties/?comment=17170#post-id-17170Hmmm. *Why* is my way of passing the function not appropriate? As far as I understand things, the lambda is a python construct, which sage will treat as a black box. The expression, on the other hand, is something sage can work with, e.g. by turning it into some optimized fastcall thingy. So why would a black box be more appropriate than a symbolic expression?Fri, 09 Aug 2013 15:00:32 +0200https://ask.sagemath.org/question/10429/finding-root-with-difficulties/?comment=17170#post-id-17170