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.

Contents

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

Published with MATLAB® R2023a