Ask Your Question
1

find_root of function involving maximal eigenvalue

asked 2020-11-28 05:06:26 +0100

ericrozon gravatar image

I have defined a function returning the leading eigenvalue of a matrix L as follows (the first function is to handle what I think are numerical errors, where what should be a real eigenvalue is seen as x + 0.?e-80*I, for instance):

def is_essentially_real(x):
    if x.imag() == 0:
        return(True)
    else:
        return(False)

def get_leading_eigenvalue(L):
    evals = L.eigenvalues()
    moduli = [e.n() for e in evals if is_essentially_real(e)]
    moduli = [e for e in moduli if e >= 0]
    r = max(moduli)
    return(r)

I want to solve equations involving this function using find_root. As a toy example, I have the following:

def mat(b):
    M = matrix(RR, 2, 2, [1, b, 1, 2])
    return(M)

def f(b):
    M = mat(b)
    r = get_leading_eigenvalue(M)
    return(r)

b = var(b)
find_root(f(b) - 3, 1,3 ) #produces error

I know that f(2) = 3, so if I understand correctly, my last line of code should return the number 2. Instead, I get the following message:

TypeError: Cannot evaluate symbolic expression to a numeric value.

My best guess is that the method matrix.eigenvalues() does not play well with the method find_roots(), but I am not sure how to get around this. I have tried changing the matrix field to both of QQ and SR, to no avail. Any feedback is appreciated. Also, this is my first question on the site, so if I've committed any sins, please let me know.

edit retag flag offensive close merge delete

1 Answer

Sort by ยป oldest newest most voted
1

answered 2020-11-28 21:43:42 +0100

dsejas gravatar image

Hello, @ericrozon! Before answering your question, I would like to make you note that return is a sentence, not a function. That is, you should write something like return x instead of return(x). That being said, I can't think of any potential problem that could arise from using return(x), except that it can cause a momentarily raised eyebrow on somebody reading your code (like it happened to me.) A second observation is that b = var(b) is not syntactically correct, since the correct form actually is b = var('b'). That WILL cause problems. Finally, One last suggestion would be to avoid using the ring RR to define your matrix, since that can cause rounding errors.

Anyway, concerning your question, the problem seems to be the following: Sage is interpreting your variable b not as the independent variable of your function f, but as a symbolic constant. Basically, the problem comes from your function mat: feeding it b tries to create a symbolic matrix, which is impossible, since you declared its elements to be on the numerical ring RR. If you delete the RR as I suggested above, mat will try to create a symbolic matrix (that is allowed in that case, because Sage defaults to the symbolic ring SR), which will have symbolic eigenvalues, and you cannot compute the max of symbolic values

You can solve your problem simply doing a few modifications. Define the function

def g(x):
    return f(x) - 3

and call find_root on it. Alternatively, you could also define your f to return r - 3. I am going to use this second approach for this code:

def is_essentially_real(x):
    if x.imag() == 0:
        return True # changed
    else:
        return False # changed

def get_leading_eigenvalue(L):
    evals = L.eigenvalues()
    moduli = [e.n() for e in evals if is_essentially_real(e)]
    moduli = [e for e in moduli if e >= 0]
    r = max(moduli)
    return r # changed

def mat(b):
    M = matrix(2, 2, [1, b, 1, 2]) # changed: removed the RR ring
    return M # changed

def f(b):
    M = mat(b)
    r = get_leading_eigenvalue(M)
    return r - 3 # changed: subtracted 3

find_root(f, 1, 3) # no more error

This returns the root 1.9999999999999996 (pretty close!). You will notice that I marked every change I made to your code.

I hope this helps!

edit flag offensive delete link more

Comments

This is an incredibly detailed and generous answer - thanks a lot!! It works perfectly. One last question for future asks: how do I highlight numbers and sage code within a question the way you do? (Eg f looks nice in yours, but not in mine)

ericrozon gravatar imageericrozon ( 2020-11-29 02:40:38 +0100 )edit

To display blocks of code or error messages in Ask Sage, skip a line above and below, and do one of the following (all give the same result):

  • indent all code lines with 4 spaces
  • select all code lines and click the "code" button (the icon with '101 010')
  • select all code lines and hit ctrl-K

To display inline code fragments, use backticks.

For instance, typing

If we define `f` by

    def f(x, y, z):
        return x * y * z

then `f(2, 3, 5)` returns `30` but `f(2*3*5)` gives:

    TypeError: f() takes exactly 3 arguments (1 given)

produces:

If we define f by

def f(x, y, z):
    return x * y * z

then f(2, 3, 5) returns 30 but f(2*3*5) gives:

TypeError: f() takes exactly 3 arguments (1 given)
slelievre gravatar imageslelievre ( 2020-11-29 05:06:55 +0100 )edit

@ericrozon, you are welcome! I am very happy to help. Welcome to Ask SageMath!

dsejas gravatar imagedsejas ( 2020-11-29 06:07:31 +0100 )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: 2020-11-28 05:06:26 +0100

Seen: 217 times

Last updated: Nov 28 '20