Sympy second order ode

Question

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

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)), which also has y_n = c1.cos(√(−k/m)t)+c2.sin(√(−k/m)t).

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

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}

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

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.

Answer 1

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