Solution to harmonic oscillator with periodic forcing

Question

Using matlab symbolic processor, I can get the homogenous and particular solution to a harmonic oscillator with periodic forcing. I'm trying to write the particular solution in a compact form, with limited success. Does anyone have the general solution to a harmonic oscillator under periodic forcing?

Here's the equation I want to solve

$$ \frac{d^2}{dt^2} h(t) + 2 \zeta \omega_n \frac{d}{dt} h(t) + \omega_n^2 h(t) = \frac{C_i}{m} \sin(\omega_f t + \phi) . $$

Here's one approach

syms h(t) zeta omega_n h0 dh0 Cfi wfi phif m 
Dh(t)=diff(h(t),t);
eqs = diff(h(t), t, t) == -2*zeta*omega_n*Dh-omega_n*omega_n*h(t)+Cfi/m*sin(wfi*t+phif);
disp('equation')
pretty(eqs)
sol=dsolve(eqs, h(0)==h0, Dh(0)==dh0);
pretty(sol)
chk=simplify(diff(sol,t,t)+2*zeta*omega_n*diff(sol,t)+omega_n*omega_n*sol-Cfi/m*sin(wfi*t+phif));
disp('check solution, should be zero');
disp(chk)

Answer 1

If you have $h(t)$ and want to find the forcing function, you can Fourrier analyze it. For each term $a \sin (\psi t)$ in the expansion, the left side becomes $a (\omega^2 -\psi^2)\sin (\psi t) +2a \zeta \omega \cos(\psi t)$ That matches the right side.

Answer 2

The solution to $$ h''(t)+2\zeta\omega_nh'(t)+\omega_n^2h(t)=\frac {C_i} m cos(\Omega_i t + \phi_i) $$ is better written as $$ h''(t)+2\zeta\omega_nh'(t)+\omega_n^2h(t)=\frac {C_i} m cos(\Omega_i t) + \frac {S_i} m sin(\Omega_i t) $$ and solved using

syms zeta omega_n t A B w Cfi Sfi m 
sol = A * cos(w*t) + B * sin(w*t);
eqn1 = diff(sol,t,t) + 2 * zeta * omega_n * diff(sol,t) + omega_n^2 * sol - Sfi/m*sin(w*t)-Cfi/m*cos(w*t);
eqn1 = collect(eqn1,{cos(t*w),sin(t*w)});
expr_1 = children(eqn1);
expr_11 = children(expr_1(1));
expr_12 = children(expr_1(2));
solAB=solve({expr_11(1),expr_12(1)},{A,B});
disp('A');
pretty(solAB.A);
disp('B');
pretty(solAB.B);
disp('A cos + B sin')
pretty(solAB.A*cos(Omega_i*t)+solAB.B*sin(Omega_i*t));
chksoln = (cos(Omega_i*t)*(Cfi*(omega_n^2-Omega_i^2)-2*Sfi*zeta*Omega_i*omega_n) ...
    +sin(Omega_i*t)*(Sfi*(omega_n^2-Omega_i^2)+2*Cfi*zeta*Omega_i*omega_n)) ...
    /(m*((omega_n^2-Omega_i^2)^2+(2*zeta*Omega_i*omega_n)^2));
% chksoln = solAB.A*cos(Omega_i*t)+solAB.B*sin(Omega_i*t)
disp('this should be zero');
chkeqn1 = simplify(diff(chksoln,t,t) + 2 * zeta * omega_n * diff(chksoln,t) + omega_n^2 * chksoln - Sfi/m*sin(Omega_i*t)-Cfi/m*cos(Omega_i*t))

as

$$ h_P(t) = \frac{ \cos \left(\Omega_i t\right) \left( C_i \left( {\omega_n}^2- {\Omega_i}^2 \right) - 2 S_i \zeta \Omega_i \omega_n \right) +\sin \left(\Omega_i t\right) \left( S_i \left( {\omega_n}^2-{\Omega_i}^2 \right) + 2 C_i \zeta \Omega_i \omega_n \right) } {m \left({\left( {\omega_n}^2- {\Omega_i}^2 \right)}^2 + { \left( 2 \Omega_i \omega_n \zeta \right) }^2 \right)} $$