SymPy TR8 Trigonometric Linearization Not Working as Expected

I'm using SymPy version 1.13.3 to integrate an expression involving trigonometric functions. My goal is to linearize trigonometric products using the TR8 function from sympy.simplify.fu. However, the results with and without linearization should be identical, but they are not.

Here's a simplified version of my code:

import sympy as sp
import numpy as np
import matplotlib.pyplot as plt
from sympy.simplify.fu import TR8 as linearize_trigo_expr
from sympy import cos, sin

# Definitions
JMax = 3
nb_vals_t = 100
t_vals = np.linspace(0, 0.2, nb_vals_t)
s, t = sp.symbols('s t')
u = -12.33203125*sin(t) + 5.6064453125*sin(4*t) - 2.2708740234375*sin(9*t) + sin(12*t) - 0.103253126144409*sin(16*t) - 1.39863395690918*cos(t) + 1.81103515625*cos(4*t) - 0.4464111328125*cos(9*t) + 0.0336227416992188*cos(16*t)
integrand = u * (u.subs(t, s) * sum(sp.sin(j**2 * (t - s)) / j**4 for j in range(1, JMax + 1)))

# Without trigonometric linearization
integral = sp.integrate(integrand, (s, 0, t))
integral_values = []
for t_val in t_vals:
    integral_values.append(integral.subs(t, t_val))

plt.plot(t_vals, integral_values)
plt.xlabel('t')
plt.ylabel('Integral')
plt.title("Test without trigonometric linearization")
plt.grid(True)
plt.show()

# With trigonometric linearization
integrand_simple = linearize_trigo_expr(integrand)
integral_simple = sp.integrate(integrand_simple, (s, 0, t))

integral_simple_values = []
for t_val in t_vals:
    integral_simple_values.append(integral_simple.subs(t, t_val))

plt.plot(t_vals, integral_simple_values)
plt.xlabel('t')
plt.ylabel('Integral')
plt.title("Test with trigonometric linearization")
plt.grid(True)
plt.show()

Both graphs should be identical, but they differ. Does anyone know why the trigonometric linearization isn't working as expected, and if there is a solution to make it work properly?

NB: this is probably not a MWE, but I have tried to find a simple example showing the mistakes that the TR8 function does.

Here are the graphs:

The problem is that you are using float numbers with symbolic integration, which sometimes leads to wrong results. When doing symbolic integration, always use exact numbers. Use sympy's nsimplify to convert floats to rational.

integrand_2 = integrand.nsimplify()
integrand_3 = linearize_trigo_expr(integrand_2)
r2 = integrand_2.integrate((s, 0, t))
r3 = integrand_3.integrate((s, 0, t))

sp.plot(r2, r3, (t, 0, 0.2))

The results are the same.