SymPy中TR8三角函数线性化功能不符合预期问题排查

问题背景 #

使用SymPy 1.13.3对含三角函数的表达式积分时，尝试通过sympy.simplify.fu.TR8将三角函数乘积线性化，理论上线性化前后积分结果应完全一致，但实际生成的图表存在明显差异：

• 未线性化：图表呈现平滑波动趋势
• 已线性化：图表起始段出现偏移和异常波动

简化代码如下：

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()

原因分析 #

1. TR8的局限性：TR8仅处理顶层的三角函数乘积，不会递归展开嵌套的乘积项。你的integrand是多层嵌套的乘积结构（u * u(s) * sum(...)），部分深层的三角函数乘积未被线性化，导致表达式与原表达式不等价。
2. 积分器的处理差异：SymPy积分器对未完全线性化的表达式和部分线性化的表达式可能采用不同的积分策略，引入了额外的常数项或计算误差，最终导致数值结果差异。

解决方案 #

方案1：使用expand_trig替代TR8 #

sp.expand_trig()会递归展开所有三角函数乘积为和差形式，比TR8更彻底：

# 替换原线性化代码
integrand_simple = sp.expand_trig(integrand)
integral_simple = sp.integrate(integrand_simple, (s, 0, t))

方案2：先展开表达式再应用TR8 #

先通过sp.expand()将所有乘积展开为单项之和，确保所有三角函数乘积都暴露在顶层，再用TR8处理：

# 替换原线性化代码
expanded_integrand = sp.expand(integrand)
integrand_simple = linearize_trigo_expr(expanded_integrand)
integral_simple = sp.integrate(integrand_simple, (s, 0, t))

方案3：验证表达式等价性 #

修改后可以通过以下代码验证线性化前后的表达式是否等价：

diff = sp.simplify(integrand - integrand_simple)
print(diff)  # 输出应为0或恒等于0的表达式

内容的提问来源于stack exchange，提问作者Thomas Guegamian