CODEWITHSUNDEEP
pythonfunctionscipynumerical-integrationquad
George
George

Reputation: 232

Defining a integral of multi variable function as a second function python

Defining a integral of multi variable function as a second function python

I am using python to integrate a multivariable function over only one of the variables (it is a function of x and theta and I am integrating over theta from 0 to 2*pi, hence the result is a function of x). I have attempted the following:

import numpy as np
import scipy.integrate as inte
d=10.0
xvals=np.linspace(-d,d,1000)
def aIntegrand(theta,x):
    return 1/(2*np.pi)*np.sin(2*np.pi*x*np.sin(theta)/d)**2
def A(x):
    return (inte.quad(aIntegrand,0,2*np.pi,args=(x,))[0])**(1/2)

plt.plot(xvals,A(xvals))
plt.xlabel("x")
plt.ylabel("A(x)")
plt.show()

and I get the following error:

TypeError: only size-1 arrays can be converted to Python scalars

I assume this is because the result of the quad integrator is an array with two elements and python doesn't like defining the function based on an indexed array? That is a complete guess of the problem though. If anyone knows how I can fix this and could let me know, that'd be great :)

A second attempt

I have managed to get the plot of the integral successfully using the following code:

import numpy as np
import scipy.integrate as inte
import matplotlib.pyplot as plt
d=10.0
xvals=np.linspace(-d,d,1000)
thetavals=np.linspace(0.0,2*np.pi,1000)
def aIntegrand(theta,x):
    return 1/(2*np.pi)*np.sin(2*np.pi*x*np.sin(theta)/d)**2
def A(x):
    result=np.zeros(len(x))
    for i in range(len(x)):
        result[i]=(inte.quad(aIntegrand,0,2*np.pi,args=(x[i],))[0])**(1/2)
    return result
def f(x,theta):
    return x**2* np.sin(theta)

plt.plot(xvals,A(xvals))
plt.xlabel("x")
plt.ylabel("A(x)")
plt.show()

Integral Plot

But this does not give A(x) as a function, due to the way I have defined it, it requires an array form input. I need the function to be of the same form as aIntegrand, where when given parameters returns a single value & so the function can be integrated repeatedly.

Upvotes: 1

Views: 1013

Answers (2)

Nico Schlömer
Nico Schlömer

Reputation: 58821

First note that your integral can be calculated analytically. It is

0.5 * (1 - J0(4 * pi * x / d))

where J0 is the Bessel function of the first kind.

Second, you could use quadpy (one of my projects); it has fully vectorized computation.

import numpy as np
import quadpy
import matplotlib.pyplot as plt
import scipy.special

d = 10.0
x = np.linspace(-d, d, 1000)


def aIntegrand(theta):
    return (
        1
        / (2 * np.pi)
        * np.sin(2 * np.pi * np.multiply.outer(x, np.sin(theta)) / d) ** 2
    )


Ax2, err = quadpy.quad(aIntegrand, 0, 2 * np.pi)
Ax = np.sqrt(Ax2)

plt.plot(x, Ax, label="quadpy")
plt.xlabel("x")
plt.ylabel("A(x)")
plt.plot(x, np.sqrt(0.5 * (1 - scipy.special.jv(0, 4*np.pi*x/d))), label="bessel")
# plt.show()
plt.savefig("out.png", transparent=True, bbox_inches="tight")

enter image description here

Upvotes: 1

Patol75
Patol75

Reputation: 4547

I do not think that what you seek exists within Scipy. However, you have at least two alternatives.

  1. First, you can create an interpolator using interp1d. This means that the range of x values you initially give determines the range for which you will be able to call the interpolator. This interpolator can then be called for any value of x in this interval.
  2. You can perform the second integral even though you do not have a callable. The function simps only requires values and their location to provide an estimate of the integration process. Otherwise, you can perform the double integration in one call with dblquad.

Upvotes: 1

Related Questions