Download
<< Back
OutputCode
%% Question 1(c) - Harmonic Oscillator
% Comparison of Euler, Runge, and Runge-Kutta for
% $\tau=\frac12,\frac14,\frac18$ to ode23
%
% Note; the simulations are only ran for $t\in[0, 20]$ for ease of demonstration/visualisation.
%% Initial Setup
a = 0;
b = 9;
c = 10;
t0 = 0;
T = 20;
x0 = [1 0];
%% $\omega=2.5$
omega = 2.5;
oscillator = @(t,x)[0 1; -b -a]*x + [0; c*cos(omega*t)];
[t1,x1]=ode23(oscillator, [t0, T], x0);
figure;
plot(t1,x1(:,1),'-k','DisplayName','ode23');
hold on;
title("Euler");
xlabel('time');
ylabel('x');
[t,x]=eul(oscillator, t0, T, x0, 1/2); plot(t,x(:,1),'-ro','DisplayName','Euler (\tau=1/2)');
[t,x]=eul(oscillator, t0, T, x0, 1/4); plot(t,x(:,1),'-bx','DisplayName','Euler (\tau=1/4)');
[t,x]=eul(oscillator, t0, T, x0, 1/8); plot(t,x(:,1),'-m+','DisplayName','Euler (\tau=1/8)');
legend('Location','NorthWest');
ylim([-10, 10]);
% Following increase plot width
pos = get(gcf,'position');
pos(3) = pos(3)*1.5;
set(gcf,'position', pos);
figure;
plot(t1,x1(:,1),'-k','DisplayName','ode23');
hold on;
title("Runge");
xlabel('time');
ylabel('x');
[t,x]=runge(oscillator, t0, T, x0, 1/2); plot(t,x(:,1),'-ro','DisplayName','Runge (\tau=1/2)');
[t,x]=runge(oscillator, t0, T, x0, 1/4); plot(t,x(:,1),'-bx','DisplayName','Runge (\tau=1/4)');
[t,x]=runge(oscillator, t0, T, x0, 1/8); plot(t,x(:,1),'-m+','DisplayName','Runge (\tau=1/8)');
legend('Location','NorthWest');
ylim([-10, 10]);
% Following increase plot width
pos = get(gcf,'position');
pos(3) = pos(3)*1.5;
set(gcf,'position', pos);
figure;
plot(t1,x1(:,1),'-k','DisplayName','ode23');
hold on;
title("Runge-Kutta");
xlabel('time');
ylabel('x');
[t,x]=rk_classical(oscillator, t0, T, x0, 1/2); plot(t,x(:,1),'-ro','DisplayName','Runge-Kutta (\tau=1/2)');
[t,x]=rk_classical(oscillator, t0, T, x0, 1/4); plot(t,x(:,1),'-bx','DisplayName','Runge-Kutta (\tau=1/4)');
[t,x]=rk_classical(oscillator, t0, T, x0, 1/8); plot(t,x(:,1),'-m+','DisplayName','Runge-Kutta (\tau=1/8)');
legend('Location','NorthWest');
ylim([-10, 10]);
% Following increase plot width
pos = get(gcf,'position');
pos(3) = pos(3)*1.5;
set(gcf,'position', pos);
%% $\omega=2.9$
omega = 2.9;
oscillator = @(t,x)[0 1; -b -a]*x + [0; c*cos(omega*t)];
[t1,x1]=ode23(oscillator, [t0, T], x0);
figure;
plot(t1,x1(:,1),'-k','DisplayName','ode23');
hold on;
title("Euler");
xlabel('time');
ylabel('x');
[t,x]=eul(oscillator, t0, T, x0, 1/2); plot(t,x(:,1),'-ro','DisplayName','Euler (\tau=1/2)');
[t,x]=eul(oscillator, t0, T, x0, 1/4); plot(t,x(:,1),'-bx','DisplayName','Euler (\tau=1/4)');
[t,x]=eul(oscillator, t0, T, x0, 1/8); plot(t,x(:,1),'-m+','DisplayName','Euler (\tau=1/8)');
legend('Location','NorthWest');
ylim([-30, 30]);
% Following increase plot width
pos = get(gcf,'position');
pos(3) = pos(3)*1.5;
set(gcf,'position', pos);
figure;
plot(t1,x1(:,1),'-k','DisplayName','ode23');
hold on;
title("Runge");
xlabel('time');
ylabel('x');
[t,x]=runge(oscillator, t0, T, x0, 1/2); plot(t,x(:,1),'-ro','DisplayName','Runge (\tau=1/2)');
[t,x]=runge(oscillator, t0, T, x0, 1/4); plot(t,x(:,1),'-bx','DisplayName','Runge (\tau=1/4)');
[t,x]=runge(oscillator, t0, T, x0, 1/8); plot(t,x(:,1),'-m+','DisplayName','Runge (\tau=1/8)');
legend('Location','NorthWest');
ylim([-30, 30]);
% Following increase plot width
pos = get(gcf,'position');
pos(3) = pos(3)*1.5;
set(gcf,'position', pos);
figure;
plot(t1,x1(:,1),'-k','DisplayName','ode23');
hold on;
title("Runge-Kutta");
xlabel('time');
ylabel('x');
[t,x]=rk_classical(oscillator, t0, T, x0, 1/2); plot(t,x(:,1),'-ro','DisplayName','Runge-Kutta (\tau=1/2)');
[t,x]=rk_classical(oscillator, t0, T, x0, 1/4); plot(t,x(:,1),'-bx','DisplayName','Runge-Kutta (\tau=1/4)');
[t,x]=rk_classical(oscillator, t0, T, x0, 1/8); plot(t,x(:,1),'-m+','DisplayName','Runge-Kutta (\tau=1/8)');
legend('Location','NorthWest');
ylim([-30, 30]);
% Following increase plot width
pos = get(gcf,'position');
pos(3) = pos(3)*1.5;
set(gcf,'position', pos);
%% $\omega=3.1$
omega = 3.1;
oscillator = @(t,x)[0 1; -b -a]*x + [0; c*cos(omega*t)];
[t1,x1]=ode23(oscillator, [t0, T], x0);
figure;
plot(t1,x1(:,1),'-k','DisplayName','ode23');
hold on;
title("Euler");
xlabel('time');
ylabel('x');
[t,x]=eul(oscillator, t0, T, x0, 1/2); plot(t,x(:,1),'-ro','DisplayName','Euler (\tau=1/2)');
[t,x]=eul(oscillator, t0, T, x0, 1/4); plot(t,x(:,1),'-bx','DisplayName','Euler (\tau=1/4)');
[t,x]=eul(oscillator, t0, T, x0, 1/8); plot(t,x(:,1),'-m+','DisplayName','Euler (\tau=1/8)');
legend('Location','NorthWest');
ylim([-30, 30]);
% Following increase plot width
pos = get(gcf,'position');
pos(3) = pos(3)*1.5;
set(gcf,'position', pos);
figure;
plot(t1,x1(:,1),'-k','DisplayName','ode23');
hold on;
title("Runge");
xlabel('time');
ylabel('x');
[t,x]=runge(oscillator, t0, T, x0, 1/2); plot(t,x(:,1),'-ro','DisplayName','Runge (\tau=1/2)');
[t,x]=runge(oscillator, t0, T, x0, 1/4); plot(t,x(:,1),'-bx','DisplayName','Runge (\tau=1/4)');
[t,x]=runge(oscillator, t0, T, x0, 1/8); plot(t,x(:,1),'-m+','DisplayName','Runge (\tau=1/8)');
legend('Location','NorthWest');
ylim([-30, 30]);
% Following increase plot width
pos = get(gcf,'position');
pos(3) = pos(3)*1.5;
set(gcf,'position', pos);
figure;
plot(t1,x1(:,1),'-k','DisplayName','ode23');
hold on;
title("Runge-Kutta");
xlabel('time');
ylabel('x');
[t,x]=rk_classical(oscillator, t0, T, x0, 1/2); plot(t,x(:,1),'-ro','DisplayName','Runge-Kutta (\tau=1/2)');
[t,x]=rk_classical(oscillator, t0, T, x0, 1/4); plot(t,x(:,1),'-bx','DisplayName','Runge-Kutta (\tau=1/4)');
[t,x]=rk_classical(oscillator, t0, T, x0, 1/8); plot(t,x(:,1),'-m+','DisplayName','Runge-Kutta (\tau=1/8)');
legend('Location','NorthWest');
ylim([-30, 30]);
% Following increase plot width
pos = get(gcf,'position');
pos(3) = pos(3)*1.5;
set(gcf,'position', pos);
%% $\omega=3$
omega = 3;
oscillator = @(t,x)[0 1; -b -a]*x + [0; c*cos(omega*t)];
[t1,x1]=ode23(oscillator, [t0, T], x0);
figure;
plot(t1,x1(:,1),'-k','DisplayName','ode23');
hold on;
title("Euler");
xlabel('time');
ylabel('x');
[t,x]=eul(oscillator, t0, T, x0, 1/2); plot(t,x(:,1),'-ro','DisplayName','Euler (\tau=1/2)');
[t,x]=eul(oscillator, t0, T, x0, 1/4); plot(t,x(:,1),'-bx','DisplayName','Euler (\tau=1/4)');
[t,x]=eul(oscillator, t0, T, x0, 1/8); plot(t,x(:,1),'-m+','DisplayName','Euler (\tau=1/8)');
legend('Location','NorthWest');
ylim([-30, 30]);
% Following increase plot width
pos = get(gcf,'position');
pos(3) = pos(3)*1.5;
set(gcf,'position', pos);
figure;
plot(t1,x1(:,1),'-k','DisplayName','ode23');
hold on;
title("Runge");
xlabel('time');
ylabel('x');
[t,x]=runge(oscillator, t0, T, x0, 1/2); plot(t,x(:,1),'-ro','DisplayName','Runge (\tau=1/2)');
[t,x]=runge(oscillator, t0, T, x0, 1/4); plot(t,x(:,1),'-bx','DisplayName','Runge (\tau=1/4)');
[t,x]=runge(oscillator, t0, T, x0, 1/8); plot(t,x(:,1),'-m+','DisplayName','Runge (\tau=1/8)');
legend('Location','NorthWest');
ylim([-30, 30]);
% Following increase plot width
pos = get(gcf,'position');
pos(3) = pos(3)*1.5;
set(gcf,'position', pos);
figure;
plot(t1,x1(:,1),'-k','DisplayName','ode23');
hold on;
title("Runge-Kutta");
xlabel('time');
ylabel('x');
[t,x]=rk_classical(oscillator, t0, T, x0, 1/2); plot(t,x(:,1),'-ro','DisplayName','Runge-Kutta (\tau=1/2)');
[t,x]=rk_classical(oscillator, t0, T, x0, 1/4); plot(t,x(:,1),'-bx','DisplayName','Runge-Kutta (\tau=1/4)');
[t,x]=rk_classical(oscillator, t0, T, x0, 1/8); plot(t,x(:,1),'-m+','DisplayName','Runge-Kutta (\tau=1/8)');
legend('Location','NorthWest');
ylim([-30, 30]);
% Following increase plot width
pos = get(gcf,'position');
pos(3) = pos(3)*1.5;
set(gcf,'position', pos);
%% $\omega=\sqrt{3}$
omega = sqrt(3);
oscillator = @(t,x)[0 1; -b -a]*x + [0; c*cos(omega*t)];
[t1,x1]=ode23(oscillator, [t0, T], x0);
figure;
plot(t1,x1(:,1),'-k','DisplayName','ode23');
hold on;
title("Euler");
xlabel('time');
ylabel('x');
[t,x]=eul(oscillator, t0, T, x0, 1/2); plot(t,x(:,1),'-ro','DisplayName','Euler (\tau=1/2)');
[t,x]=eul(oscillator, t0, T, x0, 1/4); plot(t,x(:,1),'-bx','DisplayName','Euler (\tau=1/4)');
[t,x]=eul(oscillator, t0, T, x0, 1/8); plot(t,x(:,1),'-m+','DisplayName','Euler (\tau=1/8)');
legend('Location','NorthWest');
ylim([-5, 5]);
% Following increase plot width
pos = get(gcf,'position');
pos(3) = pos(3)*1.5;
set(gcf,'position', pos);
figure;
plot(t1,x1(:,1),'-k','DisplayName','ode23');
hold on;
title("Runge");
xlabel('time');
ylabel('x');
[t,x]=runge(oscillator, t0, T, x0, 1/2); plot(t,x(:,1),'-ro','DisplayName','Runge (\tau=1/2)');
[t,x]=runge(oscillator, t0, T, x0, 1/4); plot(t,x(:,1),'-bx','DisplayName','Runge (\tau=1/4)');
[t,x]=runge(oscillator, t0, T, x0, 1/8); plot(t,x(:,1),'-m+','DisplayName','Runge (\tau=1/8)');
legend('Location','NorthWest');
ylim([-5, 5]);
% Following increase plot width
pos = get(gcf,'position');
pos(3) = pos(3)*1.5;
set(gcf,'position', pos);
figure;
plot(t1,x1(:,1),'-k','DisplayName','ode23');
hold on;
title("Runge-Kutta");
xlabel('time');
ylabel('x');
[t,x]=rk_classical(oscillator, t0, T, x0, 1/2); plot(t,x(:,1),'-ro','DisplayName','Runge-Kutta (\tau=1/2)');
[t,x]=rk_classical(oscillator, t0, T, x0, 1/4); plot(t,x(:,1),'-bx','DisplayName','Runge-Kutta (\tau=1/4)');
[t,x]=rk_classical(oscillator, t0, T, x0, 1/8); plot(t,x(:,1),'-m+','DisplayName','Runge-Kutta (\tau=1/8)');
legend('Location','NorthWest');
ylim([-5, 5]);
% Following increase plot width
pos = get(gcf,'position');
pos(3) = pos(3)*1.5;
set(gcf,'position', pos);
