% ### EXhoResonance.m ### 2017.02.04 CB % Code to solve the damped (sinusoidally-) driven harmonic oscillator (DDHO) % for a variety of driving freqs. so to buildup the "resonance curve" via % computation of the mag/phase of the Fourier transform of the steady-state % response. Furthermore, the analytic solution for the DDHO as well as the % transfer function are shown to be equivalent (Fig.1) % Damped driven Harmonic Oscillator (DDHO) % d^2xdt^2 = -((P.wo)^2)*x - P.gamma*dx/dt + (P.A)*sin(P.w*t) % Reqs: % EXhoResonanceFunc.m (re ode45), rfft.m % Notes % o To solve this numerically, need to turn 2nd order ODE into series of 1st order ODEs: % dx/dt= y % dy/dt= -P.wo^2*x - P.gamma*y + (P.A)*sin(P.w*t) % o For autonomous case (i.e., no drive), can rewrite in matrix form such that % A= [0 1;-P.wo^2 -P.gamma]; straight-forward to find associated eigenvalues (see below) % o via P.solveType, user can solve either via ode45 or a hard-coded RK4 % (both should yield the same solution!); Note that (surprisingly) ode45 actually seems % slower than the RK4 (the slowest is ode45 w/ the fixed step-size), possibly due to % the passing of the large-ish structure P; also note that the default ode45 routine % (i.e., adaptive step-size) introduces harmonic disortions in the spectra % due to its nonlinear nature % o For the analytic solution (below manifest as magT and phaseT), the % expression used below, as derived in French (1971) for the steady-state, % is exactly the same as if one simply put in the Fourier transform and % solved for the resulting magnitude and phase [confirmed on the back of an % envelope; let x(t)= X(w)exp(i*w*t) and plug back in, solving for X(w); note then % that magT=abs(X) and phaseT= angle(X)] % o There are a few minor kludges below [e.g., vertical adjustment of the % analytic solution so to match the (arbitrary?) ref. phase of the numeric % solution) % o Impedance (Z) for DDHO is (by definition) the complex ratio of the driving % force and the (steady-state) velocity (see 4080W2016L10REF.pdf). Real part of Z (resistance) % describes energy loss while imaginary part (reactance) describes energy storage % o Comparison of the mags. for the transfer function and admittance % (Fig.5, top right) are a bit kludgy (unsure why fliplr was needed) and % off (worser overlap as you move away from wo) clear; % ================================================================ % --- % Oscillator params. and ICs P.p0 = 0.0; % Initial position {0} P.v0 = 0.0; % Initial velocity {0} P.wo= 10; % resonant (angular) freq {10} P.gamma= 0.5; % damping coefficient {0.5} % --- % Sinusoidal driving term params. P.A = 10; % Driving force amplitude {10} P.wDrive= [5 15]; % start and end angular drive freqs. {[5 15]} P.wDriveN= 25; % # of drive freqs. to run {25} % --- P.tmax = 200; % Maximum time to solve [s; arb] {200} P.SR= 150; % sample rate for time step [Hz; arb] {150} P.Npoints= 8192; % Number of points in time series for FFT, must be 2^n {8192} % --- P.plotN= 1; % boolean re plotting the waveform and spectra for one driving freq. {1} P.plotNnum= round(P.wDriveN/2); % driving freq. index to plot {round(P.wDriveN/2)} % --- P.solveType= 1; % 0-ode45, 1-hard-coded RK4 {1} P.stepF= 0; % boolean re using a fixed step-size for ode45 {0} % ================================================================ % --- dt= 1/P.SR; % spacing of time steps init0 = [P.p0 P.v0]'; % Column vector of initial conditions. tspan = [0:dt:P.tmax]; % time interval for entire computation tW=[0:1/P.SR:(P.Npoints-1)/P.SR]; % (shorter/later) time interval for FFT window L = length(tspan); TW = L-(P.Npoints-1); % create offset point extracting FFT window % --- % create relevant freq. arrays (e.g., for FFT bin labeling) freq= [0:P.Npoints/2]; % Note: these values are not angular (i.e., [freq]= 1/s, not rads/s) freq= P.SR*freq./P.Npoints; df = P.SR/P.Npoints; % freq. spacing between bins wDT= linspace(P.wDrive(1),P.wDrive(2),500); % create ang. freq. array for plotting analytic solution % ++++++++++++++ % various relevant derived quantities (used post- main for loop) Q= P.wo/P.gamma; % "quality factor" (Note: tau=1/P.gamma=Q/P.wo, where tau is time const. of build-up) lambdaP= 0.5*(-P.gamma+ sqrt(P.gamma^2-4*P.wo^2)); % Eigenvalues, for x=0 (undriven) lambdaM= 0.5*(-P.gamma- sqrt(P.gamma^2-4*P.wo^2)); % Note - Can also get eigenvalues via command: eig([0 1;-P.wo^2 -P.gamma]) Z= P.gamma+ i*(wDT- P.wo^2./wDT); % impedance (see notes above; assumes mass is unity) Y= 1./Z; % admittance (reciprocal of impedance) % --- % grabbing driving freqs. from freq array %%indx= find(freq>=P.fDrive(1) & freq<=P.fDrive(2)); % find relevant indicies indx= find(freq>=P.wDrive(1)/(2*pi) & freq<=P.wDrive(2)/(2*pi)); % find relevant indicies indxB= round(linspace(indx(1),indx(end),P.wDriveN)); % one means to get the desired subset freqD= 2*pi*freq(indxB); % array of driving angular freqs % --- for mm=1:numel(freqD) % --- P.w= freqD(mm); % extract driving freq. %P.w= ceil(freqD(mm)/df)*df; % quantized natural freq. % --- % *** Solve in one of two ways *** if P.solveType==0 % use Matlab's ode45 % --- % tell it to actually use the specified step-size if(P.stepF==1), options = odeset('MaxStep',1/P.SR); else options=[]; end [t,y] = ode45(@EXhoResonanceFunc,tspan,init0,options,P); else % use 4th order Runge-Kutta code xPoints(1) = P.p0; vPoints(1) = P.v0; % initialize ICs into dummy arrays x= P.p0; v= P.v0; % kludge dt= 1/P.SR; % time step for nn=1:(length(tspan)-1) % --- t = tspan(nn); % Current time. % --- % step1 xk1= v; vk1= -((P.wo)^2)*x - P.gamma*v + (P.A)*sin(P.w*t); % --- % step 2 xk2 = v + (dt/2)*vk1; vk2= -((P.wo)^2)*(x + (dt/2)*xk1) - P.gamma*(v + (dt/2)*vk1)... + (P.A)*sin(P.w*(t+(dt/2))); % --- % step 3 xk3 = v + (dt/2)*vk2; vk3= -((P.wo)^2)*(x + (dt/2)*xk2) - P.gamma*(v + (dt/2)*vk2)... + (P.A)*sin(P.w*(t+(dt/2))); % --- % step 4 xk4 = v + dt*vk3; vk4= -((P.wo)^2)*(x + (dt)*xk3) - P.gamma*(v + dt*vk3)... + (P.A)*sin(P.w*(t+(dt/2))); % --- % apply RK4 weighting x = x + (dt/6)*(xk1 + 2*xk2 + 2*xk3 + xk4); v = v + (dt/6)*(vk1 + 2*vk2 + 2*vk3 + vk4); % --- % store away position and velocity xPoints(nn+1) = x; vPoints(nn+1) = v; end y(:,1)= xPoints'; y(:,2)= vPoints'; % repackage output end % --- ySPEC= y(TW:TW+P.Npoints-1,1); % steady-state portion of waveform for FFT sigSPEC= rfft(ySPEC); % --- wDrive(mm)= 2*pi*freq(indxB(mm)); % store away driving freqs. mag(mm)= abs(sigSPEC(indxB(mm))); % store away SS mag. % --- % need to correct the phase re the duration of the window allowed for settling into steady-state tPhase= angle(sigSPEC(indxB(mm))); % extract the phase tPhase= angle(exp(i*(tPhase- wDrive(mm)*tspan(TW)))); % correct phase re onset phase(mm)= tPhase; % --- % visualize relevant bits for one of the drive freqs. if mm==P.plotNnum % --- % integrated waveform and segment extracted for spectral analysis figure(2); clf; h1= plot(tspan,y(:,1)); hold on; grid on; xlabel('Time'); ylabel('Position'); title('Time Waveform of integrated solution to damped driven HO equation') L = length(tspan); TW = L-(P.Npoints-1); % create offset point ySPEC= y(TW:TW+P.Npoints-1,1); % steady-state portion of waveform for FFT h2= plot(tspan(TW:TW+P.Npoints-1),ySPEC,'r.','MarkerSize',3); legend([h1 h2],'Entire waveform','Steady-state portion (used for FFT)') % --- % phase space for waveform (entire and steady-state) figure(3); clf; hPS1= plot(y(:,1),y(:,2)); hold on; grid on; hPS2= plot(y(TW:TW+P.Npoints-1,1),y(TW:TW+P.Npoints-1,2),'r.-'); xlabel('Position'); ylabel('Velocity'); title('Phase plane'); legend([hPS1 hPS2],'Entire waveform','Steady-state portion (used for FFT)') % --- % plot spectra of steady-state waveform figure(4); clf; hS1= plot(2*pi*freq,db(sigSPEC)); hold on; grid on; xlabel('Freq [rads/s]'); ylabel('Spectral amplitude [dB]'); hS2= plot(2*pi*freq(indxB(mm)),db(mag(mm)),'rs'); % indicate extracted freq. legend([hS1 hS2],'Steady-state spectra','Driving freq.'); end % --- disp([num2str(100*mm/numel(freqD)),'% done']); end % ++++++++++++++ % [Fig.1] ** Mags/phases extracted from the numeric steady-state responses ** figure(1); clf; subplot(211); hh1= plot(wDrive/P.wo,mag,'ko','MarkerSize',6,'LineWidth',2); hold on; grid on; ylabel('Magnitude'); subplot(212); hh2= plot(wDrive/P.wo,unwrap(phase)/(2*pi),'ko','MarkerSize',6,'LineWidth',2); hold on; grid on; xlabel('Normalized (angular) angular freq (w/wo)'); ylabel('(unwrapped) Phase [cycs]'); % ++++++++++++++ % [Fig.1]** Analytic solution ** (see French, 1971; as noted above, these expressions are % equivalent to using Fourier transforms, which implicitly assume sinusoidal steady-state, to % solving the main ODE) magT= P.A./sqrt((P.wo^2-wDT.^2).^2 + ((P.gamma*wDT).^2)); % mag (theory) phaseT= atan((P.gamma*wDT)./(-P.wo^2+wDT.^2)); % phase (theory; note sign change in denom. re convention) phaseT= phaseT+ phase(1)+ abs(phaseT(1)); % (kludge) correct for (arb?) phase offset in numeric solution figure(1); subplot(211); hh3= plot(wDT/P.wo,magT,'r-','LineWidth',2); subplot(212); hh4= plot(wDT/P.wo,unwrap(2*phaseT)/(4*pi),'r-','LineWidth',2); % kludge to get unwrapping working % ++++++++++++++ % [Fig.1] ** "Transfer function" ** re linear systems theory (i.e., the Fourier transform of the % impulse response of the DHO) init0 = [0 10]'; % set ICs such that there is an "impulse" at t=0 P.w= 0; % make sure to "turn off" drive options= []; [t,yI] = ode45(@EXhoResonanceFunc,tspan,init0,options,P); specI= rfft(yI(1:P.Npoints)); magI= abs(specI); magI= magI* (max(mag)/max(magI)); % scale impulse mag. re max. value of driven case phaseI= angle(specI); phaseI= phaseI+ phase(1); % (kludge) correct for (arb?) phase offset in numeric solution figure(1); subplot(211); hh5= plot(2*pi*freq/P.wo,magI,'b--','LineWidth',2); xlim([wDT(1) wDT(end)]/P.wo); subplot(212); hh6= plot(2*pi*freq/P.wo,unwrap(2*phaseI)/(4*pi),'b--','LineWidth',2); % kludge to get unwrapping working xlim([wDT(1) wDT(end)]/P.wo); % ++++++++++++++ % [Fig.1] Make a legend to put it all together (re Fig.1) figure(1); subplot(211); legend([hh1 hh3 hh5],'Numeric solution re steady-state FFT',... 'Analytic solution','Transfer function'); % ++++++++++++++ % [Fig.5] Plot the impulse response and comparison to admittance figure(5); clf; subplot(221); plot(tW,yI(1:P.Npoints)); hold on; grid on; xlabel('Time [s]'); ylabel('x'); title('Impulse response (no drive; P.w=0, P.p0=0, P.v0=10'); xlim([0 tW(round(numel(tW)/3))]); subplot(222); hI2= plot(freq,db(specI),'LineWidth',2); grid on; hold on; ylabel('Amplitude [dB]'); title('Transfer function (mag. of FFT of IR)'); xlim(P.wDrive/(2*pi)); subplot(224); hI3= plot(freq,angle(specI)/(2*pi),'LineWidth',2); grid on; hold on; xlabel('Frequency [Hz]'); ylabel('Phase [cycles]'); title('Transfer function (phase of FFT of IR)'); xlim(P.wDrive/(2*pi)); subplot(223); hZa= plot(wDT,abs(Y),'k-'); grid on; hold on; hZb= plot(wDT,abs(Z),'r.'); grid on; hold on; ylabel('Amplitude'); xlabel('Ang. requency [rad/s]'); legend([hZa hZb],'admittance','impedance'); % --- % for reference, also include (scaled) admittance to indicate (near?) equivalence offset= max(db(Y))- max(db(specI)); % scaling (in dB) to match up %magY= fliplr(db(Y)- offset); % kludge magY= (db(Y)- offset); subplot(222); hI2b= plot(wDT/(2*pi),magY,'r--'); legend([hI2 hI2b],'Transf. func.','(scaled) Admittance','Location','SouthWest'); angleY= angle(Y)/(2*pi) - angle(Y(1))/(2*pi); % there will be a slight vert. ofseet re angle(specI)/(2*pi) subplot(224); hI3b= plot(wDT/(2*pi),angleY,'r--'); % % ++++++++++++++ % % [Fig.6] Plot energy dynamics for steady-state portion of final driven run % % (assumes mass=1) % figure(6); clf; % U= 0.5* P.wo^2* y(TW:TW+P.Npoints-1,1).^2; % potential energy % K= 0.5* y(TW:TW+P.Npoints-1,2).^2; % potential energy % plot3(tW,U,K,'ko-'); hold on; grid on; % xlabel('Time [s]'); ylabel('Potential energy [J arb.]'); zlabel('Kinetic energy [J arb.]'); % ++++++++++++++ % display some relevant #s to screen disp(['Quality factor (P.wo/P.gamma)= ',num2str(Q)]); disp(['Eigenvalues (for x=0, undriven case): ',num2str(lambdaP),' and ',num2str(lambdaM)]);