Sympy二阶ode [英] Sympy second order ode

问题描述

我有一个简单的二阶ODE的同类解决方案，当我尝试使用Sympy求解初始值时，它返回相同的解决方案.它应该替代y(0)和y'(0)并产生没有常数的解，但不是.这是建立方程式的代码(这是一个弹簧平衡方程式，k =弹簧常数，m =质量).很抱歉，我在其他地方使用了一些多余的符号.

I have a homogeneous solution to a simple second-order ODE, which when I try to solve for initial values using Sympy, returns the same solution. It should substitute for y(0) and y'(0) and yield a solution without constants, but does not. Here is the code to set up the equation (it is a spring balance equation, k = spring constant and m = mass). Some redundant symbols which I use elsewhere, sorry.

%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}}} $

The result is (correctly): $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).

y(t)=C1e^(−t√(−k/m))+C2e^(t√(−km)), which also has y_n = c1.cos(√(−k/m)t)+c2.sin(√(−k/m)t).

当对该方程进行解析求解时，使用正弦和余弦以omega = sqrt(-k/m)转换为解，则c1 = y(0)和c2 = y'(0)/omega.因此，尽管解决方案部分地涉及复数，但是dsolve会简单地返回上述原始的齐次方程.我在y(0)和y'(0)处评估ODE的代码是:

When this equation is solved analytically, and converting to a solution using sines and cosines with omega = sqrt(-k/m), then c1 = y(0) and c2 = y'(0)/omega. So, while the solution is partially involving complex numbers, of course, dsolve simply returns the original homogeneous equation as above. My code to evaluate the ODE at y(0) and y'(0) is:

Eq1_soln_IVP =dsolve(Eq1,y(t),x0=0, ics={y(0): a, y(t).diff(t).subs(t, 0): a})

我很欣赏dsolve可能无法解析地处理此IVP，但是如果这基于其其他功能，我将感到惊讶.我们将非常感谢您提供有关如何解决此问题以及由此解决的其他解析二阶问题的任何帮助.问题的核心是:

I appreciate that dsolve simply may not be able to handle this IVP analytically, but I would be surprised if this is so based on its other capacity. Any help as to how this problem and therefore other analytic second order problems can be solved will be much appreciated. The nub of the question is around the :

ics={y(0): a, y(t).diff(t).subs(t, 0): a}

Dietrich确认的我尝试的解决方案是:

So the solution I tried, which Dietrich confirms, was :

 #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)])

虽然这是一个"a"解决方案，但这似乎是一种非常复杂的查找y(t)= y(0)cos(omega * t-phi)的方式...它回答了有关此求解器某些局限性的隐含问题以及有关如何解决ics arg的直接问题.谢谢.

While it is "a" solution, this seems a very convoluted way of finding y(t) = y(0)cos(omega*t - phi)...which answers the implicit question about some limitations of this solver and the direct question about how the ics arg is resolved. Thanks.

推荐答案

dsolve()中的参数ics不能真正起作用(问题4720 )，因此您必须手动进行替换.您可以尝试:

The Parameter ics in dsolve() does not really work (Issue 4720), so you have to do the substitutions manually. You could try:

from IPython.display import display
import sympy as sy

sy.init_printing()  # LaTeX-like pretty printing for IPython

t = sy.Symbol("t", real=True)
m, k = sy.symbols('m k', real=True)  # gives C_1 Exp() + C_2 Exp() solution
# m, k = sy.symbols('m k', positive=True)  # gives C_1 sin() + C_2 cos() sol.
a0, b0 = sy.symbols('a0, b0', real=True)
y = sy.Function('y')

Eq1 = sy.Eq(m*sy.diff(y(t), t, 2) + k*y(t))
print("ODE:")
display(Eq1)

print("Generic solution:")
y_sl0 = sy.dsolve(Eq1, y(t)).rhs  # take only right hand side
display(sy.Eq(y(t), y_sl0))

# Initial conditions:
cnd0 = sy.Eq(y_sl0.subs(t, 0), a0)  # y(0) = a0
cnd1 = sy.Eq(y_sl0.diff(t).subs(t, 0), b0)  # y'(0) = b0

#  Solve for C1, C2:
C1, C2 = sy.symbols("C1, C2")  # generic constants
C1C2_sl = sy.solve([cnd0, cnd1], (C1, C2))

# Substitute back into solution:
y_sl1 = sy.simplify(y_sl0.subs(C1C2_sl))
print("Solution with initial conditions:")
display(sy.Eq(y(t), y_sl1))

这篇关于Sympy二阶ode的文章就介绍到这了，希望我们推荐的答案对大家有所帮助，也希望大家多多支持IT屋！