SymPy中如何对表达式的每个因子单独进行泰勒展开？

问题描述 #

使用SymPy对形如P(x,y)*Q(x,y)*.../(K(x,y)*L(x,y)*....)的复杂分式做x和y的泰勒展开时，程序无法正常运行（其中所有大写字母均为x和y的多项式）。计划先对每个多项式因子单独做泰勒展开，完成乘除运算后再进行二次泰勒展开，以此解决程序运行异常的问题。

待处理表达式 #

D*(-c3*(k**3*p*((sin_theta*cos_phi)*sin_psi + cos_psi*cos_theta) + 2*k**2*p**2 + k*p**3*((sin_theta*cos_phi)*sin_psi + cos_psi*cos_theta)) + (c1*(k**2 + 2*k*p*((sin_theta*cos_phi)*sin_psi + cos_psi*cos_theta) + p**2) + 2*c2*k*p*((sin_theta*cos_phi)*sin_psi + cos_psi*cos_theta))*(k**2 + 2*k*p*((sin_theta*cos_phi)*sin_psi + cos_psi*cos_theta) + p**2))*(-c3*(-2*cos_psi*k**2*p*q - cos_theta*k**3*q - cos_theta*k*q*(-2*cos_psi*p*q + p**2 + q**2) + k**3*p*((sin_theta*cos_phi)*sin_psi + cos_psi*cos_theta) + k**2*p**2 + k**2*q**2 + k**2*(-2*cos_psi*p*q + p**2 + q**2) + k*p*((sin_theta*cos_phi)*sin_psi + cos_psi*cos_theta)*(-2*cos_psi*p*q + p**2 + q**2))/2 + (c1*(-2*cos_psi*p*q - 2*cos_theta*k*q + k**2 + 2*k*p*((sin_theta*cos_phi)*sin_psi + cos_psi*cos_theta) + p**2 + q**2) - 2*c2*cos_theta*k*q + 2*c2*k*p*((sin_theta*cos_phi)*sin_psi + cos_psi*cos_theta))*(-2*cos_psi*p*q - 2*cos_theta*k*q + k**2 + 2*k*p*((sin_theta*cos_phi)*sin_psi + cos_psi*cos_theta) + p**2 + q**2)/2)*(2*k**4*kappa + kappa*(k**2 + 2*k*p*((sin_theta*cos_phi)*sin_psi + cos_psi*cos_theta) + p**2)**2 + kappa*(-2*cos_psi*p*q - 2*cos_theta*k*q + k**2 + 2*k*p*((sin_theta*cos_phi)*sin_psi + cos_psi*cos_theta) + p**2 + q**2)**2)*(c1*q**2/2 + c2*cos_psi*p*q + c2*cos_theta*k*q - c2*k**2 - 2*c2*k*p*((sin_theta*cos_phi)*sin_psi + cos_psi*cos_theta) - c2*p**2 + c3*cos_psi**2*p**2 + 2*c3*cos_psi*cos_theta*k*p - c3*cos_psi*p*q/2 + c3*cos_theta**2*k**2 - c3*cos_theta*k*q/2 - c3*k**2/2 - c3*k*p*((sin_theta*cos_phi)*sin_psi + cos_psi*cos_theta) - c3*p**2/2)/(k**2*kappa*(k**4*kappa + kappa*(k**2 + 2*k*p*((sin_theta*cos_phi)*sin_psi + cos_psi*cos_theta) + p**2)**2)*(k**4*kappa + kappa*(-2*cos_psi*p*q - 2*cos_theta*k*q + k**2 + 2*k*p*((sin_theta*cos_phi)*sin_psi + cos_psi*cos_theta) + p**2 + q**2)**2)*(kappa*(k**2 + 2*k*p*((sin_theta*cos_phi)*sin_psi + cos_psi*cos_theta) + p**2)**2 + kappa*(-2*cos_psi*p*q - 2*cos_theta*k*q + k**2 + 2*k*p*((sin_theta*cos_phi)*sin_psi + cos_psi*cos_theta) + p**2 + q**2)**2))

可通过from IPython.display import display调用display函数展示该表达式。

解决方案思路与代码示例 #

步骤1：拆分分子分母的多项式因子 #

将原表达式拆分为分子的各个多项式因子和分母的各个多项式因子，分别定义为SymPy符号表达式。

步骤2：对每个因子单独做泰勒展开 #

使用SymPy的series函数，对每个多项式因子在指定展开点（例如(0,0)）展开到二次项，并用.removeO()截断高阶项。

步骤3：组合展开后的因子为新分式 #

将分子展开后的因子相乘，分母展开后的因子相乘，得到简化后的分式表达式。

步骤4：对组合后的分式做二次泰勒展开 #

对简化后的分式再次执行泰勒展开，得到最终结果。

示例代码框架：

import sympy as sp

# 定义所有符号变量
D, c1, c2, c3, k, p, q, kappa, sin_theta, cos_phi, sin_psi, cos_psi, cos_theta = sp.symbols('D c1 c2 c3 k p q kappa sin_theta cos_phi sin_psi cos_psi cos_theta')
# 替换为你实际需要展开的变量，比如x,y
expand_var1, expand_var2 = sp.symbols('x y')  

# 拆分并定义各个多项式因子，替换为你实际的因子表达式
numerator_factor1 = ...
numerator_factor2 = ...
denominator_factor1 = ...
denominator_factor2 = ...

# 设置展开点与展开阶数（阶数设为3是为了保留到二次项）
expand_point1, expand_point2 = 0, 0
expand_order = 3

# 对每个因子做二次泰勒展开
nf1_expanded = sp.series(numerator_factor1, expand_var1, expand_point1, expand_order).removeO()
nf1_expanded = sp.series(nf1_expanded, expand_var2, expand_point2, expand_order).removeO()

nf2_expanded = sp.series(numerator_factor2, expand_var1, expand_point1, expand_order).removeO()
nf2_expanded = sp.series(nf2_expanded, expand_var2, expand_point2, expand_order).removeO()

df1_expanded = sp.series(denominator_factor1, expand_var1, expand_point1, expand_order).removeO()
df1_expanded = sp.series(df1_expanded, expand_var2, expand_point2, expand_order).removeO()

df2_expanded = sp.series(denominator_factor2, expand_var1, expand_point1, expand_order).removeO()
df2_expanded = sp.series(df2_expanded, expand_var2, expand_point2, expand_order).removeO()

# 组合成新分式并展开
new_numerator = nf1_expanded * nf2_expanded
new_denominator = df1_expanded * df2_expanded
simplified_expr = new_numerator / new_denominator

final_expansion = sp.series(simplified_expr, expand_var1, expand_point1, expand_order).removeO()
final_expansion = sp.series(final_expansion, expand_var2, expand_point2, expand_order).removeO()

# 展示结果
from IPython.display import display
display(final_expansion)

注意事项 #

• 拆分因子时必须保证每个因子都是纯多项式，不能包含分式，否则单独展开会报错。
• 若实际展开变量不是x和y，需替换为你需要的变量（比如k、p等）。
• 展开阶数设置为3是因为SymPy的series函数第三个参数是展开的最高阶数+1，要保留二次项则设为3。

内容的提问来源于stack exchange，提问作者Nikos Pap