Integration problems in high harmonic generation model

Question

I’m running this code, to provide numerical solution for process of high harmonic generation.

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


X = np.zeros(100, dtype="complex")
E = 1  # амплитуда поля
w = 1  # частота поля
Up = 20  # np.power(E / 2 / w, 2)  # пондеромоторная энергия
Ip = 13.6  # потенциал ионизации атома
a = 2 * Ip  # параметр


def C(t):
    return np.sin(t) - 2 * 1j / a * np.cos(t) + 2 * 1j / (1 / a + 1j * t / 2) * \
           np.power(np.cos(t / 2) / a + 1j * np.sin(t / 2), 2)


def B(t):
    return -1 / 2 / np.power(1 / a + 1j * t / 2, 2) * np.power(np.cos(t / 2) / a + 1j * np.sin(t / 2), 2) + \
           np.cos(t / 2) / (1 / a + 1j * t / 2) * (np.cos(t / 2) / a + 1j * np.sin(t / 2)) - 0.5


def D(t):
    return -2 * B(t) - 1 + np.cos(t) + 1 / 4 / (1 / a + 1j * t / 2) / Up


def Fk(k, t):
    return (Ip + Up - k) * t - 2 * 1j * Up / a + 2 * 1j * Up / (1 / a + 1j * t / 2) * \
           np.power(np.cos(t / 2) / a + 1j * np.sin(t / 2), 2)


def int_func(k, t):
    return np.power(np.pi / 1 / a + 1j * t / 2, 1.5) * np.exp(-1j*Fk(k, t)) * \
           (B(t) * mp.besselj(k + 2, Up * C(t)) +
            1j * (B(t) * np.exp(1j * t) + D(t)) * mp.besselj(k + 1, Up * C(t)) +
            (B(t) + D(t) * np.exp(1j * t)) * mp.besselj(k, Up * C(t)) -
            1j * B(t) * np.exp(1j * t) * mp.besselj(k - 1, Up * C(t)))


def int_func_real(k, t):
    return np.real(int_func(k, t))


def int_func_imag(k, t):
    return np.imag(int_func(k, t))


def Integr(k):
    return quad(int_func_real, 0, np.inf, args=k)[0] + \
           1j * quad(int_func_imag, 0, np.inf, args=k)[0]

But it runs out with error:

IntegrationWarning: The maximum number of subdivisions (50) has been achieved.
If increasing the limit yields no improvement it is advised to analyze
the integrand in order to determine the difficulties. If the position of a
local difficulty can be determined (singularity, discontinuity) one will
probably gain from splitting up the interval and calling the integrator
on the subranges. Perhaps a special-purpose integrator should be used.
return quad(int_func_real, 0, np.inf, args=k)[0] + \

IntegrationWarning: The integral is probably divergent, or slowly convergent.
1j * quad(int_func_imag, 0, np.inf, args=k)[0]

How can I solve it?

I tried to change integration method to nquad but it didn’t help.

Leave a Comment Cancel reply