% RKdamposc.m % Numerical solution to the differential equation of motion for a % damped oscillator, using the Runge-Kutta algorithm % % Raghuveer Parthasarathy % Oct. 24, 2007 % Differential Equation % x'' = - gamma x' - w0^2 x % System of first-order equations % x' = y = f(y) % y' = -gamma y - w0^2 x = g(x,y) clear all close all % Physical Parameters k = 0.1; % Newtons / meter m = 1.0; % kilograms b = 0.03; % Newtons-seconds / meter gamma = b/m; w0 = sqrt(k/m); % Timestep Deltat = (2*pi/w0)/25.0; % 25 pts per period % Total time Tfinal = 10*2*pi/w0; % 10 cycles % Initial conditions x(1) = 6.0; % initial position, meters y(1) = 0.0; % initial velocity, m/s t(1) = 0.0; % initial time, seconds for j=2:round(Tfinal/Deltat) t(j) = Deltat*(j-1); kn1 = y(j-1); ln1 = -gamma*y(j-1) - w0*w0*x(j-1); kn2 = y(j-1) + Deltat*ln1/2; ln2 = -gamma*(y(j-1) + Deltat*ln1/2) - w0*w0*(x(j-1)+Deltat*kn1/2); kn3 = y(j-1) + Deltat*ln2/2; ln3 = -gamma*(y(j-1) + Deltat*ln2/2) - w0*w0*(x(j-1)+Deltat*kn2/2); kn4 = y(j-1) + Deltat*ln3; ln4 = -gamma*(y(j-1) + Deltat*ln3) - w0*w0*(x(j-1)+Deltat*kn3); x(j) = x(j-1) + (Deltat/6.0)*(kn1 + 2*kn2 + 2*kn3 + kn4); y(j) = y(j-1) + (Deltat/6.0)*(ln1 + 2*ln2 + 2*ln3 + ln4); end % Analytic solution -- see Problem Set 4 w = sqrt(w0*w0-gamma*gamma/4); A = x(1)*sqrt(w*w+gamma*gamma/4)/w; phi = atan(gamma/2/w); ta = t(1):Deltat/10:Tfinal; xa = A*exp(-gamma*ta/2).*cos(w*ta-phi); % "Error" calculation % evaluate the analytic solution at the same time points as the numerical % solution xa2 = A*exp(-gamma*t/2).*cos(w*t-phi); chi2 = mean((xa2-x).*(xa2-x)); fs = sprintf('Delta_t = %.3f, chi^2 = %.3e', Deltat, chi2); disp(fs) figure; plot(t,x, 'o-', 'Color', [0 0.2 1.0]) xlabel('time (s)') ylabel('x (m)') title('Damped oscillator'); hold on plot(ta, xa, 'k-');