Ask Your Question
1

How do I get solve () to use floats?

asked 2017-09-09 23:57:17 -0500

Lesparagus gravatar image

I want to solve a system of two (nonlinear) equations in floats. What the Sage documentation suggests I do is this:

var('x y')

eq1 = (x-x0)^2 +(y-y0)^2 == d0^2

eq2 = (x-x1)^2 + (y-y1)^2 == d1^2

solns = solve([eq1, eq2], x, y)

[[s[x].n(30), s[y].n(30)] for s in solns]

x0, x1, y0, y1, d0, and d1 are constants, by the way.

But this just solves it exactly and then approximates it. I don't want that! I need to solve this equation hundreds to thousands of times quickly, and the exact solving takes a really long time, creating a short delay to even solve one equation. I want to just solve in floats from the beginning.

I could solve the equations symbolically myself and plug them into the program as expressions in x and y. But the expressions are really long and not only is that a huge error-prone pain to type up (the expressions would take up several lines for each variable, not to mention that there are two solutions), but I have to solve several different equations that will be like this. I want to avoid doing this if at all possible.

I've tried using find_root, but I don't know how it works on systems, I've tried just putting floats into the equation, I've tried setting domain=float, RR, and RDF in the solve function, and a bunch of other random things that I don't even remember. How do I do this?

edit retag flag offensive close merge delete

1 answer

Sort by » oldest newest most voted
3

answered 2017-09-10 12:13:57 -0500

mforets gravatar image

updated 2017-09-10 12:38:13 -0500

In this case you could try to substitute the values in the exact solution. An example of this sort is given in pages 46-50 of Sage for power users (1).

Here is a minimal working example:

var('x, y, x0, y0, x1, y1, d0, d1')

# two circles
eq1 = (x-x0)^2 +(y-y0)^2 == d0^2
eq2 = (x-x1)^2 + (y-y1)^2 == d1^2

Let's ask Sage to solve exactly this system of equations:

%time solns = solve([eq1, eq2], x, y)
CPU times: user 7.28 s, sys: 8.01 ms, total: 7.28 s
Wall time: 7.28 s

It takes a couple of seconds, but it does a good job in finding the exact formulas,

sage: solns[0][0].rhs().variables()
(d0, d1, x0, x1, y0, y1)

Of course, the advantage is that the solutions obtained with this method are exact,

sage: solns[0][0].rhs().subs(x0=1, y0=2, x1=3, y1=4, d0=2, d1=1)
-1/8*sqrt(7) + 19/8
sage: _.n()
2.04428108611693

Let's see how much it takes to evaluate it numerically:

%timeit _ = solns[0][0].rhs().subs(x0=1, y0=2, x1=3, y1=4, d0=2, d1=1).n()
1000 loops, best of 3: 704 µs per loop

Now let's wrap the solution into a Cython function (1), (2):

%%cython
from math import sqrt
def sol1_x(double x0, double y0, double x1, double y1, double d0, double d1):
    return -0.5*(d0**2*x0 - d1**2*x0 - x0**3 + x0*x1**2 - x1**3 - x0*y0**2 - (x0 + x1)*y1**2 - (d0**2 - d1**2 - x0**2 + y0**2)*x1 + (2*x0*y0 + 2*x1*y0 + sqrt(-d0**4 - d1**4 + 2*d0**2*x0**2 - x0**4 + 4*x0*x1**3 - x1**4 - y0**4 + 4*y0*y1**3 - y1**4 + 2*(d0**2 + x0**2 + y0**2)*d1**2 + 2*(d0**2 + d1**2 - 3*x0**2 - y0**2)*x1**2 + 2*(d0**2 - x0**2)*y0**2 + 2*(d0**2 + d1**2 - x0**2 + 2*x0*x1 - x1**2 - 3*y0**2)*y1**2 - 4*(d0**2*x0 + d1**2*x0 - x0**3 - x0*y0**2)*x1 - 4*(d1**2*y0 + 2*x0*x1*y0 - x1**2*y0 - y0**3 + (d0**2 - x0**2)*y0)*y1))*y1 - sqrt(-d0**4 - d1**4 + 2*d0**2*x0**2 - x0**4 + 4*x0*x1**3 - x1**4 - y0**4 + 4*y0*y1**3 - y1**4 + 2*(d0**2 + x0**2 + y0**2)*d1**2 + 2*(d0**2 + d1**2 - 3*x0**2 - y0**2)*x1**2 + 2*(d0**2 - x0**2)*y0**2 + 2*(d0**2 + d1**2 - x0**2 + 2*x0*x1 - x1**2 - 3*y0**2)*y1**2 - 4*(d0**2*x0 + d1**2*x0 - x0**3 - x0*y0**2)*x1 - 4*(d1**2*y0 + 2*x0*x1*y0 - x1**2*y0 - y0**3 + (d0**2 - x0**2)*y0)*y1)*y0)/(x0**2 - 2*x0*x1 + x1**2 + y0**2 - 2*y0*y1 + y1**2)

Then,

%timeit _ = sol1_x(1.,2.,3.,4.,2.,1.)
100000 loops, best of 3: 7.03 µs per loop
sol1_x(1.,2.,3.,4.,2.,1.)
2.0442810861169263

That's a factor 100x of faster!

(1) if you are in the Jupyter notebook, create a new Code cell and in the first line use %%cython (and not %cython).

(2) Two observations:

  • replace the operator ^ with **, for example with the command str(solns[0][0].rhs()).replace('^', '**')
  • use 0.5 instead of 1/2

Finally, Sage does provide a function to do this automatically, (eg. operating with Sage objects defined in the global scope), it is cython_lambda. I wasn't able to apply it here.. strange.


If you have to evaluate thousands of values, then it will be even more efficient to work with arrays, see the Cython documentation for more info or ask here :)


Without passing through the exact solution obtained with solve, here is the same answer with SciPy's fsolve interface:

from scipy.optimize import fsolve

def circles(variables, *args):
    x, y = variables
    x0, y0, x1, y1, d0, d1 = args
    eq1 = (x-x0)^2 +(y-y0)^2 - d0^2
    eq2 = (x-x1)^2 + (y-y1)^2 - d1^2
    return (eq1, eq2)


%timeit _ = fsolve(circles, (-2.0, 2.0), args=(1.,2.,3.,4.,2.,1.))
1000 loops, best of 3: 1.24 ms per loop

sage: fsolve(circles, (-2.0, 2.0), args=(1.,2.,3.,4.,2.,1.))
array([ 2.04428109,  3.70571891])
edit flag offensive delete link more

Comments

Works like a charm! Thank you!

Also, how do you format code like that?

Lesparagus gravatar imageLesparagus ( 2017-09-10 18:04:53 -0500 )edit

@Lesparagus Indenting the block of code by 4 spaces should be enough to get syntax highlight. No need to surround the code by anything else. That should do the trick.

fidbc gravatar imagefidbc ( 2017-09-10 20:50:45 -0500 )edit

Alright, thank you!

Lesparagus gravatar imageLesparagus ( 2017-09-10 20:55:28 -0500 )edit

Your Answer

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

Add Answer

Question Tools

1 follower

Stats

Asked: 2017-09-09 23:49:32 -0500

Seen: 48 times

Last updated: Sep 10