Ask Your Question

How to solve numericaly with arbitrary precision

asked 2020-03-09 10:19:57 +0200

JonasA gravatar image

updated 2020-03-09 17:19:42 +0200


I want to solve a numerical equation for which I can only access by a lambda function, the M matrix is to big for the determinant to be computed on the symbolic ring.

f = lambda om: RRR(M.subs({omega:om}).change_ring(RRR).det())

But , find-root solve in the built in float type of Python which is lacking the precision I need. Is there a way in sage to solve numericaly with arbitrary precision?

Thank you


(edited with example more in line with what I am trying to achieve)

sage: x=var('x')
sage: M= Matrix(SR,[[cos(x),cosh(x)],[sin(x),sinh(x)]])
sage: RRR = RealField(200)
sage: f = lambda om: M.subs({x:om}).change_ring(RRR).det()
sage: find_root(f,1,2)
edit retag flag offensive close merge delete

1 Answer

Sort by ยป oldest newest most voted

answered 2020-03-09 14:47:59 +0200

vdelecroix gravatar image

updated 2020-03-11 10:21:05 +0200

This depends very much on the type of equation you are considering. For polynomial equations

sage: x = polygen(ZZ)
sage: (x**5 - x - 1).roots(RealField(128))

If you want a more specific answer, ask a more precise question.


Sadly, Sage has no rountines for more general arbitrary general root finding. Here is a straightforward implementation of the dichotomy

def dichotomy(f, left, right, n):
    sl = f(left).sign()
    sr = f(right).sign()
    if sl == sr:
        raise ValueError
    for _ in range(n):
        middle = (left + right) / 2
        sm = f(middle).sign()
        if sm == 0:
            return sm
        elif sm == sl:
            left = middle
            right = middle
    return (left,right)

which you can use as

sage: f = lambda x: x.sin()
sage: R = RealField(256)
sage: dichotomy(f, R(3), R(3.5), 10)
sage: dichotomy(f, R(3), R(3.5), 256)

EDITED BIS: For a non-trivial example involving functions with cos/sin/cosh/sinh

sage: f(x) = -2*cos(x)**2 + 5 * sin(x) + cos(3*x) * cosh(x) - sinh(2*x - 1)
sage: g = fast_callable(f, vars=[x], domain=R)
sage: dichotomy(g, R(-0.3), R(0), 100)

To find the appropriate initial conditions for f you need to either carry on some analysis on paper or look at its graph via

sage: plot(f)
edit flag offensive delete link more


Thank you for answering. I have a matrix of hyperbolic and trigonometric function of a variable, I want to know for which values of the variable the matrix is singular. It is analog to the physics problem of knowing the natural frequencies of an oscillator. Although my matrix is not very dense, it is dense enough so that the determinant can't be computed in the symbolic ring. I have to use a lambda function to substitute de variable by its value, change the matrix to the real ring then compute the determinant (the computation of the determinant is too slow if I stay in the symbolic ring). Which is the lambda function I write on my first message with RRR a custom Real Field with 200 bits precision.

A code example will follow in the next message

JonasA gravatar imageJonasA ( 2020-03-09 17:11:11 +0200 )edit


sage: x=var('x')
sage: M= Matrix(SR,[[cos(x),cosh(x)],[sin(x),sinh(x)]])
sage: RRR = RealField(200)


 sage: f = lambda om: RRR(M.subs({x:om}).change_ring(RRR).det())

gives f a very long runtime, find_root don't converge on a reasonable time and

sage: f = lambda om: M.subs({x:om}).change_ring(RRR).det()
sage: find_root(f,1,2)

while f has a much faster runtime, throw an exception

JonasA gravatar imageJonasA ( 2020-03-09 17:13:21 +0200 )edit

Dichotomy solve the problem although my function is too much unsmooth to be used so I will stop trying. My function was taking those value: (f(37) was returning 0, but when I used a better resolution it was more 1e67)

sage: f(36.99)
sage: f(36.999)
sage: f(36.99999)
sage: f(36.9999999)
sage: f(36.999999999)
JonasA gravatar imageJonasA ( 2020-03-12 14:07:01 +0200 )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



Asked: 2020-03-09 10:19:57 +0200

Seen: 714 times

Last updated: Mar 11 '20