6

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.

Time Lord
  • 331
  • 2
  • 5
  • 12

1 Answers1

4

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))
Dietrich
  • 5,241
  • 3
  • 24
  • 36
  • Thanks Dietrich for noting the issue 4720 and for the independent confirmation of the only workaround I could construct. Answered with my appreciation. – Time Lord Jan 13 '16 at 09:07
  • So the mathematical issue here is that the solver doesn't utilise substitutions using the Euler equations. Turning the e^{\pm i \omega t} terms in the solution into m.cos(\omega t) + n.i sin(\omega t) form is critical to finding the real and physical meaning of these equations, in this case a spring with a weight suspended on the end where y(t) is an oscillating displacement. Sympy. plot copes with the trig form but not at all with complex forms, precluding effective visualisation, also. – Time Lord Jan 13 '16 at 09:19
  • It depends on the signs of `m` and `k`. If you write `m, k = sy.symbols('m k', positive=True)`, you'll get the real (physical) solution. There are quite a few applications, where complex solutions are used. BTW, you'll have to cope with the same problem, if you use Mathematica or Maple. – Dietrich Jan 13 '16 at 19:14
  • It has been a useful exercise to look at Sympy's capacity. I guess my thought with no particular depth is that solvers would need to have a method of converting complex to trigonometric forms for some solutions. The problem of course in engineering applications is that m and k coefficients typically represent opposing forces, necessitating one of them to be negative, and hence a complex solution. Good to know that Maple and Mathematica are no different. Thanks for the input. – Time Lord Jan 14 '16 at 20:05