% ### HOode45EX.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] P.y0(2) = 1.0; % initial velocity [m/s] P.b= 0.1; % damping coefficient [kg/s] P.k= 250.0; % stiffness [N m] P.m= 0.01; % mass [kg] % sinusoidal driving term P.A= 0.0; % amplitude [N] (set to zero to turn off) fD= 1.05*sqrt(P.k/P.m)/(2*pi); % freq. (Hz) [expressed as fraction of resonant freq.] % Integration limits P.t0 = 0.0; % Start value P.tf = 3.0; % Finish value P.dt = 0.0001; % time step % ---------------------------------------------------------------------- % +++ % 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]')