Ask Your Question

Error from integration using mpmath quad

asked 2022-03-03 17:37:54 +0200

samosonto gravatar image

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 flag offensive close merge delete


Welcome to Ask Sage! Thank you for your question.

slelievre gravatar imageslelievre ( 2022-03-03 18:02:56 +0200 )edit

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...

Emmanuel Charpentier gravatar imageEmmanuel Charpentier ( 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!

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

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

Emmanuel Charpentier gravatar imageEmmanuel Charpentier ( 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.

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

1 Answer

Sort by » oldest newest most voted

answered 2022-03-05 16:29:27 +0200

Emmanuel Charpentier gravatar image

updated 2022-03-05 18:22:43 +0200

Staying with mpmath... After evaluating :

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)
    rp = mp.quad(lambda, kx_interval)
    ip = mp.quad(lambda, kx_interval)
    return complex(rp + 1j*ip)

one gets :

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

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
    return complex(mp.quad(Integrand(x, y), kx_interval))

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


edit flag offensive delete link more


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

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

Your Answer

Please start posting anonymously - your entry will be published after you log in or create a new account.

Add Answer

Question Tools

1 follower


Asked: 2022-03-03 17:02:43 +0200

Seen: 50 times

Last updated: Mar 05