Error from integration using mpmath quad

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!

Comments

Welcome to Ask Sage! Thank you for your question.

slelievre ( 3 years ago )

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 ( 3 years ago )

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 ( 3 years ago )

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

Emmanuel Charpentier ( 3 years ago )

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

samosonto ( 3 years ago )

add a comment

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)

# 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)

Comments

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

samosonto ( 3 years ago )

add a comment

Your Answer

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

Error from integration using mpmath quad

Comments

1 Answer

Comments

Your Answer

Question Tools

Stats

Related questions

Error from integration using mpmath quad savecancel

Comments

1 Answer

Comments

Your Answer

Question Tools

Stats

Related questions

Error from integration using mpmath quad