Resonance in non-homogeneous equations

We have used the method of variation of parameters to consider the non-homogeneous equation $y'' + \omega_F^2 y = cos(\omega_I t)$ when $\omega_F \neq \omega_I$ and also when $\omega_F = \omega_I$ . Somewhat remarkably, we get different behavior for the formula of the solutions in the two cases. The general solution to the associated homogeneous equation is the same $c_1 \cos(\omega_F t) + c_2 \sin(\omega_F t), \qquad c_1, c_2 \in \mathbb{R}$ but our particular solution depends on $\omega_I$ . We have $P(t) = \begin{cases} \frac{1}{\omega_F^2 - \omega_I^2} \cos(\omega_I t) & \omega_F \neq \omega_I \\ \frac{t}{2\omega_F} \sin(\omega_I t) & \omega_F = \omega_I\end{cases}.$ This notebook lets us study the solutions to an IVP $y(0) = y_0, y'(0) = y_0'$ as we vary $\omega_I$ to approach $\omega_S$ and see that, despite the formula for the solutions being different, the solutions to the IVP vary nicely as $\omega_I$ varies.

In [3]:

@interact
def plot_2nd_order_inhomogeneous(of = slider(vmin = 0.1, vmax = 5, default = 2, step_size = 0.1, label = "$\\omega_F$"),
                               oi = slider(vmin = 0.1, vmax = 6, default = 1.2, step_size = 0.02, label = '$\\omega_I$'), 
                               y0 = slider(vmin = -3, vmax = 3, default = .2, step_size = 0.05, label = '$y(0) = $'),
                               yp0 = slider(vmin = -3, vmax = 3, default = 0, step_size = 0.05, label = '$y\'(0) = $'),
                               tmax = slider(vmin = 5, vmax = 50, default = 30, step_size = 1, label = "max $t$ value"),
                               non_forced = checkbox(False, label = "Plot sol. to homog. eq."),
                               forcing_term = checkbox(True, label = "Plot forcing term"),
                               auto_update = True):
    t = var('t')
    cv1, cv2 = var('cv1,cv2')
    o_f = of.numerical_approx(digits = 3)
    o_i = oi.numerical_approx(digits = 3)
    P(t) = t/(2*o_f)*sin(o_i*t) if o_f == o_i else 1/(o_f^2 - o_i^2)*cos(o_i*t)
    solve_c = solve([cv1*cos(0) + cv2*sin(0) + P(t = 0) == y0,
                     -cv1*o_f*sin(0) + cv2*o_f*cos(0) + diff(P(t), t)(t = 0) == yp0], 
                   cv1, cv2, solution_dict = True)
    c1,c2 = solve_c[0][cv1], solve_c[0][cv2]
    ivp_sol(t) = c1*cos(o_f*t) + c2*sin(o_f*t) + P(t)
    plt = plot(ivp_sol, (t,0, tmax), 
              title = "Solution plot", axes_labels = ['$t$', '$y$'], legend_label = "$y(t)$")
    if non_forced:
        solve_c_nf = solve([cv1*cos(0) + cv2*sin(0) == y0,
                     -cv1*o_f*sin(0) + cv2*o_f*cos(0) == yp0],
                   cv1, cv2, solution_dict = True)
        c1nf,c2nf = solve_c_nf[0][cv1], solve_c_nf[0][cv2]
        plt += plot(c1nf*cos(o_f*t) + c2nf*sin(o_f*t), (t, 0, tmax), color = 'red', alpha = 0.2, thickness = 1 , legend_label = "homog. sol'n")
    if forcing_term:
        plt += plot(cos(o_i*t), (t, 0, tmax), color = 'orange', alpha = 0.3, thickness =1, legend_label = "forcing term")
    show(plt, ymin = -10, ymax = 10, figsize = 6)

In [0]: