Opened 14 years ago

Closed 13 years ago

Last modified 13 years ago

#3870 closed defect (invalid)

scipy - bug in find_root

Reported by: was Owned by: gfurnish
Priority: major Milestone: sage-duplicate/invalid/wontfix
Component: calculus Keywords: scipy
Cc: jason Merged in:
Authors: Reviewers:
Report Upstream: N/A Work issues:
Branch: Commit:
Dependencies: Stopgaps:

GitHub link to the corresponding issue

Description

It appears that find_root give the wrong result in the following
example:

----------------------------------------------------------------------
| SAGE Version 3.0.5, Release Date: 2008-07-11                       |
| Type notebook() for the GUI, and license() for information.        |
----------------------------------------------------------------------

sage: var('av')
av
sage: f=2.00000000000000*av - (10*sqrt(36.0000000000000*(av^4 - 2*av^3
+ av^2) + 207.360000000000*(-4*av^4 + 8*av^3 - 4*av^2) +
365*(51.8400000000000*(4*av^3 - 4*av) + 9.00000000000000*(2*av -
2*av^2)) + 133225*(2.25000000000000 - 51.8400000000000*av)) +
60.0000000000000*(av - av^2) + 5475.00000000000)/(720.000000000000*(2
- 2*av) + 131400.000000000)
sage: find_root(f,0,0.1)
0.0
sage: plot(f,0,0.1)

The last command will plot f and reveal that the root is somewhere
near to 0.05, but certainly not 0.0. No matter what I use for xtol or
rtol, I always get 0.0 as a result. Can anyone help?

Thanks!

Stan

I don't what is going on yet. It may be a bug in scipy.

Change History (7)

comment:1 Changed 14 years ago by was

Stan Schymanski points out that:

I just found out that the problem can be solved using maxima's
find_root:

sage: maxima.find_root(f,0,0.1)
.03134446907530818

It's great that sage offers multiple ways of doing the same thing. :)

Stan

comment:2 Changed 14 years ago by anakha

I don't know if that is the cause, but I don't think find_root is prepared to deal with complex numbers:

sage: var('av')
av
sage: f=2.00000000000000*av - (10*sqrt(36.0000000000000*(av^4 - 2*av^3
+ av^2) + 207.360000000000*(-4*av^4 + 8*av^3 - 4*av^2) +
365*(51.8400000000000*(4*av^3 - 4*av) + 9.00000000000000*(2*av -
2*av^2)) + 133225*(2.25000000000000 - 51.8400000000000*av)) +
60.0000000000000*(av - av^2) + 5475.00000000000)/(720.000000000000*(2
- 2*av) + 131400.000000000)
sage: f(0.1).n()
0.158699584011575 - 0.0475301543661689*I

comment:3 Changed 14 years ago by anakha

Continuing from the previous block:

sage: ff = sage.ext.fast_eval.fast_float(f, av)
sage: ff
<sage.ext.fast_eval.FastDoubleFunc object at 0x8a7a8f0>
sage: ff(0)
-0.082429990966576328
sage: ff(0.1)
nan
sage: ff(0.05)
nan
sage: ff(0.04)
0.02790669038508465
sage: find_root(f, 0, 0.04)
0.031344469075307836

Yeah, it is clearly not prepared. I don't know why maxima works.

comment:4 Changed 14 years ago by jason

Cc: jason added
Keywords: scipy added

comment:5 Changed 14 years ago by jason

Summary: bug in find_rootscipy - bug in find_root

comment:6 Changed 13 years ago by jason

Report Upstream: N/A
Resolution: invalid
Status: newclosed

Yep. In the scipy documentation, it assumes the function is continuous (from context, it seems over the reals). So this function fails, as allowed by the documentation.

comment:7 Changed 13 years ago by mvngu

Milestone: sage-4.3.1sage-duplicate/invalid/wontfix
Note: See TracTickets for help on using tickets.