Ask Your Question
1

Problem with find_root

asked 2012-10-08 08:28:08 +0100

mattias gravatar image

updated 2012-10-08 08:29:18 +0100

Hello, I'm running into a problem with "find_root." To give an example, when I use the following code, "find_root" returns a value that is not a root:


sage: var('q')

sage: h = ((q+.2(40 - q))^-0.5)/((40-q+.2(q))^-0.5) == ((q+.2(40 - q))^.5-(50)^.5+(50-8)^.5)/((40-q+.2(q))^0.5-(50)^.5+(50-32)^.5)

sage: h_a = find_root(h, 0, 40)

sage: print h_a

sage: #h_a = 11.76222222321

sage: t = ((h_a+.2(40 - h_a))^-0.5)/((40-h_a+.2(h_a))^-0.5)

sage: print t

sage: z = ((h_a+.2(40 - h_a))^.5-(50)^.5+(50-8)^.5)/((40-h_a+.2(h_a))^0.5-(50)^.5+(50-32)^.5)

sage: print z

40.0

0.447213595499989

3.56694309609917e13

As you can see, sage returns a root of 40, which is not a root (to check this, note that t does not equal z). My question is why find_root is doing this. There are several ways in which I could fix the problem "locally;" e.g. setting the bounds between 0 and 30 returns the correct root of 11.76, but I need to eliminate the error to apply my code to a more general problem. Please let me know if I am doing something wrong here!

Thanks!

edit retag flag offensive close merge delete

2 Answers

Sort by ยป oldest newest most voted
1

answered 2012-10-08 11:57:09 +0100

achrzesz gravatar image

Mpmath can give you the error message

f=lambda q:((q+.2*(40. - q))^-0.5)/((40.-q+.2*(q))^-0.5) -( ((q+.2*(40. - q))^.5-(50.)^.5+(50.-8.)^.5)/((40.-q+.2*(q))^0.5-(50.)^.5+(50.-32.)^.5))

from mpmath import *
mp.dps = 15; mp.pretty = True
s1=findroot(f,11)
print s1            # no error message
s1=findroot(f,40)   # error message
edit flag offensive delete link more
0

answered 2012-10-08 09:16:46 +0100

achrzesz gravatar image

updated 2012-10-08 09:20:32 +0100

You can't include 40 into interval in this case

plot(h,(q,0,45),ymin=-50,ymax=50)
edit flag offensive delete link more

Comments

Thank you achrzesz. right, is there any way to get find_root to return the root over the interval (0,40) instead of returning the `discontinuity'? To clarify, the problem is that it does return a value, but not the correct one (or an error message).

mattias gravatar imagemattias ( 2012-10-08 10:00:21 +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

Stats

Asked: 2012-10-08 08:28:08 +0100

Seen: 625 times

Last updated: Oct 08 '12