On Nov 27, 2007 9:25 AM, Ted Kosan <[EMAIL PROTECTED]> wrote:
>
> Ondrej wrote:
>
> > > Could you please clarify, what exact functionality in solve you expect
> > > in order for 1235 to be solved?
> > >
> > > Should it just run the iterative numerical solver if it cannot find
> > > the solution analytically?
>
>
> And William wrote:
>
> > I don't know. However, Ted, what do you think of the following, i.e.,
> > it is a way in Sage to solve your problem which is probably pretty
> > clean and flexible, and could certainly made a little more student
> > friendly?
> >
> > sage: var('t')
> > sage: a = .004*(8*e^(-(300*t)) - 8*e^(-(1200*t)))*(720000*e^(-(300*t))
> > - 11520000*e^(-(1200*t))) +.004*(9600*e^(-(1200*t)) -
> > 2400*e^(-(300*t)))^2
> > sage: from scipy.optimize import brentq
> > sage: # Given two points x, y such that a(x) and a(y) have different sign,
> > this
> > sage: # brentq uses "inverse quadratic extrapolation" to find a root of a
> > in the
> > sage: # interval [x,y]. It has lots of extra tolerance and other options.
> > sage: brentq(a, 0, 0.002)
> > 0.00041105140493493411
> > sage: show(plot(a,0,.002),xmin=0, xmax=.002)
> >
> > I.e., what we provide an numerical_root method so that
> > a.numerical_root(x,y)
> > would fine a numerical root of a in the interval [x,y], if it exists?
> > It could be built on brentq. The main thing we would have to add
> > is some sort of analysis to find x', y' in the interval so that a(x')
> > has different sign from a(y'), i.e., decide if there is a sign switch,
> > which could be doable for many analytically defined functions at least.
>
> Here is an excerpt from a Mathematica FAQ that I located on the Internet:
> -----
> 3.2 I've properly entered a Solve command but all Mathematica returns
> is an empty list!
> What's going on:
>
> You've asked Mathematica to solve an equation it can't solve
> analytically. So instead
> of a list of solutions, it gives you an empty list. The same thing can
> happen, incidentally, with NSolve.
>
> How to fix the problem: Try using FindRoot to solve the equation.
> First write the equation in the form
> expression = 0.
> (1)
> Then use Plot to graph the expression. Use Mathematica's coordinate
> locator to determine roughly where
> the zeros of the expression are. Feed these to FindRoot as initial guesses.
> -----
>
> It appears that Mathematica uses the same technique that you describe
> using brentq to solve this problem.
>
> Also, this recent discussion on sage-developer seems to be related to
> this issue:
>
> http://www.mail-archive.com/[EMAIL PROTECTED]/msg06571.html
>
>
> For the engineering oriented problems like the two I originally
> submitted, we are usually interested in numeric results. I am now
> thinking that having functions like nsolve() and find_root() in SAGE
> would serve our needs better than enhancing the solve() function.
>
> What is coming to mind is that nsolve() would work like Mathematica's
> NSolve function
> (http://documents.wolfram.com/mathematica/functions/NSolve) and
> find_root() would be a wrapper around brentq.
>
> Does this seem reasonable?
Ted,
As a start I've implemented find_root (and some minizing and
maximizing functions)
and posted a patch here:
http://trac.sagemath.org/sage_trac/ticket/1235
If you look at the patch you'll see examples like this:
sage: t = var('t')
sage: v = 0.004*(9600*e^(-(1200*t)) - 2400*e^(-(300*t)))
sage: v.find_root(0, 0.002)
0.0015403270679114178
sage: a = .004*(8*e^(-(300*t)) -
8*e^(-(1200*t)))*(720000*e^(-(300*t)) - 11520000*e^(-(1200*t)))
+.004*(9600*e^(-(1200*t)) - 2400*e^(-(300*t)))^2
There is a 0 very close to the origin:
sage: show(plot(a, 0, .002),xmin=0, xmax=.002)
Using solve does not work to find it:
sage: a.solve(t)
[]
However \code{find_root} works beautifully:
sage: a.find_root(0,0.002)
0.00041105140493493411
----
One interesting thing is that I made it so that
f.find_root(a,b)
works even if the sign of f(a) and f(b) are the same. In that case,
it will find a min or max of f on the interval, and use that as a new
endpoint as input to the root finding algorithm. The root finder
itself is some serious numerical code from scipy...
Also
find_root(f, a, b)
will now work for any callable python object -- that takes floats
as input and outputs floats.
William
--~--~---------~--~----~------------~-------~--~----~
To post to this group, send email to sage-support@googlegroups.com
To unsubscribe from this group, send email to [EMAIL PROTECTED]
For more options, visit this group at
http://groups.google.com/group/sage-support
URLs: http://sage.math.washington.edu/sage/ and http://sage.scipy.org/sage/
-~----------~----~----~----~------~----~------~--~---
