ASKSAGE: Sage Q&A Forum - RSS feedhttps://ask.sagemath.org/questions/Q&A Forum for SageenCopyright Sage, 2010. Some rights reserved under creative commons license.Sat, 23 Jun 2018 11:27:09 +0200Semantics of RealBallField and RealIntervalField ?https://ask.sagemath.org/question/42722/semantics-of-realballfield-and-realintervalfield/I do not understand the meaning (or the output) of RealBallField elements. An example will be clearer :
sage: foo=RealBallField(10)(pi/2);foo
[1.57 +/- 2.27e-3]
So far, so good.
sage: sin(foo)
[1.0 +/- 3.92e-3]
An interval *centered* on 1, Isn't there a problem ? Mathematically, he upper bound is 1, and, since the derivative is 0, the maximal error is intuitively expected to be much smaller than the error on the argument. Let's find out analytically. If I'm not mistaken, the maximal error is :
sage: 1-sin(pi/2-2^-10).n()
4.76837120344520e-7
So I would expect something along the lines of [1-4.76*10^-7 1].
Maybe my expectation to use directly the printed notation is misguided. Let's use the specialized methods :
sage: sin(foo).below_abs()
[0.996 +/- 9.38e-5]
sage: sin(foo).above_abs()
[1.00 +/- 3.91e-3]
Three surprises :
* I'd expect the bounds to be ("sharp") reals, not balls...
* The stated radii make no sense to me/
* The stated lower bound is way too low.
Therefore, I do not understand the results. But I may be mistaken : my last formal contact with numerical analysis of errors is about 35 years old...
Can some kind soul point me to an explanation ?
Note : the RealIntervalField gives me *different* nonsense :
sage: bar=RealIntervalField(10)(pi/2);bar
1.58?
sage: bar.str(style="brackets")
'[1.5703 .. 1.5723]'
Okay so far...
sage: sin(bar)
1.000?
Possibly okay : value at the center of the argument. But not the center of the interval...
sage: sin(bar).str(style="brackets")
'[0.99902 .. 1.0000]'
The upper bound is correct ; I beleeve that the lower bound is too low...
sage: sin(bar).center()
1.0
Again, this is not the center of the interval, but the value at the center of the interval of the *argument*...
sage: sin(bar).absolute_diameter()
0.00098
I think this is too low...
**EDIT :** My understanding problems mat be better summarized by the following :
sage: foo=RealBallField(10)(pi/2);foo
[1.57 +/- 2.27e-3]
sage: bar=RealIntervalField(10)(pi/2);bar
1.58?
sage: foo.sin().upper()
1.1
sage: bar.sin().upper()
1.0Emmanuel CharpentierSat, 23 Jun 2018 11:27:09 +0200https://ask.sagemath.org/question/42722/find roots with arbitrary precisionhttps://ask.sagemath.org/question/9486/find-roots-with-arbitrary-precision/It seems that `find_root` will always convert its parameter range to `float`, even if `a` and `b` were originally given in some arbitrary precision real number field. Is there a variant of this algorithm that can find roots to arbitrary precision?
Even greater would be some algorithm that can make use of interval arithmetic based on `RealIntervalField`, which in particular will return an interval that is guaranteed to contain the zero. I have written some code along these lines:
def bisect_interval(f, i, d):
# find zero of function f in interval i with desired diameter d
# f must be monotonically increasing and must contain a zero in i
d2 = d/2
zero = i.parent()(0)
hi = i
while hi.absolute_diameter() > d2:
l, r = hi.bisection()
c = l.intersection(r)
fc = f(c)
if fc > zero:
hi = l
else:
hi = r
lo = i
while lo.absolute_diameter() > d2:
l, r = lo.bisection()
c = l.intersection(r)
fc = f(c)
if fc < zero:
lo = r
else:
lo = l
return lo.union(hi)
Am I reproducing functionality that's already available somewhere in Sage? If not, do you consider this functionality worth adding? Should it use some better algorithm than simple bisection? Is the algorithm even correct? You don't *have* to answer all of these questions, as my core question remains the one about an arbitrary precision version of `find_root`. But additional answers would be welcome.MvGWed, 31 Oct 2012 05:44:29 +0100https://ask.sagemath.org/question/9486/interface to MPFIhttps://ask.sagemath.org/question/9431/interface-to-mpfi/Thanks to the Sage team for providing an interface to the MPFI package for arbitrary precision interval arithmetic. It looks like the MPFI version used in Sage is 1.5.1, the most recent. Reading the included info file mpfi.info, there are references to additional functionality such as extended interval division, which would be useful for implementing extended interval Newton and related numerical methods. I can't tell if this is available in MPFI 1.5.1 and does not yet have a wrapper in Sage (not complaining!), or will only be "available soon" in MPFI as the file states.
Since I can't find any information at any MPFI site (Mon français n'est pas très bon), does anyone know if or when the extended interval arithmetic functions will be implemented in MPFI? And this looks like it would be quite a job to get into Sage as well. For example, currently:
1./RIF(-1.,1.)
[-infinity .. +infinity]
The extended interval division of this would be the (sharper bound) union of two intervals:
[-infinity .. -1] <union> [1 .. +infinity]
which would seem to require a data structure to hold multiple RIF objects. A numerical algorithm can then detect this split, perform iterations on both parts, and (ideally) find all solutions to the original problem without multiple starting points. This could get messier in higher dimensions since multiple intervals now become multiple boxes, cubes, and hypercubes.
I see that MPFI or something does solve systems of equations over RIF:
A = matrix([[RIF(2,2.1), RIF(1,1.1)],[RIF(1,1.1),RIF(2,2.0)]]); A
[2.1? 1.1?]
[1.1? 2]
b = vector([RIF(1,1.1),RIF(1,1.1)]); b
(1.1?, 1.1?)
x = A\b
x[0].endpoints()
(0.230244068953746, 0.426562500000001)
x[1].endpoints()
(0.259218749999999, 0.447175285884964)
which is very nice, but also breaks down if the determinant of the matrix is an interval that includes 0, resulting in all infinities for the solution. There is an entire area of numerical linear algebra now devoted to solving these problems, which I'm interested in implementing. I might be able to handle the extended interval arithmetic with my own generalized Python classes, but it's better to use whatever standards are being developed.
burningbrightMon, 15 Oct 2012 20:15:10 +0200https://ask.sagemath.org/question/9431/