CODEWITHSUNDEEP
pythonscipy
AnSy
AnSy

Reputation: 123

Integration of oscillatory function does not converge in python

I want to evaluate the integral of a function using scipy.integrate.quad. Here is how the integrand looks like: Integrand of oscillatory function We can notice that most of the contribution from that integrand will come from .1 to 4 or 5 ish. For x = 10 and more, while the function is oscillatory (hard to tell on the picture), it is very small and keeps getting smaller.

Here is what the result of the integration from 0 to some upper bound looks like. The upper bound of the integration is on the x axis. Integration of integrand up to some upper bound Here, while the result seems to be stable for the first hundred x's, it is no longer the case after, even if I was expecting a straight line...

I am quite new to python, and I don't know what is the best course of action. My best guess for now would be to take some value smaller than 100 as an upper bound for my integral, and discard the other values because it is just a bad convergence from integrate.quad.

Edit: To plot the second graph, I used the scipy.integrate.quad function. However, if I just take the points that I generated to plot the integrand (first figure) and use this in scipy.integrate.simps and vary the maximum x up to which I integrate, I get consistent results.

Upvotes: 2

Views: 1485

Answers (1)

user6655984
user6655984

Reputation:

When the integrand has an important feature that is much smaller than the range of integration, it may get "overlooked" by the adaptive quad routine. In constrast, simps does not miss them if you use fine enough mesh, but may take longer to evaluate. I'll describe two ways to deal with this, the second being more practical.

points parameter

You can use the points parameter of quad to make sure this does not happen. Here is an example, where I integrate Gaussian function exp(-x**2) on the interval [-1000, 5000]. This function is localized near 0; pretty much all of it is in the interval [-5,5] so I include points=[-5, 5] to make sure this range is not overlooked. (The points are required to be within the range of integration, hence the appearance of if).

import numpy as np
from scipy.integrate import quad
import matplotlib.pyplot as plt

f = lambda x: np.exp(-x**2) 
numpoints = 1000
t = np.linspace(-1000, 5000, numpoints)
y = np.zeros((numpoints,))

for i in range(numpoints):
    y[i] = quad(f, -1000, t[i])[0]      # integration without points
plt.plot(t, y, 'r') 

for i in range(numpoints):            # using points if upper bound is above 5 
    y[i] = quad(f, -1000, t[i], points=[-5,5])[0] if t[i] > 5 else quad(f, -1000, t[i])[0]
plt.plot(t, y, 'b') 
plt.show()

The red curve is the output without points, the blue one is with points. The latter behaves as it should: rises from 0 to pi/2 and stays there.

results

Reusing previous calculation

A much more efficient approach to calculating antiderivative at multiple points is to use the previously computed value, adding to it the contribution of the interval that was not integrated over. 

y = np.zeros((numpoints,))
for i in range(1, numpoints):
    y[i] = y[i-1] + quad(f, t[i-1], t[i])[0]
plt.plot(t, y, 'g')

This has the same output as the blue curve above.

Upvotes: 6

Related Questions