How to numerically solve a non-Markovian integral equation?

Question

In the book Breuer, Heinz-Peter, and Francesco Petruccione. The Theory of Open Quantum Systems. Oxford: Clarendon Press, 2009., there is equation(10.17)

$$ \frac{\mathrm{d}}{\mathrm{d}t}c_1(t) = - \int_0^t \mathrm{d}t_1 f(t-t_1) c_1(t_1) $$ where $$ f(t) = \frac{1}{2} \gamma_0\lambda e^{-\lambda |t|} $$ $\gamma_0$ and $\lambda$ are parameters, for example, we can set $\gamma_0=1, \lambda=10$. For this $f(t)$, we can solve it analytically, the solution is $$ c_1(t) = c_1(0) e^{-\lambda t/2} \left[\cosh\left(\frac{dt}{2}\right) + \frac{\lambda}{d} \sinh\left(\frac{dt}{2}\right) \right] $$ where $d = \sqrt{\lambda^2 - 2\gamma_0\lambda}$.

My question is, how to numerically solve (for example, using python function scipy.integrate.solve_ivp, which is a function solving OEDs by Runge-Kutta method or other methods) this equation?

I want to numerically solve it, because for other form of $f(t)$, we may not solve it analytically.

Answer 1

As a first approach you can do something like the following to match the analytical solution

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

gamma_0 = 1
lambda_ = 10

d = np.sqrt(lambda_**2 - 2 * gamma_0 * lambda_)

c1_0 = 1.0
t_span = (0, 10.0)
t_eval = np.linspace(t_span[0], t_span[1], 500)

def create_differential_eq_function(gamma_0, lambda_):
    c1_history = []
    t_history = []

    def f(t):
        return 0.5 * gamma_0 * lambda_ * np.exp(-lambda_ * np.abs(t))

    def differential_eq(t, c1):
        c1_history.append(c1[0])
        t_history.append(t)

        if len(t_history) < 2:
            return 0

        # Calculate the integral using trapezoidal rule
        integral_value = 0
        for i in range(1, len(t_history)):
            if t_history[i] <= t:
                dt = t_history[i] - t_history[i - 1]
                average_c1 = 0.5 * (c1_history[i] + c1_history[i - 1])
                integral_value += f(t - t_history[i]) * average_c1 * dt
        return -integral_value

    return differential_eq

differential_eq_function = create_differential_eq_function(gamma_0, lambda_)

solution = solve_ivp(differential_eq_function, t_span, [c1_0], t_eval=t_eval, method='RK23',rtol=1e-6)

solution.t, solution.y[0]  

And you can compare to the analytical solution using

def analytical_solution(t, lambda_, d):
    return np.exp(-lambda_ * t / 2) * (np.cosh(d * t / 2) + (lambda_ / d) * np.sinh(d * t / 2))

analytical_values = analytical_solution(solution.t, lambda_, d)

plt.plot(solution.t, solution.y[0], label='Numerical Solution', color='blue')
plt.plot(solution.t, analytical_values, label='Analytical Solution', linestyle='--', color='red')
plt.xlabel('Time')
plt.ylabel('$c_1(t)$')
plt.title('Comparison of Numerical and Analytical Solutions')
plt.legend()
plt.grid(True)
plt.show()

Which will show