Problem while using SymPY's dsolve to solve differential equation of criticaly damped oscilator

Question

I have been trying to write a program which will take a set of values and initial conditions for harmonic oscillator, solve the differential equation and display the function of the phi(t).

import sympy as sp
from matplotlib import pyplot as plt
import numpy as np

t = sp.Symbol("t")
phi = sp.Function('phi')

#Values of the oscilator
beta = 0.4
omega_0 = 0.4
f0 = 0
omega = 0

#Initial conditions
phi_0 = 3.14/4
v_phi_0 = 0

#phi''(t) + 2*beta*phi'(t) + omega_0^2*phi(t) - f0*sin(omega*t) = 0
eq = sp.diff(phi(t),t,2)+2*beta*sp.diff(phi(t),t)+omega_0**2*phi(t)-f0*sp.sin(omega*t)
print(sp.dsolve(eq))
equation = sp.dsolve(eq,ics={phi(0):phi_0,sp.diff(phi(t),t).subs(t,0):v_phi_0})
print(equation)
wzor = sp.lambdify(t,equation.rhs,"numpy")
xvals = np.arange(0,30,.1)
yvals = wzor(xvals)

# make figure
fig, ax = plt.subplots(1,1)
ax.plot(xvals, yvals)
ax.set_xlabel('t')
ax.set_ylabel('phi(t)')
plt.show()

Everything seems to be working fine, but when the values of beta and omega_0 are equal (so called critical damping) and are lower than 1, the program appears to be solving the equation wrongly.

Eq(phi(t), (C1*sin(3.65002571393103e-9*t) + C2*sin(3.6500259065127e-9*t) + C3*cos(3.65002571393103e-9*t) + C4*cos(3.6500259065127e-9*t))*exp(-0.4*t))

In other cases, including when beta and omega_0 are equal and higher than 1, there are only 2 constants C1 and C2 after solving the differential equation, but in this case there are 4, which shouldn't be happening. Because of this, the program can't continue and draw the plot, because of these 2 additional constants that I can't solve for.

ax.plot(xvals, yvals)

TypeError: can't convert expression to float

Edit: Also, when a drive force is applied, so f0 and omega are different from 0, following error occurs

NotImplementedError: Cannot find 2 solutions to the homogeneous equation necessary to apply undetermined coefficients to 0.16*phi(t) - 5*sin(0.4*t) + 0.8*Derivative(phi(t), t) + Derivative(phi(t), (t, 2)) (number of terms != order)

It's probably connected to the main problem, but I thought it's good to mention.

Edit 2: Other problems begin to appear when beta and omega_0 are equal and bigger than 3, but are also decimal. No problems appears when they are whole numbers.

Answer 1

Symbolic computation and floating point constants do not go well or at all together. One simple fix is to replace all the constants by rational numbers using the sympy Rational data type,

#Values of the oscilator
beta = sp.Rational(2,5)
omega_0 = sp.Rational(2,5)
f0 = 0
omega = 0

#Initial conditions
phi_0 = sp.Rational(314,400)
v_phi_0 = 0

With only these changes, the output is correctly

Eq(phi(t), (C1 + C2*t)*exp(-2*t/5))
Eq(phi(t), (157*t/500 + 157/200)*exp(-2*t/5))

and the plot is produces as