% ### EXharmOSCode45.m ### 09.19.14 % Numerically integrate the damped/driven harmonic oscillator % m*x''+ b*x' + k*x = A*sin(wt) clear % ----------------------------------------------------- % User input (Note: All paramters are stored in a structure) P.y0(1) = 0.0; % initial position [m] {0} P.y0(2) = 1.0; % initial velocity [m/s] {1} P.b= 0.05; % damping coefficient [kg/s] {0.05} P.k= 250.0; % stiffness [N m] {250} P.m= 0.01; % mass [kg] {0.01} % sinusoidal driving term P.A= 1.0; % amplitude [N] (set to zero to turn off) {1} fD= 1.001*sqrt(P.k/P.m)/(2*pi); % freq. (Hz) [expressed as fraction of resonant freq.] % Integration limits P.t0 = 0.0; % Start value {0} P.tf = 3.0; % Finish value {3} P.dt = 0.0001; % time step {0.0001} % ---------------------------------------------------------------------- % +++ % spit back out some basic derived quantities P.wr= 2*pi*fD; % convert to angular freq. disp(sprintf('Resonant frequency ~%g [Hz]', sqrt(P.k/P.m)/(2*pi))); Q = (sqrt(P.k/P.m))/(P.b/P.m); % quality factor disp(sprintf('Q-value = %g', Q)); % +++ % use built-in ode45 to solve [t y] = ode45('HOfunction', [P.t0:P.dt:P.tf],P.y0,[],P); % ------------------------------------------------------ % visualize figure(1); clf; plot(t,y(:,1)); hold on; grid on; xlabel('t [s]'); ylabel('x(t) [m]') % Phase plane figure(2); clf; plot(y(:,1), y(:,2)); hold on; grid on; xlabel('x [m]'); ylabel('dx/dt [m/s]') % ------------------------------------------------------ % compute Hilbert transform to extract inst. mag. and phase if 1==0 wH= hilbert(y(:,1)); figure(3); clf; subplot(211); hold on; grid on; title('Analytic signal representation') plot(t,abs(wH)); xlabel('t [s]'); ylabel('Magnitude') subplot(212); hold on; grid on; title('Analytic signal representation') ang= cycs(wH)- fD*t; % subtract accumulating phase settling into driving freq. (better way to do this?) plot(t,ang); xlabel('t [s]'); ylabel('Phase [cyc]') %theta= angle(wH)-angle(exp(-i*P.wr*t)); %plot(t,(theta)); xlabel('t [s]'); ylabel('Phase [cyc]') end