#8321: numerical integration with arbitrary precision
-------------------------+--------------------------------------------------
Reporter: burcin | Owner:
Type: defect | Status: needs_work
Priority: major | Milestone: sage-4.7.2
Component: symbolics | Keywords: numerics,integration, sd32
Work_issues: | Upstream: N/A
Reviewer: | Author: Stefan Reiterer
Merged: | Dependencies:
-------------------------+--------------------------------------------------
Comment(by fredrik.johansson):
OK, here's a version also avoiding type conversions:
{{{
from mpmath.calculus.quadrature import TanhSinh
from mpmath import quad, mp
from sage.libs.mpmath.all import call, sage_to_mpmath, mpmath_to_sage
class TanhSinhRel(TanhSinh):
instance = None
def __new__(cls, ctx):
if not cls.instance:
cls.instance = object.__new__(cls)
cls.instance.ctx = ctx
cls.instance.standard_cache = {}
cls.instance.transformed_cache = {}
cls.instance.interval_count = {}
return cls.instance
def estimate_error(self, results, prec, epsilon):
mag = abs(results[-1])
if len(results) == 2:
return abs(results[0]-results[1]) / mag
try:
if results[-1] == results[-2] == results[-3]:
return self.ctx.zero
D1 = self.ctx.log(abs(results[-1]-results[-2])/mag, 10)
D2 = self.ctx.log(abs(results[-1]-results[-3])/mag, 10)
if D2 == 0:
D2 = self.ctx.inf
else:
D2 = self.ctx.one / D2
except ValueError:
print factorial(1000000)
return epsilon
D3 = -prec
D4 = min(0, max(D1**2 * D2, 2*D1, D3))
return self.ctx.mpf(10) ** int(D4)
class TanhSinhRel2(TanhSinhRel):
instance = None
def __init__(self, ctx):
pass
class TanhSinhRel3(TanhSinhRel2):
instance = None
def transform_nodes(self, nodes, a, b, verbose=False):
nodes = TanhSinh.transform_nodes(self, nodes, a, b, verbose)
prec = self.ctx.prec
nodes = [(mpmath_to_sage(x, prec), mpmath_to_sage(w, prec)) for
(x, w) in nodes]
return nodes
def sum_next(self, f, nodes, degree, prec, previous, verbose=False):
h = self.ctx.mpf(2)**(-degree)
if previous:
S = previous[-1]/(h*2)
else:
S = self.ctx.zero
s = 0
for (x, w) in nodes:
s += w * f(x)
return h * (S + sage_to_mpmath(s, prec))
def f(x):
return exp(-x**2) * log(x)
def g(x):
x = mpmath_to_sage(x, mp.prec)
y = exp(-x**2) * log(x)
return sage_to_mpmath(y, mp.prec)
mp.prec = 100
timeit("mp.quad(g, [17, 42], method=TanhSinhRel)")
timeit("mp.quad(g, [17, 42], method=TanhSinhRel2)")
timeit("mp.quad(f, [17, 42], method=TanhSinhRel3)")
print mp.quad(g, [17, 42], method=TanhSinhRel)
print mp.quad(g, [17, 42], method=TanhSinhRel2)
print mp.quad(f, [17, 42], method=TanhSinhRel3)
}}}
This prints:
{{{
5 loops, best of 3: 82.1 ms per loop
5 loops, best of 3: 72.5 ms per loop
5 loops, best of 3: 38.2 ms per loop
2.5657285005610514829173563961e-127
2.5657285005610514829173563961e-127
2.5657285005610514829173563961e-127
}}}
The estimate_error method still needs some work, though...
--
Ticket URL: <http://trac.sagemath.org/sage_trac/ticket/8321#comment:45>
Sage <http://www.sagemath.org>
Sage: Creating a Viable Open Source Alternative to Magma, Maple, Mathematica,
and MATLAB
--
You received this message because you are subscribed to the Google Groups
"sage-trac" group.
To post to this group, send email to [email protected].
To unsubscribe from this group, send email to
[email protected].
For more options, visit this group at
http://groups.google.com/group/sage-trac?hl=en.
