0

I am trying to calculate an integral with mpmath.quad. I basically have to calculate three moments of a distribution, call it f (pseudocode):

integrate(f(x)/x, 0, infinity)
integrate(f(x)*x, 0, infinity)
integrate(f(x)ln(x)/x, 0, infinity)

As far as I understood, the tanh-sinh algorithm from mpmath applies a coordinate transformation to the integration interval, uses it to find suitable nodes xk-s and weights wk-s and returns (pseudocode):

sum f(xk)*wk

with N the number of nodes. Since in my original problem the calculation takes a long time, I was wishing to reuse the values of f calculated by the first integration for the other integrals, i.e., computing them from discrete samples with something like scipy.integrate.trapezoid or simpson. By employing a structure Simulator I simplified from an answer on this forum, I managed to cache the xk-s and f(xk)-s, but not the weights. I checked that the xk-s from different integrations are the same.

Now, if I naively apply a standard quadrature from the scipy module to my samples, say trapezoid(fs, xs), the result I get is different from that calculated by mpmath.quad, especially if the samples are few. While this fact does not surprise me, I would like to find out how to retrieve the weights mpmath.quad uses in the first calculation, say the one for f(x)/x, so I could avoid running the time consuming algorithm thrice.

I cannot understand the documentation as for this point. mpmath documentation gives a lot of examples, but none concerning the retrieval of calculated nodes, although it states several times the nodes are "cached". Where are they? mpmath.quad only returns the integral result!

So I would like to know: is what I'm trying to achieve sensible at all? And if it is, how could I accomplish the task?

Below is a code that reproduces the behaviour. Any help is very appreciated.

import numpy as np
import mpmath as mp
import matplotlib.pyplot as plt
from scipy.integrate import trapezoid, simpson

class Simulator:
    def __init__(self, func, storex:np.ndarray=np.array([]), storef:np.ndarray=np.array([])):
        self.func = func
        self.storex = storex
        self.storef = storef
        
    def simulate(self, x, *args):
        result = self.func(x, *args)
        self.storex = np.append(self.storex, x)
        self.storef = np.append(self.storef, result)
        
        return result

def lorentz(x):
    x = mp.mpf(x)
    return mp.mpf(1)/(mp.power(x-mp.mpf(1), 2) + mp.mpf(1))


simratio = simproduct = Simulator(lorentz)

integralratio = mp.quad(lambda x: simratio.simulate(x)/x, [0, 1, mp.inf])
integralproduct = mp.quad(lambda x: simproduct.simulate(x)*x, [0, 1, mp.inf])

ratiox, ratioy = np.transpose(sorted([(x,y) for x, y in zip(simratio.storex, simratio.storef)])) #we get the sampled xs and f(x)s

integralproduct_trap = trapezoid(ratioy*ratiox, ratiox) 


print("int[f(x)/x], quad: ", integralratio)
print("int[f(x)*x], quad: ", integralproduct)
print("int[f(x)*x], scipy.trapezoid: ", integralproduct_trap)
EXB-over-mu0
  • 1

0 Answers0