3

I am new to stackoverflow and also quite new to Python. So, I hope to ask my question in an appropriate manner. I am running a Python code similar to this minimal example with an example function that is a product of a lorentzian with a cosinus that I want to numerically integrate:

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

#minimal example:
omega_loc = 15
gamma = 5

def Lorentzian(w):
    #print(w)
    return (w**3)/((w/omega_loc) + 1)**2*(gamma/2)/((w-omega_loc)**2+(gamma/2)**2)

def intRe(t):
    return quad(lambda w: w**(-2)*Lorentzian(w)*(1-np.cos(w*t)),0,np.inf,limit=10000)[0]


plt.figure(1)
plot_range = np.linspace(0,100,1000)
plt.plot(plot_range, [intRe(t) for t in plot_range])

Independent on the upper limit of the integration I never get the code to run and to give me a result. When I enable the #print(w) line it seems like the code just keeps on probing the integral at random different values of w in an infinite loop (?). Also the console gives me a detection of a roundoff error. Is there a different way for numerical integration in Python that is better suited for this kind of function than the quad function or did I do a more fundamental error?

Frederik
  • 35
  • 4
  • 2
    If you integrate to `100` instead of `np.inf`, the integration result will be much quicker and it'll not be that far off from the `np.inf` result. Not sure about the reason but it might be the detection of the infinity conversion is not working well for your function. I'll leave this to the experts. – hesham_EE Jan 10 '22 at 18:42

1 Answers1

2

Observations

  1. Close to zero (1 - cos(w*t)) / w**2 tends to 0/0. We can take the taylor expansion t**2(1/2 - (w*t)**2/24).

  2. Going to the infinity the Lorentzian is a constant and the cosine term will cause the output to oscillate indefinitely, the integral can be approximated by multiplying that term by a slowly decreasing term.

  3. You are using a linearly spaced scale with many points. It is easier to visualize with w in log scale.

The plot looks like this before damping the cosine term

enter image description here

I introduced two parameters to tune the attenuation of the oscilations

def cosinus_term(w, t, damping=1e4*omega_loc):
    return np.where(abs(w*t) < 1e-6,  t**2*(0.5 - (w*t)**2/24.0), (1-np.exp(-abs(w/damping))*np.cos(w*t))/w**2)
def intRe(t, damping=1e4*omega_loc):
    return quad(lambda w: cosinus_term(w, t)*Lorentzian(w),0,np.inf,limit=10000)[0]

Plotting with the following code

plt.figure(1)
plot_range = np.logspace(-3,3,100)
plt.plot(plot_range, [intRe(t, 1e2*omega_loc) for t in plot_range])
plt.plot(plot_range, [intRe(t, 1e3*omega_loc) for t in plot_range])
plt.xscale('log')

It runs in less than 3 minutes here, and the two results are close to each other, specially for large w, suggesting that the damping doesn't affect too much the result.

enter image description here

Bob
  • 13,867
  • 1
  • 5
  • 27