github-featured/mpmath--mpmath
2
0
Fork 0
Code Issues Pull Requests Deployments Models Agent Notes Packages Projects Releases Wiki Activity

naive quadrature using adaptive subdivision

Browse Source
This commit is contained in:
fredrik
2021-02-26 13:05:54 +01:00
parent c6a35f9ee7
commit 7123a67a1d
5 changed files with 117 additions and 5 deletions
Download Patch File Download Diff File Expand all files Collapse all files
doc/source/calculus/integration.txt
+5
View File
@@ -6,6 +6,11 @@ Standard quadrature (``quad``)
.. autofunction:: mpmath.quad
Quadrature with subdivision (``quadsubdiv``)
............................................
.. autofunction:: mpmath.quadsubdiv
Oscillatory quadrature (``quadosc``)
....................................
mpmath/__init__.py
+1
View File
@@ -103,6 +103,7 @@ quad = mp.quad
quadgl = mp.quadgl
quadts = mp.quadts
quadosc = mp.quadosc
quadsubdiv = mp.quadsubdiv
invertlaplace = mp.invertlaplace
invlaptalbot = mp.invlaptalbot
mpmath/calculus/quadrature.py
+104
View File
@@ -718,6 +718,8 @@ class QuadratureMethods(object):
>>> quad(f, [-100, 0, 100]) # Also good
3.12159332021646
For such integr
**References**
1. http://mathworld.wolfram.com/DoubleIntegral.html
@@ -1006,6 +1008,108 @@ class QuadratureMethods(object):
s += ctx.nsum(term, [n, ctx.inf])
return s
def quadsubdiv(ctx, f, interval, tol=None, maxintervals=None, **kwargs):
"""
Computes the integral of *f* over the interval or path specified
by *interval*, using :func:`~mpmath.quad` together with adaptive
subdivision of the interval.
This function gives an accurate answer for some integrals where
:func:`~mpmath.quad` fails::
>>> mp.dps = 15; mp.pretty = True
>>> quad(lambda x: abs(sin(x)), [0, 2*pi])
3.99900894176779
>>> quadsubdiv(lambda x: abs(sin(x)), [0, 2*pi])
4.0
>>> quadsubdiv(sin, [0, 1000])
0.437620923709297
>>> quadsubdiv(lambda x: 1/(1+x**2), [-100, 100])
3.12159332021646
>>> quadsubdiv(lambda x: ceil(x), [0, 100])
5050.0
>>> quadsubdiv(lambda x: sin(x+exp(x)), [0,8])
0.347400172657248
The argument *maxintervals* can be set to limit the permissible
subdivision::
>>> quadsubdiv(lambda x: sin(x**2), [0,100], maxintervals=5, error=True)
(-5.40487904307774, 5.011)
>>> quadsubdiv(lambda x: sin(x**2), [0,100], maxintervals=100, error=True)
(0.631417921866934, 1.10101120134116e-17)
Subdivision does not guarantee a correct answer since, the error
estimate on subintervals may be inaccurate::
>>> quadsubdiv(lambda x: sech(10*x-2)**2 + sech(100*x-40)**4 + sech(1000*x-600)**6, [0,1], error=True)
(0.209736068833883, 1.00011000000001e-18)
>>> mp.dps = 20
>>> quadsubdiv(lambda x: sech(10*x-2)**2 + sech(100*x-40)**4 + sech(1000*x-600)**6, [0,1], error=True)
(0.21080273550054927738, 2.200000001e-24)
The second answer is correct. We can get an accurate result at lower
precision by forcing a finer initial subdivision::
>>> quadsubdiv(lambda x: sech(10*x-2)**2 + sech(100*x-40)**4 + sech(1000*x-600)**6, linspace(0,1,5))
0.210802735500549
The following integral is too oscillatory for convergence, but we can get a
reasonable estimate::
>>> v, err = fp.quadsubdiv(lambda x: fp.sin(1/x), [0,1], error=True)
>>> round(v, 6), round(err, 6)
(0.504067, 1e-06)
>>> sin(1) - ci(1)
0.504067061906928
"""
queue = []
for i in range(len(interval)-1):
queue.append((interval[i], interval[i+1]))
total = ctx.zero
total_error = ctx.zero
if maxintervals is None:
maxintervals = 10 * ctx.prec
count = 0
quad_args = kwargs.copy()
quad_args["verbose"] = False
quad_args["error"] = True
if tol is None:
tol = +ctx.eps
orig = ctx.prec
try:
ctx.prec += 5
while queue:
a, b = queue.pop()
s, err = ctx.quad(f, [a, b], **quad_args)
if kwargs.get("verbose"):
print("subinterval", count, a, b, err)
if err < tol or count > maxintervals:
total += s
total_error += err
else:
count += 1
if count == maxintervals and kwargs.get("verbose"):
print("warning: number of intervals exceeded maxintervals")
if a == -ctx.inf and b == ctx.inf:
m = 0
elif a == -ctx.inf:
m = min(b-1, 2*b)
elif b == ctx.inf:
m = max(a+1, 2*a)
else:
m = a + (b - a) / 2
queue.append((a, m))
queue.append((m, b))
finally:
ctx.prec = orig
if kwargs.get("error"):
return +total, +total_error
else:
return +total
if __name__ == '__main__':
import doctest
doctest.testmod()
mpmath/functions/elliptic.py
+2 -2
View File
@@ -586,8 +586,8 @@ def RJ_calc(ctx, x, y, z, p, r, integration):
N += margin
F = lambda t: 1/(ctx.sqrt(t+x)*ctx.sqrt(t+y)*ctx.sqrt(t+z)*(t+p))
if integration == 2:
return 1.5 * ctx.quad(F, [0, N, ctx.inf])
initial_integral = 1.5 * ctx.quad(F, [0, N])
return 1.5 * ctx.quadsubdiv(F, [0, N, ctx.inf])
initial_integral = 1.5 * ctx.quadsubdiv(F, [0, N])
x += N; y += N; z += N; p += N
xm,ym,zm,pm = x,y,z,p
A0 = Am = (x + y + z + 2*p)/5
mpmath/tests/test_elliptic.py
+5 -3
View File
@@ -654,9 +654,11 @@ def test_elliptic_integrals():
assert elliprj(0.3068, -4.037+0.0632j, 1.654, -0.9609).ae(-1.17157627949475577 - 0.069182614173988811j, abs_eps=1e-8)
assert elliprj(0.3068, -4.037+0.00632j, 1.654, -0.9609).ae(-1.17337595670549633 - 0.0623069224526925j, abs_eps=1e-8)
# these don't work because of inaccurate integration
# assert elliprj(0.3068, -4.037-0.0632j, 1.654, -0.9609).ae(1.77940452391261626 + 0.0388711305592447234j, abs_eps=1e-8)
# assert elliprj(0.3068, -4.037-0.00632j, 1.654, -0.9609).ae(1.77806722756403055 + 0.0592749824572262329j, abs_eps=1e-8)
# these require accurate integration
assert elliprj(0.3068, -4.037-0.0632j, 1.654, -0.9609).ae(1.77940452391261626 + 0.0388711305592447234j)
assert elliprj(0.3068, -4.037-0.00632j, 1.654, -0.9609).ae(1.77806722756403055 + 0.0592749824572262329j)
# issue #571
assert ellippi(2.1 + 0.94j, 2.3 + 0.98j, 2.5 + 0.01j).ae(-0.40652414240811963438 + 2.1547659461404749309j)
assert ellippi(2.0-1.0j, 2.0+1.0j).ae(1.8578723151271115 - 1.18642180609983531j)
assert ellippi(2.0-0.5j, 0.5+1.0j).ae(0.936761970766645807 - 1.61876787838890786j)