Error from integration using mpmath quad

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

samosonto
13 ●3

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

Comments

Welcome to Ask Sage! Thank you for your question.

slelievre ( 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 ( 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 ( 2022-03-03 22:33:39 +0200 )edit

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

Emmanuel 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 ( 2022-03-05 18:26:54 +0200 )edit

add a comment

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

Emmanuel Charpentier
7057 ●6 ●46 ●137

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

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)
    rp = mp.quad(lambda kx:mp.re(itg(kx)), kx_interval)
    ip = mp.quad(lambda kx:mp.im(itg(kx)), kx_interval)
    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
    return complex(mp.quad(Integrand(x, y), kx_interval))

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

HTH,

edit flag offensive delete link

Comments

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

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

add a comment

Error from integration using mpmath quad

Comments

1 Answer

Comments

Your Answer

Question Tools

Stats

Related questions

Error from integration using mpmath quad edit

Comments

1 Answer

Comments

Your Answer

Question Tools

Stats

Related questions

Error from integration using mpmath quad