Sympy - dsolve 函数计算能力的限制？
Sympy - limitations in the computation ability of dsolve function?
我正在尝试使用 Sympy 求解四阶微分方程组。我使用的方程式如图所示，并写在下面的代码中：
latex_equations :
from sympy import *
x = symbols('x')
EI1,EI2,EI3,a1,a2,a3,Qh,Mecc = symbols('EI1 EI2 EI3 a1 a2 a3 Qh Mecc')
u1,u2,u3 = symbols('u1 u2 u3', cls=Function)
eq = (Eq(EI1*diff(u1(x),x,x,x,x)+EI2*diff(u2(x),x,x,x,x)+EI3*diff(u3(x),x,x,x,x), Qh),Eq(a1*EI1*diff(u1(x),x,x,x,x)+a2*EI2*diff(u2(x),x,x,x,x)+a3*EI3*diff(u3(x),x,x,x,x),Mecc),Eq((u1(x)+u3(x))/2,u2(x)))
dsolve(eq)
我收到以下错误：
AttributeError Traceback (most recent call last)
<ipython-input-1-63c42d2751be> in <module>
6 eq = (Eq(EI1*diff(u1(x),x,x,x,x)+EI2*diff(u2(x),x,x,x,x)+EI3*diff(u3(x),x,x,x,x), Qh),Eq(a1*EI1*diff(u1(x),x,x,x,x)+a2*EI2*diff(u2(x),x,x,x,x)+a3*EI3*diff(u3(x),x,x,x,x),Mecc),Eq((u1(x)+u3(x))/2,u2(x)))
7
----> 8 dsolve(eq)
~\Anaconda3\lib\site-packages\sympy\solvers\ode.py in dsolve(eq, func, hint, simplify, ics, xi, eta, x0, n, **kwargs)
583 """
584 if iterable(eq):
--> 585 match = classify_sysode(eq, func)
586 eq = match['eq']
587 order = match['order']
~\Anaconda3\lib\site-packages\sympy\solvers\ode.py in classify_sysode(eq, funcs, **kwargs)
1528 if isinstance(func, list):
1529 for func_elem in func:
-> 1530 if len(func_elem.args) != 1:
1531 raise ValueError("dsolve() and classify_sysode() work with "
1532 "functions of one variable only, not %s" % func)
AttributeError: 'list' object has no attribute 'args'
我尝试使用 dsolve 求解更简单的方程组，结果很好：
from sympy import *
x = symbols('x')
EI1,EI2 = symbols('EI1 EI2')
u1,u2 = symbols('u1 u2', cls=Function)
eq = (Eq(EI1*diff(u1(x),x), 12*x*u1(x) + 8*u2(x)), Eq(EI2*diff(u2(x),x), 21*u1(x) + 7*x*u2(x)))
dsolve(eq)
这两种情况我用的格式是一样的，只是一个解决了，一个失败了。我知道第一个方程组有解，因为我已经在Maple中解过了。
是我的代码出错了，还是 Sympy dsolve 根本无法求解如此复杂的方程组？在 dsolve 无法再求解之前，方程组的复杂程度是否有限制？非常感谢任何对此问题的帮助或见解。
谢谢！
SymPy 的 ODE 系统代码需要大量工作，目前无法处理这种系统。
首先，第三个方程是一个代数方程（这意味着这实际上是一个 DAE 而不是 ODE 系统）但这不是一个大问题：只需对它进行 4 次微分。那么 SymPy 应该真的能够解决生成的系统，但不幸的是不能，因为它还没有实现。不过，您可以自己手动完成：
# Your code first:
from sympy import *
x = symbols('x')
EI1,EI2,EI3,a1,a2,a3,Qh,Mecc = symbols('EI1 EI2 EI3 a1 a2 a3 Qh Mecc')
u1,u2,u3 = symbols('u1 u2 u3', cls=Function)
eq = (Eq(EI1*diff(u1(x),x,x,x,x)+EI2*diff(u2(x),x,x,x,x)+EI3*diff(u3(x),x,x,x,x), Qh),Eq(a1*EI1*diff(u1(x),x,x,x,x)+a2*EI2*diff(u2(x),x,x,x,x)+a3*EI3*diff(u3(x),x,x,x,x),Mecc),Eq((u1(x)+u3(x))/2,u2(x)))
# Differentiate eq3 4 times:
eq1, eq2, eq3 = eq
eq3 = Eq(eq3.lhs.diff(x, 4).doit(), eq3.rhs.diff(x, 4).doit())
eq = eq1, eq2, eq3
# Isolate equations for u1(x).diff(x, 4) etc
derivs = [u(x).diff(x, 4) for u in (u1, u2, u3)]
(sol,) = solve(eq, derivs, dict=True)
eq = [Eq(du, sol[du]) for du in derivs]
# Each can be solved separately:
for e in eq:
pprint(dsolve(e))
这给出了解决方案：
4
2 3 x ⋅(EI₂⋅Mecc - EI₂⋅Qh⋅a₂ + 2⋅EI₃⋅Mecc - 2⋅EI₃⋅Qh⋅a₃)
u₁(x) = C₁ + C₂⋅x + C₃⋅x + C₄⋅x + ────────────────────────────────────────────────────────────────────────────────────
24⋅(EI₁⋅EI₂⋅a₁ - EI₁⋅EI₂⋅a₂ + 2⋅EI₁⋅EI₃⋅a₁ - 2⋅EI₁⋅EI₃⋅a₃ + EI₂⋅EI₃⋅a₂ - EI₂⋅EI₃⋅a₃)
4
2 3 x ⋅(-EI₁⋅Mecc + EI₁⋅Qh⋅a₁ + EI₃⋅Mecc - EI₃⋅Qh⋅a₃)
u₂(x) = C₁ + C₂⋅x + C₃⋅x + C₄⋅x + ────────────────────────────────────────────────────────────────────────────────────
24⋅(EI₁⋅EI₂⋅a₁ - EI₁⋅EI₂⋅a₂ + 2⋅EI₁⋅EI₃⋅a₁ - 2⋅EI₁⋅EI₃⋅a₃ + EI₂⋅EI₃⋅a₂ - EI₂⋅EI₃⋅a₃)
4
2 3 x ⋅(-2⋅EI₁⋅Mecc + 2⋅EI₁⋅Qh⋅a₁ - EI₂⋅Mecc + EI₂⋅Qh⋅a₂)
u₃(x) = C₁ + C₂⋅x + C₃⋅x + C₄⋅x + ────────────────────────────────────────────────────────────────────────────────────
24⋅(EI₁⋅EI₂⋅a₁ - EI₁⋅EI₂⋅a₂ + 2⋅EI₁⋅EI₃⋅a₁ - 2⋅EI₁⋅EI₃⋅a₃ + EI₂⋅EI₃⋅a₂ - EI₂⋅EI₃⋅a₃)
请注意，尽管每个方程的积分常数应该不同，因此第一个方程中的 C1 与第二个方程中的 C2 不同。 ODE 系统应该有 12 个不同的积分常数。实际上，由于这是一个DAE系统，因此应该有8个独立常数，因为三个因变量之间存在代数关系。
编辑：在 SymPy 1.8 中，dsolve 函数现在可以解决这个问题，前提是第三个方程微分：
In [67]: # Your code first:
...: from sympy import *
...: x = symbols('x')
...: EI1,EI2,EI3,a1,a2,a3,Qh,Mecc = symbols('EI1 EI2 EI3 a1 a2 a3 Qh Mecc')
...: u1,u2,u3 = symbols('u1 u2 u3', cls=Function)
...:
...: eq = (Eq(EI1*diff(u1(x),x,x,x,x)+EI2*diff(u2(x),x,x,x,x)+EI3*diff(u3(x),x,x,x,x), Qh),Eq(a1*EI1*diff(u1(x),x,x,x,x)+a2*EI2*diff(u2(x),x,x,x,x)+a3*EI3
...: *diff(u3(x),x,x,x,x),Mecc),Eq((u1(x)+u3(x))/2,u2(x)))
...:
...: # Differentiate eq3 4 times:
...: eq1, eq2, eq3 = eq
...: eq3 = Eq(eq3.lhs.diff(x, 4).doit(), eq3.rhs.diff(x, 4).doit())
...: eq = eq1, eq2, eq3
In [68]: dsolve(eq)
Out[68]:
⎡ 2 3 4 2 3 4
⎢ C₃⋅x C₄⋅x x ⋅(EI₂⋅(Mecc - Qh⋅a₂) + 2⋅EI₃⋅(Mecc - Qh⋅a₃)) C₇⋅x C₈⋅x x ⋅(EI₁⋅(M
⎢u₁(x) = C₁ + C₂⋅x + ───── + ───── + ──────────────────────────────────────────────────────────────, u₂(x) = C₅ + C₆⋅x + ───── + ───── - ───────────────────
⎣ 2 6 24⋅(EI₁⋅(EI₂⋅(a₁ - a₂) + 2⋅EI₃⋅(a₁ - a₃)) + EI₂⋅EI₃⋅(a₂ - a₃)) 2 6 24⋅(EI₁⋅(EI₂⋅(a₁ -
2 3 4 ⎤
ecc - Qh⋅a₁) - EI₃⋅(Mecc - Qh⋅a₃)) C₁₁⋅x C₁₂⋅x x ⋅(2⋅EI₁⋅(Mecc - Qh⋅a₁) + EI₂⋅(Mecc - Qh⋅a₂)) ⎥
───────────────────────────────────────────, u₃(x) = C₁₀⋅x + ────── + ────── + C₉ - ──────────────────────────────────────────────────────────────⎥
a₂) + 2⋅EI₃⋅(a₁ - a₃)) + EI₂⋅EI₃⋅(a₂ - a₃)) 2 6 24⋅(EI₁⋅(EI₂⋅(a₁ - a₂) + 2⋅EI₃⋅(a₁ - a₃)) + EI₂⋅EI₃⋅(a₂ - a₃))⎦
我正在尝试使用 Sympy 求解四阶微分方程组。我使用的方程式如图所示，并写在下面的代码中：
latex_equations :
from sympy import *
x = symbols('x')
EI1,EI2,EI3,a1,a2,a3,Qh,Mecc = symbols('EI1 EI2 EI3 a1 a2 a3 Qh Mecc')
u1,u2,u3 = symbols('u1 u2 u3', cls=Function)
eq = (Eq(EI1*diff(u1(x),x,x,x,x)+EI2*diff(u2(x),x,x,x,x)+EI3*diff(u3(x),x,x,x,x), Qh),Eq(a1*EI1*diff(u1(x),x,x,x,x)+a2*EI2*diff(u2(x),x,x,x,x)+a3*EI3*diff(u3(x),x,x,x,x),Mecc),Eq((u1(x)+u3(x))/2,u2(x)))
dsolve(eq)
我收到以下错误：
AttributeError Traceback (most recent call last)
<ipython-input-1-63c42d2751be> in <module>
6 eq = (Eq(EI1*diff(u1(x),x,x,x,x)+EI2*diff(u2(x),x,x,x,x)+EI3*diff(u3(x),x,x,x,x), Qh),Eq(a1*EI1*diff(u1(x),x,x,x,x)+a2*EI2*diff(u2(x),x,x,x,x)+a3*EI3*diff(u3(x),x,x,x,x),Mecc),Eq((u1(x)+u3(x))/2,u2(x)))
7
----> 8 dsolve(eq)
~\Anaconda3\lib\site-packages\sympy\solvers\ode.py in dsolve(eq, func, hint, simplify, ics, xi, eta, x0, n, **kwargs)
583 """
584 if iterable(eq):
--> 585 match = classify_sysode(eq, func)
586 eq = match['eq']
587 order = match['order']
~\Anaconda3\lib\site-packages\sympy\solvers\ode.py in classify_sysode(eq, funcs, **kwargs)
1528 if isinstance(func, list):
1529 for func_elem in func:
-> 1530 if len(func_elem.args) != 1:
1531 raise ValueError("dsolve() and classify_sysode() work with "
1532 "functions of one variable only, not %s" % func)
AttributeError: 'list' object has no attribute 'args'
我尝试使用 dsolve 求解更简单的方程组，结果很好：
from sympy import *
x = symbols('x')
EI1,EI2 = symbols('EI1 EI2')
u1,u2 = symbols('u1 u2', cls=Function)
eq = (Eq(EI1*diff(u1(x),x), 12*x*u1(x) + 8*u2(x)), Eq(EI2*diff(u2(x),x), 21*u1(x) + 7*x*u2(x)))
dsolve(eq)
这两种情况我用的格式是一样的，只是一个解决了，一个失败了。我知道第一个方程组有解，因为我已经在Maple中解过了。
是我的代码出错了，还是 Sympy dsolve 根本无法求解如此复杂的方程组？在 dsolve 无法再求解之前，方程组的复杂程度是否有限制？非常感谢任何对此问题的帮助或见解。
谢谢！
SymPy 的 ODE 系统代码需要大量工作，目前无法处理这种系统。
首先，第三个方程是一个代数方程（这意味着这实际上是一个 DAE 而不是 ODE 系统）但这不是一个大问题：只需对它进行 4 次微分。那么 SymPy 应该真的能够解决生成的系统，但不幸的是不能，因为它还没有实现。不过，您可以自己手动完成：
# Your code first:
from sympy import *
x = symbols('x')
EI1,EI2,EI3,a1,a2,a3,Qh,Mecc = symbols('EI1 EI2 EI3 a1 a2 a3 Qh Mecc')
u1,u2,u3 = symbols('u1 u2 u3', cls=Function)
eq = (Eq(EI1*diff(u1(x),x,x,x,x)+EI2*diff(u2(x),x,x,x,x)+EI3*diff(u3(x),x,x,x,x), Qh),Eq(a1*EI1*diff(u1(x),x,x,x,x)+a2*EI2*diff(u2(x),x,x,x,x)+a3*EI3*diff(u3(x),x,x,x,x),Mecc),Eq((u1(x)+u3(x))/2,u2(x)))
# Differentiate eq3 4 times:
eq1, eq2, eq3 = eq
eq3 = Eq(eq3.lhs.diff(x, 4).doit(), eq3.rhs.diff(x, 4).doit())
eq = eq1, eq2, eq3
# Isolate equations for u1(x).diff(x, 4) etc
derivs = [u(x).diff(x, 4) for u in (u1, u2, u3)]
(sol,) = solve(eq, derivs, dict=True)
eq = [Eq(du, sol[du]) for du in derivs]
# Each can be solved separately:
for e in eq:
pprint(dsolve(e))
这给出了解决方案：
4
2 3 x ⋅(EI₂⋅Mecc - EI₂⋅Qh⋅a₂ + 2⋅EI₃⋅Mecc - 2⋅EI₃⋅Qh⋅a₃)
u₁(x) = C₁ + C₂⋅x + C₃⋅x + C₄⋅x + ────────────────────────────────────────────────────────────────────────────────────
24⋅(EI₁⋅EI₂⋅a₁ - EI₁⋅EI₂⋅a₂ + 2⋅EI₁⋅EI₃⋅a₁ - 2⋅EI₁⋅EI₃⋅a₃ + EI₂⋅EI₃⋅a₂ - EI₂⋅EI₃⋅a₃)
4
2 3 x ⋅(-EI₁⋅Mecc + EI₁⋅Qh⋅a₁ + EI₃⋅Mecc - EI₃⋅Qh⋅a₃)
u₂(x) = C₁ + C₂⋅x + C₃⋅x + C₄⋅x + ────────────────────────────────────────────────────────────────────────────────────
24⋅(EI₁⋅EI₂⋅a₁ - EI₁⋅EI₂⋅a₂ + 2⋅EI₁⋅EI₃⋅a₁ - 2⋅EI₁⋅EI₃⋅a₃ + EI₂⋅EI₃⋅a₂ - EI₂⋅EI₃⋅a₃)
4
2 3 x ⋅(-2⋅EI₁⋅Mecc + 2⋅EI₁⋅Qh⋅a₁ - EI₂⋅Mecc + EI₂⋅Qh⋅a₂)
u₃(x) = C₁ + C₂⋅x + C₃⋅x + C₄⋅x + ────────────────────────────────────────────────────────────────────────────────────
24⋅(EI₁⋅EI₂⋅a₁ - EI₁⋅EI₂⋅a₂ + 2⋅EI₁⋅EI₃⋅a₁ - 2⋅EI₁⋅EI₃⋅a₃ + EI₂⋅EI₃⋅a₂ - EI₂⋅EI₃⋅a₃)
请注意，尽管每个方程的积分常数应该不同，因此第一个方程中的 C1 与第二个方程中的 C2 不同。 ODE 系统应该有 12 个不同的积分常数。实际上，由于这是一个DAE系统，因此应该有8个独立常数，因为三个因变量之间存在代数关系。
编辑：在 SymPy 1.8 中，dsolve 函数现在可以解决这个问题，前提是第三个方程微分：
In [67]: # Your code first:
...: from sympy import *
...: x = symbols('x')
...: EI1,EI2,EI3,a1,a2,a3,Qh,Mecc = symbols('EI1 EI2 EI3 a1 a2 a3 Qh Mecc')
...: u1,u2,u3 = symbols('u1 u2 u3', cls=Function)
...:
...: eq = (Eq(EI1*diff(u1(x),x,x,x,x)+EI2*diff(u2(x),x,x,x,x)+EI3*diff(u3(x),x,x,x,x), Qh),Eq(a1*EI1*diff(u1(x),x,x,x,x)+a2*EI2*diff(u2(x),x,x,x,x)+a3*EI3
...: *diff(u3(x),x,x,x,x),Mecc),Eq((u1(x)+u3(x))/2,u2(x)))
...:
...: # Differentiate eq3 4 times:
...: eq1, eq2, eq3 = eq
...: eq3 = Eq(eq3.lhs.diff(x, 4).doit(), eq3.rhs.diff(x, 4).doit())
...: eq = eq1, eq2, eq3
In [68]: dsolve(eq)
Out[68]:
⎡ 2 3 4 2 3 4
⎢ C₃⋅x C₄⋅x x ⋅(EI₂⋅(Mecc - Qh⋅a₂) + 2⋅EI₃⋅(Mecc - Qh⋅a₃)) C₇⋅x C₈⋅x x ⋅(EI₁⋅(M
⎢u₁(x) = C₁ + C₂⋅x + ───── + ───── + ──────────────────────────────────────────────────────────────, u₂(x) = C₅ + C₆⋅x + ───── + ───── - ───────────────────
⎣ 2 6 24⋅(EI₁⋅(EI₂⋅(a₁ - a₂) + 2⋅EI₃⋅(a₁ - a₃)) + EI₂⋅EI₃⋅(a₂ - a₃)) 2 6 24⋅(EI₁⋅(EI₂⋅(a₁ -
2 3 4 ⎤
ecc - Qh⋅a₁) - EI₃⋅(Mecc - Qh⋅a₃)) C₁₁⋅x C₁₂⋅x x ⋅(2⋅EI₁⋅(Mecc - Qh⋅a₁) + EI₂⋅(Mecc - Qh⋅a₂)) ⎥
───────────────────────────────────────────, u₃(x) = C₁₀⋅x + ────── + ────── + C₉ - ──────────────────────────────────────────────────────────────⎥
a₂) + 2⋅EI₃⋅(a₁ - a₃)) + EI₂⋅EI₃⋅(a₂ - a₃)) 2 6 24⋅(EI₁⋅(EI₂⋅(a₁ - a₂) + 2⋅EI₃⋅(a₁ - a₃)) + EI₂⋅EI₃⋅(a₂ - a₃))⎦