Finding root with difficulties

asked 2013-08-09 03:46:22 +0200

MvG
457 ●15 ●29 ●42

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.

answered 2013-08-09 05:16:08 +0200

Eviatar Bach
748 ●15 ●37 ●52

updated 2013-08-09 19:15:26 +0200

Unfortunately, the current root-finding only works with real numbers (see #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!

edit flag offensive delete link

Comments

Hmmm. *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?

MvG ( 2013-08-09 15:00:32 +0200 )edit

Okay, 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) nan

Eviatar Bach ( 2013-08-09 19:17:05 +0200 )edit

I 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.

Eviatar Bach ( 2013-08-09 19:21:48 +0200 )edit

add a comment

Your Answer

Please start posting anonymously - your entry will be published after you log in or create a new account.

Finding root with difficulties

1 Answer

Comments

Your Answer

Question Tools

Stats

Related questions

Finding root with difficulties edit

1 Answer

Comments

Your Answer

Question Tools

Stats

Related questions

Finding root with difficulties