Sympy drugiego rzędu od

Question

Mam homogeniczne rozwiązanie prostego rozwiązania ODE o drugiej kolejności, które, gdy próbuję rozwiązać początkowe wartości za pomocą Sympy, zwraca to samo rozwiązanie. Powinien zastąpić y (0) i y '(0) i dać rozwiązanie bez stałych, ale nie. Oto kod do ustawienia równania (jest równanie równowagi wiosennej, k = stała sprężyny i m = masa). Niektóre zbędne symbole, których używam gdzie indziej, przepraszam.

%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 

W rezultacie (prawidłowo) $ y {\ lewo (T \ prawej)} = 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)), który również ma y_n = c1.cos (√ (-k/m) t) + c2.sin (√ (-k/m) t).

Kiedy to równanie zostanie rozwiązane analitycznie i przekształcone w rozwiązanie za pomocą sinusów i cosinusów z omega = sqrt (-k/m), wówczas c1 = y (0) i c2 = y '(0)/omega. Tak więc, chociaż rozwiązanie częściowo dotyczy liczb zespolonych, oczywiście, funkcja dsolve po prostu zwraca oryginalne równanie homogeniczne jak powyżej. Mój kod ocenić ODE w y (0) i y '(0) wynosi:

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

Doceniam to dsolve po prostu może nie być w stanie obsłużyć tego IVP analitycznie, ale byłbym zaskoczony, gdyby tak jest w oparciu o jego inną pojemność. Pomocne będzie rozwiązanie tego problemu, a więc i innych analitycznych problemów drugiego rzędu. Sedno tej kwestii jest wokół:

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

więc rozwiązanie próbowałem, który Dietrich potwierdza, było:

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

Choć jest „a” rozwiązanie to wydaje się być bardzo zawiłe drogi znalezienie y (t) = y (0) cos (omega * t - phi) ... które odpowiada na ukryte pytanie o pewne ograniczenia tego solwera i bezpośrednie pytanie o to, jak rozwiązuje się ics arg. Dzięki.

Answer 1

Parameter ics w dsolve() naprawdę nie działa (Issue 4720), więc trzeba wykonać podstawienia ręcznie. Możesz spróbować:

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