我有一个简单二阶常微分方程的齐次解,尝试使用 Sympy 求解初始值时,返回相同的解。应该将其代入 y(0) 和 y'(0) 并得到没有常数的解,但实际上并不是这样。以下是建立方程的代码(这是弹簧平衡方程,k=弹簧常数,m=质量)。抱歉使用了一些其他地方也用到的冗余符号。
%matplotlib inline
from sympy import *
m,k, x,t,s,T, omega,A,B = symbols('m k x t s T omega A B',float=True)
a = symbols('a', positive=True)
f,y,g,H,delta,S=symbols('f y g H delta S',cls=Function)
Eq1 = Eq(m*diff(y(t),t,2)+k*y(t))
Eq1
结果是(正确的): $y{\left (t \right )} = C_{1} e^{- t \sqrt{- \frac{k}{m}}} + C_{2} e^{t \sqrt{- \frac{k}{m}}}$
其中 y(t)=C1e^(−t√(−k/m))+C2e^(t√(−km)),也可以写成y_n = c1.cos(√(−k/m)t)+c2.sin(√(−k/m)t)。
当此方程被解析求解,并转换为使用正弦和余弦函数的解时,其中 omega = sqrt(-k/m),那么 c1 = y(0),c2 = y'(0)/omega。因此,尽管该解的一部分涉及复数,但 dsolve 仅返回原始齐次方程如上所述。我计算 ODE 在 y(0) 和 y'(0) 处的代码为:
Eq1_soln_IVP =dsolve(Eq1,y(t),x0=0, ics={y(0): a, y(t).diff(t).subs(t, 0): a})
我理解dsolve可能无法对此IVP进行解析处理,但基于它的其他能力,如果它真的做不到,我会感到惊讶。如果有任何帮助来解决这个问题以及其他解析二阶问题,将不胜感激。问题的核心在于:
ics={y(0): a, y(t).diff(t).subs(t, 0): a}
所以我尝试的解决方案,Dietrich 也确认了,是:
#Create IVP for y(0)
expr = Eq(Eq1_soln_IVP.rhs.subs(sqrt(-k/m),I*omega),y(0))
#Create IVP for y'(0)
expr2 = Eq(diff(y(t),t).subs(t,0),expr.lhs.diff(t))
#Maps all free variables and solves for each where t = 0.
solve([expr.subs(t,0),expr2.subs(t,0)])
虽然这是一种解决方案,但似乎这是一种非常复杂的找到y(t) = y(0)cos(omega*t - phi)的方法......这回答了关于此求解器某些限制的隐含问题以及有关如何解决ics arg的直接问题。谢谢。
m
和k
的符号。如果你写成m, k = sy.symbols('m k', positive=True)
,你将得到实际(物理)解。有相当多的应用程序使用复杂解。顺便说一句,如果你使用Mathematica或Maple,你也将面临同样的问题。 - Dietrich