1

I am trying to integrate an oscillatory expression with mpmath.quad. The integral only has one variable but I need to evaluate it in a 2D space. I need the resulting function from the integral along kx to be evaluated for a group of coordinated in x and y that defined a plane. That's why I evaluate integral(x,y) as a test in a single point of that 2D array.

This is the code I am using:

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], maxdegree=10)[0]  + 1j*quad(integrandImag, [-100*k, 100*k], maxdegree=10)[0]
G = integral(0.1,0.1)

I am getting the following errors when I try to evaluate the integral in an arbitrary point:

Traceback (most recent call last):

  File "<ipython-input-87-901e6e5dc630>", line 1, in <module>
    integral(1,1)

  File "/d/dfg/Documents/ImpedanceModel/GreenImpedance.py", line 55, in <lambda>
    integral = lambda x, y: quad(integrandReal, [-100*k, 100*k], maxdegree=10)[0]  + 1j*quad(integrandImag, [-100*k, 100*k], maxdegree=10)[0]

  File "/opt/tools/anaconda/2020.11/lib/python3.8/site-packages/mpmath/calculus/quadrature.py", line 742, in quad
    v, err = rule.summation(f, points[0], prec, epsilon, m, verbose)

  File "/opt/tools/anaconda/2020.11/lib/python3.8/site-packages/mpmath/calculus/quadrature.py", line 232, in summation
    results.append(self.sum_next(f, nodes, degree, prec, results, verbose))

  File "/opt/tools/anaconda/2020.11/lib/python3.8/site-packages/mpmath/calculus/quadrature.py", line 304, in sum_next
    S += self.ctx.fdot((w,f(x)) for (x,w) in nodes)

  File "/opt/tools/anaconda/2020.11/lib/python3.8/site-packages/mpmath/ctx_mp_python.py", line 944, in fdot
    for a, b in A:

  File "/opt/tools/anaconda/2020.11/lib/python3.8/site-packages/mpmath/calculus/quadrature.py", line 304, in <genexpr>
    S += self.ctx.fdot((w,f(x)) for (x,w) in nodes)

TypeError: <lambda>() missing 2 required positional arguments: 'x' and 'y'

I thought the error was because I wasn't defining the arguments in the integration as it has to be done with scipy.integrate.quad (as you can see in one of my last posts) but with mpmath I don't know what is generating this error.

Thanks in advance for any help!

  • Your lambdas (e.g. `integrandReal`) require 3 input arguments and you are passing only one in `quad` (you can see quad's syntax [here](https://mpmath.org/doc/current/calculus/integration.html)) – nonDucor Mar 03 '22 at 13:03
  • Also, what do you mean with "The integral only has one variable but I need to evaluate it in a 2D space". If the integral has only one variable, it cannot be evaluated over a 2D space. – nonDucor Mar 03 '22 at 13:03
  • @nonDucor I will try to explain it better. I need the resulting integral to be evaluated for a group of coordinated in x and y that defined a plane. That's why I evaluate integral (x,y) as a test in a single point of that 2D array. – Monosgaitnas Mar 03 '22 at 13:06
  • That's why I define first my whole function in terms of (kx, x, y), then I integrate in kx and then I evaluate the resulting function in (x,y). It's not a multivariable integral. – Monosgaitnas Mar 03 '22 at 13:12
  • With `scipy quad` parameters like this are either global to the function, or passed in via the `args` tuple. If the mpmath doesn't provide `args` then you need to go the global route. – hpaulj Mar 03 '22 at 15:56
  • Yes! It was solved with the modification that @nonDucor made to the code below. Now the problem is in terms of numerical results. Thanks! – Monosgaitnas Mar 03 '22 at 16:35

1 Answers1

1

You need to pass x and y to your lambdas (integrandReal and integrandImag). Also, there is another issue with your example as you try to subscript the mpc object returned by quad. I removed it to get a result, it shows you how to pass x and y through the calls, but I cannot assure you the numerical result is correct:

integral = lambda x, y: quad(integrandReal, [-100*k, 100*k], [x], [y], maxdegree=10)  + 1j*quad(integrandImag, [-100*k, 100*k], [x], [y], maxdegree=10)
G = integral(0.1,0.1)
# Returns: mpc(real='0.0', imag='0.0')
nonDucor
  • 2,057
  • 12
  • 17
  • Thanks! There are no errors when I run the code, nevertheless for any point I evaluate in integral(x,y) I get mpc(real='0.0', imag='0.0') :/ – Monosgaitnas Mar 03 '22 at 15:31