% 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-');