Lets's do the examples we did in the theory lecture.

Consider the simple pendulum, which we showed in the theory lecture is modeled as

If we were to use forward Euler, we would get the following scheme:

omega = 0;

theta = 1;

dt = 0.1;

for i = 1:100

theta = theta + dt*omega;

omega = omega + dt*(-sin(theta));

end

format long

theta

Now lets try it with ode45

omega = 0;

theta = 1;

[T,Y] = ode45(@(t, y) SimplePendulum(t,y),[0 10],[theta omega]);

format long

Y(end,1)

We can use this to plot our phase plane. We first we need a bunch of initial points, and we let them run.

omega = [-2:0.5:2];

theta = [-3*pi:6*pi/5:3*pi];

for i = 1:length(omega)

for j = 1:length(theta)

[T,Y] = ode45(@(t, y) SimplePendulum(t,y),[0 10],[theta(j) omega(i)]);

plot(Y(:, 1), Y(:, 2))

hold on

end

end

hold off

axis([-3*pi 3*pi -3 3])

xlabel('\theta')

ylabel('\omega')

We can even simulate it as a video. Try it for different initial points.

omega = 0;

theta = 1;

[T,Y] = ode45(@(t, y) SimplePendulum(t,y),[0 10],[theta omega]);

y = -cos(Y(:,1));

x = sin(Y(:,1));

for t = 1:length(y)

h = plot([0 x(t)],[0 y(t)],'k',x(t),y(t),'.',0,0,'.',x(1:t),y(1:t));

set(h(1),'linewidth',2);

set(h(2),'MarkerSize',50);

set(h(3),'MarkerSize',20);

axis([-1.5 1.5 -1.5 1.5])

pause(0.1)

end

function dy = SimplePendulum(t,y) %%% This usually goes at the end of the code.

dy = zeros(2,1);

dy(1) = y(2);

dy(2) = -sin(y(1));

end