[Sympy处理长表达式时给出错误的导数

Question

我有一个ODE系统，我想找到jacobi，但是sympy对派生给出错误的答案。

printing.init_printing(use_latex=True)
from sympy import *
import numpy as np
import matplotlib.pyplot as plt
import sympy as sp

N,x1,x2,x3,x4,y1,y2,r1,r2,r3,r4,eta1,eta2,eta3,eta4,R,c1,c2,c3,c4,a11,a12,a21,a22,a31,a32,a41,a42,b12,h,h11,h12,h21,h22,h31,h32,h41,h42,s1,s2,s3,s4,epsilon1,epsilon2,omega1,omega2,K11,K22,beta11,beta21,beta31,beta41,beta12,beta22,beta32,beta42,gamma12=sp.symbols('x1,x2,x3,x4,y1,y2,r1,r2,r3,r4,eta1,eta2,eta3,eta4,N,R,c1,c2,c3,c4,a11,a12,a21,a22,a31,a32,a41,a42,b12,h,h11,h12,h21,h22,h31,h32,h41,h42,s1,s2,s3,s4,epsilon1,epsilon2,omega1,omega2,K11,K22,beta11,beta21,beta3a,beta41,beta12,beta22,beta32,beta42,gamma12')
F1=R-c1*x1-c2*x2-c3*x3-c4*x4
F2=x1*(r1*(1-(eta1*x1+eta2*x2+eta3*x3+eta4*x4)/N)-(a11*y1)/(y1+a11*h11*x1)-(a12*y2)/(y2+a12*h12*x1))+s1
F3=x2*(r2*(1-(eta1*x1+eta2*x2+eta3*x3+eta4*x4)/N)-(a21*y1)/(y1+a21*h21*x2)-(a22*y2)/(y2+a22*h22*x2))+s2
F4=x3*(r3*(1-(eta1*x1+eta2*x2+eta3*x3+eta4*x4)/N)-(a31*y1)/(y1+a31*h31*x3)-(a32*y2)/(y2+a32*h32*x3))+s3
F5=x4*(r4*(1-(eta1*x1+eta2*x2+eta3*x3+eta4*x4)/N)-(a41*y1)/(y1+a42*h41*x4)-(a42*y2)/(y2+a42*h42*x4))+s4
F6=y1*(-epsilon1*(1+(y1+omega2*y2)/K22)-(b12*y2)/(y2+b12*h*y1)\
              +beta11*(a11*x1)/(y1+a11*h11*x1)\
              +beta21*(a21*x2)/(y1+a21*h21*x2)\
              +beta31*(a31*x3)/(y1+a31*h31*x3)\
              +beta41*(a41*x4)/(y1+a41*h41*x4))
F7=y2*(-epsilon2*(1+(omega1*y1+y2)/K11)-gamma12*(b12*y1)/(y2+b12*h*y1)\
              +beta12*(a12*x1)/(y2+a12*h12*x1)\
              +beta22*(a22*x2)/(y2+a22*h22*x2)\
              +beta32*(a32*x3)/(y2+a32*h32*x3)\
              +beta42*(a42*x4)/(y2+a42*h42*x4))               

dF2dN=sp.diff(F2,N)
print(dF2dN)

例如，当我尝试dF2dN时，我得到：

r2 * x2 *（N * y1 + eta2 * x2 + eta3 * x3 + eta4 * x4）/ x1 ** 2

正确的答案应该是：x1 * r1 *（eta1 * x1 + eta2 * x2 + eta3 * x3 + eta4 * x4）/ N ** 2）

有没有一种方法可以正确找到衍生产品，或者sympy无法做到？

Answer 1

请确保符号名称和变量名称匹配。您有

N,x1,...=sp.symbols('x1,x2,...')

因此有错误的（看起来）导数。理想情况下，避免使用像这样的质量符号定义，而是将它们分组，例如

x1, x2, x3, x4 = sp.symbols('x1, x2, x3, x4')
y1, y2 = sp.symbols('y1 y2')
N = sp.symbols('N')
etc.