# Error from integration using mpmath quad

Hello all,

I am trying to solve an integral with one variable kx and then I need the result from that integral to be evaluated in a point defined by x,y coordinates.

The code I am using is the following:

import numpy as np
from numpy.lib.scimath import sqrt
from mpmath import *
import matplotlib.pyplot as plt
mp.dps = 15

f = 8500
rho = 1.225
c0 = 343
Omega = 2*np.pi*f
k = Omega/c0
Z = -426

integrandReal = lambda kx, x, y: np.real(((2*sqrt(k**2 - kx**2)*Z)/(sqrt(k**2 - kx**2)*Z + Omega*rho))*((np.exp(1j*(kx*x + sqrt(k**2 - kx**2)*y)))/(sqrt(k**2 - kx**2))))
integrandImag = lambda kx, x, y: np.imag(((2*sqrt(k**2 - kx**2)*Z)/(sqrt(k**2 - kx**2)*Z + Omega*rho))*((np.exp(1j*(kx*x + sqrt(k**2 - kx**2)*y)))/(sqrt(k**2 - kx**2))))

integral = lambda x, y: quad(integrandReal, [-100*k, 100*k], [x], [y])  + 1j*quad(integrandImag, [-100*k, 100*k], [x], [y])
G = integral(0,0)


My problem is that for any value of x and ythat I use I get:

mpc(real='0.0', imag='0.0')


I appreciate any help to solve this!

edit retag close merge delete

( 2022-03-03 18:02:56 +0200 )edit
1

The problem may reside in the third and fourth arguments of quad(integrandImag, [-100*k, 100*k], [x], [y]). As far as I understand it, mpmath.quad expects a function as its first argument and integration intervals (possibly subdivided) as arguments after the first ; the intervals are lists of values. The documentation defines what is done for lists of length superior to 2, but not for length one intervals.

It is possible that mpmath tries ti integrate on the interval [x x], which is of null measure ; since your function is finite, this integral is 0.

"L'enfer est pavé de bonnes intentions." ("Hell is paved with good intents")...

You may try to pass a call to a two-argument functions x, y returning a one-argument xk function as your first argument...

( 2022-03-03 22:02:24 +0200 )edit

It makes sense what you say. I was trying to apply something similar to how scipy works. The problem of scipy is that the integration methods it uses don't let me achieve good numerical precision, but mpmath evidently works differently. I will try what you suggest. Thanks!

( 2022-03-03 22:33:39 +0200 )edit

Mixing numpy, scimath and mpmath seems to be unsupported. Why numpy ?

( 2022-03-05 16:19:11 +0200 )edit

My mistake! I was just playing before with other expressions and I left numpy in the code.

( 2022-03-05 18:26:54 +0200 )edit

Sort by » oldest newest most voted

Staying with mpmath... After evaluating :

reset()
import mpmath as mp

mp.dps = 15

f = 8500
rho = 1.225
c0 = 343
Omega = 2*mp.pi*f
k = Omega/c0
Z = -426

def Integrand(x, y):
# return a one-parameter anonymous function of kx at point x, y.
# This function returns an mpc object
return lambda kx: ((2*mp.sqrt(k**2 - kx**2)*Z)/(mp.sqrt(k**2 - kx**2)*Z + Omega*rho))*((mp.exp(1j*(kx*x + mp.sqrt(k**2 - kx**2)*y)))/(mp.sqrt(k**2 - kx**2)))

def integration(kx_interval, x, y):
# builds the one-parameter integrand and integrates them.
# Return the result as a Python complex
itg = Integrand(x, y)
return complex(rp + 1j*ip)


one gets :

sage: integration([-100*k, 100*k], 0, 0)
(356.37828188993564-16.743674967114078j)


Is that what you aimed at ?

Note that you do not need to separately integrate the real and imaginary parts :

# Shorter and cleaner
def integration1(kx_interval, x, y):
# builds the one-parameter integrand and integrate it.
# Return the result as a Python complex

sage: integration1([-100*k, 100*k], 0, 0)
(356.37828188993564-16.743674967114078j)


HTH,

more

This is exactly what I aimed at. Thank you very much!

( 2022-03-05 18:28:15 +0200 )edit