Solving Ode Numerically With Python

I am solving an ODE for an harmonic oscillator numerically with Python. When I add a driving force it makes no difference, so I'm guessing something is wrong with the code. Can any

Solution 1:

The problem here is with the term np.cos(Wd*i). It should be np.cos(Wd*i*dT), that is note that dT has been added into the correct equation, since t = i*dT.

If this correction is made, the simulation looks reasonable. Here's a version with F0=0.001. Note that the driving force is clear in the continued oscillations in the damped condition.

The problem with the original equation is that np.cos(Wd*i) just jumps randomly around the circle, rather than smoothly moving around the circle, causing no net effect in the end. This can be best seen by plotting it directly, but the easiest thing to do is run the original form with F0 very large. Below is F0 = 10 (ie, 10000x the value used in the correct equation), but using the incorrect form of the equation, and it's clear that the driving force here just adds noise as it randomly moves around the circle.

Solution 2:

Note that your ODE is well behaved and has an analytical solution. So you could utilize sympy for an alternate approach:

import sympy as sy    
sy.init_printing()  # Pretty printer for IPython

t,k,m,b,F0,Wd = sy.symbols('t,k,m,b,F0,Wd', real=True)  # constants

consts = {k:  0.1, # values
          m:  0.01,
          b:  0.0,
          F0: 0.01,
          Wd: 7.0}

x = sy.Function('x')(t)  # declare variables
dx = sy.Derivative(x, t)
d2x = sy.Derivative(x, t, 2)

# the ODE:
ode1 = sy.Eq(m*d2x + b*dx + k*x, F0*sy.cos(Wd*t))

sl1 = sy.dsolve(ode1, x) # solve ODE
xs1 = sy.simplify(sl1.subs(consts)).rhs # substitute constants# Examining the solution, we note C3 and C4 are superfluous
xs2 = xs1.subs({'C3':0, 'C4':0})
dxs2 = xs2.diff(t)

print("Solution x(t) = ")
print(xs2)
print("Solution x'(t) = ")
print(dxs2)

gives

Solution x(t) = 
C1*sin(3.16227766016838*t) + C2*cos(3.16227766016838*t) - 0.0256410256410256*cos(7.0*t)
Solution x'(t) = 
3.16227766016838*C1*cos(3.16227766016838*t) - 3.16227766016838*C2*sin(3.16227766016838*t) + 0.179487179487179*sin(7.0*t)

The constants C1,C2 can be determined by evaluating x(0),x'(0) for the initial conditions.