Download
<< Back
OutputCode
%% Question 1(b) - Pendulum
% Comparison of Euler, Runge, and Runge-Kutta for
% $\tau=\frac12,\frac14,\frac18$ to ode23
%% Initial Setup
k=1;
pendulum = @(t,x)[x(2); -k*sin(x(1))];
t0 = 0;
T = 6*pi;
%% $x_0=(-1.5, 0)$
x0 = [-1.5 0];
[t,x]=ode23(pendulum, [t0, T], x0);
fig = figure;
plot(t,x(:,1),'-k','DisplayName','ode23');
hold on;
xlabel('time');
ylabel('x');
[t,x]=eul(pendulum, t0, T, x0, 1/2); plot(t,x(:,1),'-ro','DisplayName','Euler (\tau=1/2)');
[t,x]=eul(pendulum, t0, T, x0, 1/4); plot(t,x(:,1),'-rx','DisplayName','Euler (\tau=1/4)');
[t,x]=eul(pendulum, t0, T, x0, 1/8); plot(t,x(:,1),'-r+','DisplayName','Euler (\tau=1/8)');
[t,x]=runge(pendulum, t0, T, x0, 1/2); plot(t,x(:,1),'-bo','DisplayName','Runge (\tau=1/2)');
[t,x]=runge(pendulum, t0, T, x0, 1/4); plot(t,x(:,1),'-bx','DisplayName','Runge (\tau=1/4)');
[t,x]=runge(pendulum, t0, T, x0, 1/8); plot(t,x(:,1),'-b+','DisplayName','Runge (\tau=1/8)');
[t,x]=rk_classical(pendulum, t0, T, x0, 1/2); plot(t,x(:,1),'-mo','DisplayName','Runge-Kutta (\tau=1/2)');
[t,x]=rk_classical(pendulum, t0, T, x0, 1/4); plot(t,x(:,1),'-mx','DisplayName','Runge-Kutta (\tau=1/4)');
[t,x]=rk_classical(pendulum, t0, T, x0, 1/8); plot(t,x(:,1),'-m+','DisplayName','Runge-Kutta (\tau=1/8)');
legend('Location','NorthEastOutside');
ylim([-5, 5]);
%%
% Due to the oscillations in the solution Euler (especially for low
% $\tau$) quickly diverges from the solution. Runge is better, but towards
% the end it is still slightly divergent; whereas, Runge-Kutta seems good
% for all $\tau$.
%% $x_0=(-3, 0)$
x0 = [-3 0];
[t,x]=ode23(pendulum, [t0, T], x0);
fig = figure;
plot(t,x(:,1),'-k','DisplayName','ode23');
hold on;
xlabel('time');
ylabel('x');
[t,x]=eul(pendulum, t0, T, x0, 1/2); plot(t,x(:,1),'-ro','DisplayName','Euler (\tau=1/2)');
[t,x]=eul(pendulum, t0, T, x0, 1/4); plot(t,x(:,1),'-rx','DisplayName','Euler (\tau=1/4)');
[t,x]=eul(pendulum, t0, T, x0, 1/8); plot(t,x(:,1),'-r+','DisplayName','Euler (\tau=1/8)');
[t,x]=runge(pendulum, t0, T, x0, 1/2); plot(t,x(:,1),'-bo','DisplayName','Runge (\tau=1/2)');
[t,x]=runge(pendulum, t0, T, x0, 1/4); plot(t,x(:,1),'-bx','DisplayName','Runge (\tau=1/4)');
[t,x]=runge(pendulum, t0, T, x0, 1/8); plot(t,x(:,1),'-b+','DisplayName','Runge (\tau=1/8)');
[t,x]=rk_classical(pendulum, t0, T, x0, 1/2); plot(t,x(:,1),'-mo','DisplayName','Runge-Kutta (\tau=1/2)');
[t,x]=rk_classical(pendulum, t0, T, x0, 1/4); plot(t,x(:,1),'-mx','DisplayName','Runge-Kutta (\tau=1/4)');
[t,x]=rk_classical(pendulum, t0, T, x0, 1/8); plot(t,x(:,1),'-m+','DisplayName','Runge-Kutta (\tau=1/8)');
legend('Location','NorthEastOutside');
ylim([-5, 5]);
%%
% Runge is now also diverging from the solution, and even Runge-Kutta
% has some error.
%% $x_0=(-\pi, 1)$
x0 = [-pi 1];
[t,x]=ode23(pendulum, [t0, T], x0);
fig = figure;
plot(t,x(:,1),'-k','DisplayName','ode23');
hold on;
xlabel('time');
ylabel('x');
[t,x]=eul(pendulum, t0, T, x0, 1/2); plot(t,x(:,1),'-ro','DisplayName','Euler (\tau=1/2)');
[t,x]=eul(pendulum, t0, T, x0, 1/4); plot(t,x(:,1),'-rx','DisplayName','Euler (\tau=1/4)');
[t,x]=eul(pendulum, t0, T, x0, 1/8); plot(t,x(:,1),'-r+','DisplayName','Euler (\tau=1/8)');
[t,x]=runge(pendulum, t0, T, x0, 1/2); plot(t,x(:,1),'-bo','DisplayName','Runge (\tau=1/2)');
[t,x]=runge(pendulum, t0, T, x0, 1/4); plot(t,x(:,1),'-bx','DisplayName','Runge (\tau=1/4)');
[t,x]=runge(pendulum, t0, T, x0, 1/8); plot(t,x(:,1),'-b+','DisplayName','Runge (\tau=1/8)');
[t,x]=rk_classical(pendulum, t0, T, x0, 1/2); plot(t,x(:,1),'-mo','DisplayName','Runge-Kutta (\tau=1/2)');
[t,x]=rk_classical(pendulum, t0, T, x0, 1/4); plot(t,x(:,1),'-mx','DisplayName','Runge-Kutta (\tau=1/4)');
[t,x]=rk_classical(pendulum, t0, T, x0, 1/8); plot(t,x(:,1),'-m+','DisplayName','Runge-Kutta (\tau=1/8)');
legend('Location','NorthEastOutside');
ylim([-5, 5]);
%%
% All solutions are fairly good.
