%%% 3-inverter ring oscillator dynamics
%%% BENG 221 example

A = [0 0 -1; -1 0 0; 0 -1 0] % LTI ODE matrix for ring oscillator
[X, s] = eig(A) % eigenvalues s, and eigenvectors X

V0 = [1; 0; 0]; % initial conditions
% V0 = [1; 1; 1]; % IC exciting the first, decaying eigenmode
% V0 = [1; -0.5; -0.5]; % IC exciting the other, increasing oscillatory eigenmodes

tstep = 0.1; % time step
trange = 10; % range of time
t = 0:tstep:trange; % vector of time points

% analytic solution
V = X * diag(X' * V0) * exp(diag(s) * t);
figure(1)
plot(t, V, 'LineWidth', 2)
legend('v1', 'v2', 'v3', 'Location', 'Northwest')
xlabel('t')
title(['IC: (', num2str(V0'), ')'])

% Euler integration
hold on % continue plotting on the same graph
Ve = V0; % Euler approximation, initialize to IC
Vs = V0; % Euler state variable, initialize
for te = tstep:tstep:trange
    Vs = Vs + A * Vs * tstep;
    Ve = [Ve, Vs];
end
plot(t, Ve, '+', 'MarkerSize', 3)
hold off

% % Crank-Nicolson integration
% hold on % continue plotting on the same graph
% G = (eye(3) - A * tstep / 2) \ (eye(3) + A * tstep / 2); % CN update matrix
% Vc = V0; % CN approximation, initialize to IC
% Vs = V0; % CN state variable, initialize
% for te = tstep:tstep:trange
%     Vs = G* Vs;
%     Vc = [Vc, Vs];
% end
% plot(t, Vc, 'x', 'MarkerSize', 3)
% hold off
% 
% % forward second-order integration
% hold on % continue plotting on the same graph
% G = (eye(3) + A * tstep + A * A * tstep^2 / 2); % CN update matrix
% Vf = V0; % forward 2nd-order approximation, initialize to IC
% Vs = V0; % forward 2nd-order state variable, initialize
% for te = tstep:tstep:trange
%     Vs = G* Vs;
%     Vf = [Vf, Vs];
% end
% plot(t, Vf, 'o', 'MarkerSize', 3)
% hold off
% 
% % matlab built-in ODE solver
% hold on % continue plotting on the same graph
% OdeOptions = odeset('RelTol', 1e-9); % accuracy, please
% dVdt = @(t,V) A * V; 
% [tm, Vm] = ode23(dVdt, [0 trange], V0);   
% plot(tm, Vm, '*', 'MarkerSize', 3)
% hold off
% 
% print -depsc 'ring_100.eps'
