Sympy second order ode

Question

Ho una soluzione omogenea a un semplice ODE di secondo ordine, che quando provo a risolvere i valori iniziali usando Sympy, restituisce la stessa soluzione. Dovrebbe sostituire y (0) e y '(0) e fornire una soluzione senza costanti, ma non lo fa. Ecco il codice per impostare l'equazione (è un'equazione di bilancio di primavera, k = costante di primavera e m = massa). Alcuni simboli ridondanti che uso altrove, mi dispiace.

%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 

Il risultato è (correttamente): $ 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)), che ha anche y_n = c1.cos (√ (-k/m) t) + c2.sin (√ (-k/m) t).

Quando questa equazione viene risolta analiticamente e si converte in una soluzione utilizzando seni e coseni con omega = sqrt (-k/m), quindi c1 = y (0) e c2 = y '(0)/omega. Quindi, mentre la soluzione sta parzialmente coinvolgendo numeri complessi, naturalmente, dsolve restituisce semplicemente l'equazione omogenea originale come sopra. Il mio codice per valutare l'ODE a Y (0) e y '(0) è:

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

mi rendo conto che dsolve semplicemente non può essere in grado di gestire questa IVP analiticamente, ma sarei sorpreso se questo è così in base alla sua altra capacità. Sarà molto apprezzato qualsiasi aiuto su come questo problema e quindi altri problemi analitici del secondo ordine possano essere risolti. Il nocciolo della questione è tutto il:

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

Quindi la soluzione che ho provato, che Dietrich conferma, era:

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

Mentre è "una" soluzione, questo sembra un modo molto contorto di trovare y (t) = y (0) cos (omega * t - phi) ... che risponde alla domanda implicita su alcune limitazioni di questo risolutore e la domanda diretta su come l'argomento ics viene risolto. Grazie.

Answer 1

Il parametro ics in dsolve() in realtà non funziona (Issue 4720), quindi bisogna fare le sostituzioni manualmente. Si potrebbe provare:

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