如何解决如何使用 sympy
我试图通过解析求解这个微分方程来超越数值分析课（不要担心求解；我只是想让你们知道我来自哪里）：
这里，v 是因变量，t 是我们的自变量，所以我想明确地将 v 求解为 t 的函数。给出一个初始条件：v(0) = 0。
这是我得到的微分方程的解。我确定这是正确的，但如果您想亲自尝试，请成为我的客人。
desmos 中这个系列的情节是这样的：
进行分析计算的唯一方法是使用幂级数展开式。所以基本上，左边的表达式是一个很大的多项式作为 f(v)，但这不是我想要的。
我希望能够使用 matplotlib 绘制此函数，并且我希望能够获得任意 t 值下的速度。
我对 sympy 很陌生，所以我还不熟悉所有内容，但这里有一些代码至少可以得到等式左侧的表达式：
import sympy as sp
if __name__ == '__main__':
drag_coefficient = 1.85
g = 9.81
m = 100
rho = 1.2
A = 30
a = (drag_coefficient * rho * A)/(2*m)
v = sp.Symbol('v')
f_ana = (a**0 * v**(2*0 + 1))/(g**(0 + 1) * (2*0 + 1))
for i in range(1,101):
f_ana = f_ana + (a**i * v**(2*i + 1))/(g**(i + 1) * (2*i + 1))
print(f_ana)
f_ana 的输出符合预期：
6.0625809125107e-151*v**201 + 1.80395335070212e-149*v**199 + 5.36830184251333e-148*v**197 + 1.59769293782798e-146*v**195 + 4.75549130283346e-145*v**193 + 1.41561158253442e-143*v**191 .....
然后我尝试通过执行以下操作来计算时间 (t) 的速度：
sp.solve(f_ana + 1,v) # calculating at t = 1s
程序不会完成。我假设这是因为这个表达式实际上有无限多个根。
是否有任何想法可以将这个方程作为 t 的函数？我想过可能会在表达式的左侧找到函数的反函数，然后取 -t 的反函数，但我不知道从哪里开始。但是如果我能得到有效表达式作为 t 的函数，那么我应该可以毫无问题地对函数进行参数化并在 matplotlib 中绘制它。
编辑
我使用 Runge-Kutta 中点法来逼近 v(t)，这就是图的样子：
这是上面幂级数方程的形状，但我注意到如果我在 desmos 中进一步扩展级数，函数会变得更清晰，更接近于我的近似图的形状。
解决方法
您尝试求解的方程是一个具有浮点系数的多项式方程，并且其系数范围超过 150 个数量级，条件极其恶劣。有了精确的系数，sympy 就可以使用例如计算根Poly(f_ana).nroots() 虽然您可能需要使用非常高的精度，并且由于条件差，计算速度会很慢。
您说这只能通过幂级数方法解决，但实际上替换 x = v' 给出了可分离的 ODE，SymPy 可以解决这个问题：
In [1]: v = Function('v')
In [2]: m,c_d,rho,A,g,t = symbols('m,t')
In [3]: eq = Eq(m*Derivative(v(t),t),c_d*rho*A*v(t)**2/2 - m*g)
In [4]: eq
Out[4]:
2
d A⋅c_d⋅ρ⋅v (t)
m⋅──(v(t)) = ───────────── - g⋅m
dt 2
In [5]: dsolve(eq,v(t))
Out[5]:
_____________
╱ 1
√2⋅g⋅m⋅ ╱ ───────────
╲╱ A⋅c_d⋅g⋅m⋅ρ
v(t) = ─────────────────────────
⎛ √2⋅t ⎞
⎜ C₁ - ──── ⎟
⎜ 2 ⎟
tanh⎜───────────────────⎟
⎜ _____________⎟
⎜ ╱ 1 ⎟
⎜m⋅ ╱ ─────────── ⎟
⎝ ╲╱ A⋅c_d⋅g⋅m⋅ρ ⎠
编辑：有许多不同的方法可以编写此 ODE 的解决方案。由于您现在已经简化了 ODE 本身，我将展示如何从 sympy 获得更简单的答案：
In [86]: a,g = symbols('a,g',positive=True)
In [87]: v = Function('v')
In [88]: eq = Eq(v(t).diff(t),a*v(t)**2 - g)
In [89]: sol = dsolve(eq,v(t))
In [90]: sol
Out[90]:
√g
v(t) = ───────────────────────
√a⋅tanh(√a⋅√g⋅(C₁ - t))
这是 ODE 的解决方案，因为您可以通过将其代回 ODE 来轻松验证。问题是这与您提出的形式有何不同，答案是有多种方法可以编写相同的函数，而且您还需要记住积分常数通常可能很复杂。如果我们采用初始条件 v(0) = 0，那么我们可以求解积分常数：
In [91]: C1 = Symbol('C1')
In [92]: e1,e2 = solve(sol.rhs.subs(t,0),C1)
In [93]: e1
Out[93]:
-ⅈ⋅π
───────
2⋅√a⋅√g
In [94]: sol.subs(C1,e1)
Out[94]:
√g
v(t) = ─────────────────────────────
⎛ ⎛ ⅈ⋅π ⎞⎞
√a⋅tanh⎜√a⋅√g⋅⎜-t - ───────⎟⎟
⎝ ⎝ 2⋅√a⋅√g⎠⎠
In [95]: sol.subs(C1,e1).simplify()
Out[95]:
-√g⋅tanh(√a⋅√g⋅t)
v(t) = ──────────────────
√a
您找到的形式需要不同的积分常数值，但一旦找到特定的解决方案，结果也会相同：
In [97]: sol = Eq(v(t),2*sqrt(g/a) / (C1*exp(2*t*a*sqrt(g/a)) + 1) - sqrt(g/a))
In [98]: sol
Out[98]:
√g 2⋅√g
v(t) = - ── + ──────────────────────
√a ⎛ 2⋅√a⋅√g⋅t ⎞
√a⋅⎝C₁⋅ℯ + 1⎠
In [99]: solve(sol.rhs.subs(t,C1)
Out[99]: [1]
In [100]: e,= solve(sol.rhs.subs(t,C1)
In [101]: sol.subs(C1,e)
Out[101]:
√g 2⋅√g
v(t) = - ── + ───────────────────
√a ⎛ 2⋅√a⋅√g⋅t ⎞
√a⋅⎝ℯ + 1⎠
In [102]: sol.subs(C1,e).rewrite(tanh)
Out[102]:
√g 2⋅√g
v(t) = - ── + ──────────────────────────
√a ⎛ tanh(√a⋅√g⋅t) + 1⎞
√a⋅⎜1 + ─────────────────⎟
⎝ 1 - tanh(√a⋅√g⋅t)⎠
In [103]: sol.subs(C1,e).rewrite(tanh).simplify()
Out[103]:
-√g⋅tanh(√a⋅√g⋅t)
v(t) = ──────────────────
√a
我的第一个答案中的方法在技术上也是正确的解决方案，但还有另一种方法可以得到这个微分方程的解决方案。它是可分的，可以通过部分分数实现积分来解决。最终的解决方案是这样的：
可以在此处找到有关此解决方案如何产生的更多详细信息： https://github.com/gabemorris12/Miscellaneous/blob/main/Engineering%20Analysis/Parachute%20Project%20Merged.pdf
版权声明：本文内容由互联网用户自发贡献，该文观点与技术仅代表作者本人。本站仅提供信息存储空间服务，不拥有所有权，不承担相关法律责任。如发现本站有涉嫌侵权/违法违规的内容， 请发送邮件至 dio@foxmail.com 举报，一经查实，本站将立刻删除。