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

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.

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 