How do i replace Euler's method with ode45?

Hi, I'm trying to replace the Section of the code that uses Euler's methos to solve the ode by ode45,

and i have the individual codes, but i'm unsure how to insert one into another:

First, I have my main code, which uses Euler's method to solve the ode( TBH, i'm not sure if it even solves an ode, because there isn't

an explicitly defined ode in the code, more like a rotation matrix, but it gives me the output of an oscillating pendulum)

(Basically, i do not understand how the motion of the pendulum is plotted using a rotation matrix)

function Animation
% A script to animate the motion of the simple pendulum
clc
clear
close all;
disp('Specify the initial angle of the pendulum in degrees, e.g., 50')
disp('or press ENTER for the default value.')
theta=input('Initial angular displacement of the pendulum= ');
if isempty(theta)
    theta=50
else
    theta;
end
theta=theta*pi/180; % convert from degrees to radians
disp('Specify the final time, e.g. tfinal = 25')
disp('or press ENTER for default value.')
tfinal=input('Final time tfinal = ');
if isempty(tfinal)
    tfinal=10
else
    tfinal;
end
data_init=[0 0; -4 0]; %pendulum length
% rotation matrix
R=[cos(theta) -sin(theta);sin(theta) cos(theta)];
data=R*data_init; % initial pendulum position
bar=line('xdata',[0 data(1)],'ydata',...
[0 data(2)]','linewidth',3);
mass=line('xdata',data(1),'ydata',data(2),'marker','o',...
'markersize',15,'MarkerFaceColor','b');
hinge=line('xdata',0,'ydata',0,'marker','o',...
'markersize',7);
axis([-5 5 -5 5])
grid % comment this out if you do like the grid
set(gca,'Fontsize',14)
set(gca,'dataaspectratio',[1 1 1])
box on
dt=0.03; % step−size for solving differential
% equations can be arbitrarily selected
t=0; % initial time
thetadot=0; % initial angular speed
disp('Can maximize the display so you can see the action better,')
disp('then press ENTER.')
disp('If you are happy with the display size, press ENTER.')
pause
 
v = VideoWriter('pendulum4.avi');
open(v);
% solve the diff eqns using the Euler's method
while(t<tfinal)
    t=t+dt;
    theta = theta + thetadot*dt;
    thetadot=thetadot - 5*sin(theta)*dt; %−0.01*thetadot;
    % you can add some friction
    R=[cos(theta) -sin(theta);sin(theta) cos(theta)];
    datanew=R*data_init;
    % change the property values of the bar and hinge objects
    set(bar,'xdata',[0 datanew(1)],'ydata',[0 datanew(2)]);
    set(mass,'xdata',datanew(1),'ydata',datanew(2));
    set(hinge,'xdata',0,'ydata',0)
    set(gca,'dataaspectratio',[1 1 1])
    drawnow;
    % get the frame and write to the file
    frame = getframe(gcf);
    writeVideo(v,frame);
end
% close the video writer
close(v);
end

Next, I have my code which uses ode45 to solve the equation of motion of the pendulum, which solves both the linear and nonlinear equations of motion using two functions, pend_l ( for linear) and pend_n ( for nonlinear)

clear all;
clc;
clf;
tic;
tspan = 0:0.01:20;
a=pi/8;
b=0;
x0 = [a; b];
[t1,x] = ode45(@pend_l,tspan,x0);
X1 = x(:,1);
X2 = x(:,2);
y0 = [a ; b];
[t2,y] = ode45(@pend_n,tspan,y0);
Y1 = y(:,1);
Y2 = y(:,2);

%figure(1);
subplot(2,2,1);
 plot(t1,X1);    %linear disp vs time
 xlabel('Time (s)');
 ylabel('Displacement (rad)');
 hold on;
 grid on;
 plot(t2,Y1);    %nonlinear disp vs time
 legend('Linear','Non Linear');
 
%figure(2);
subplot(2,2,2);
 plot(t1,X2);    %linear velocity vs time
 xlabel('Time (s)');
 ylabel('Velocity (rad/s)');
 hold on;
 grid on;
 plot(t2,Y2);    %nonlinear velocity vs time
 
 subplot(2,2,3);
 plot(X1,X2);     %linear vel vs disp
 hold on;
 plot(Y1,Y2);     %nonlinear vel vs disp
 xlabel('Displacement (rad)');
 ylabel('Velocity (rad/s)');
 grid on;

toc;

function ydot = pend_n(t2,y)
wsq=3.5;         % same value of wsq in both cases is required
ydot = [y(2); -wsq*sin(y(1))];
end

function xdot = pend_l(t1,x)
wsq=3.5;        
xdot = [x(2); -wsq*x(1)];
end

Answer 1

[Alan Stevens]

Here's one way. I've simplified the data input - you can add the bells and whistes!

% A script to animate the motion of the simple pendulum
theta0=deg2rad(50);
tfinal=10;
thetadot0=0; % initial angular speed
tspan = [0 tfinal];
IC = [theta0 thetadot0];

[t, Y] = ode45(@odefn, tspan,IC);
theta = Y(:,1);
thetadot = Y(:,2);

data_init=[0 0; -4 0]; %pendulum length
R=[cos(theta(1)) -sin(theta(1));sin(theta(1)) cos(theta(1))];
data=R*data_init; % initial pendulum position
bar=line('xdata',[0 data(1)],'ydata',...
[0 data(2)]','linewidth',3);
mass=line('xdata',data(1),'ydata',data(2),'marker','o',...
'markersize',15,'MarkerFaceColor','b');
hinge=line('xdata',0,'ydata',0,'marker','o',...
'markersize',7);
axis([-5 5 -5 5])
grid % comment this out if you do like the grid
set(gca,'Fontsize',14)
set(gca,'dataaspectratio',[1 1 1])
box on

for i=2:numel(t)
    R=[cos(theta(i)) -sin(theta(i));sin(theta(i)) cos(theta(i))];
    datanew=R*data_init;
    set(bar,'xdata',[0 datanew(1)],'ydata',[0 datanew(2)]);
    set(mass,'xdata',datanew(1),'ydata',datanew(2));
    set(hinge,'xdata',0,'ydata',0)
    set(gca,'dataaspectratio',[1 1 1])
    drawnow;
    frame = getframe(gcf);
end

function dYdt = odefn(~,Y)

        theta = Y(1);
        thetadot = Y(2);
        dYdt = [thetadot;
               -5*sin(theta)];
end

[KLETECH MOTORSPORTS]

is that the same as this ? i tried to incorporate my ode45 function and script more directly into the animation code. but it is not oscillating, and it is not getting saved.

function SimplePendAnimation
% A script to animate the motion of the simple pendulum
clc
clear
close all;
disp('Specify the initial angle of the pendulum in degrees, e.g., 50')
disp('or press ENTER for the default value.')
a=input('Initial angular displacement of the pendulum= ');
if isempty(a)
    a=50;
else
    a;
end
a=a*pi/180; % convert from degrees to radians

data_init=[0 0; -4 0]; %pendulum length
% rotation matrix
R=[cos(a) -sin(a);sin(a) cos(a)];
data=R*data_init; % initial pendulum position
bar=line('xdata',[0 data(1)],'ydata',...
[0 data(2)]','linewidth',3);
mass=line('xdata',data(1),'ydata',data(2),'marker','o',...
'markersize',15,'MarkerFaceColor','b');
hinge=line('xdata',0,'ydata',0,'marker','o',...
'markersize',70);
axis([-5 5 -5 5])
grid % comment this out if you do like the grid
set(gca,'Fontsize',14)
set(gca,'dataaspectratio',[1 1 1])
box on

b=0; % initial angular speed
disp('Can maximize the display so you can see the action better,')
disp('then press ENTER.')
disp('If you are happy with the display size, press ENTER.')
pause
 
v = VideoWriter('pendulum4.avi');
open(v);
% solve the diff eqns using the ode45 method
while (t1<100)
tspan = 0:0.01:100;
x0 = [a; b];
[t1,x] = ode45(@pend_l,tspan,x0);
X1 = x(:,1);
X2 = x(:,2);
    R=[cos(a) -sin(a);sin(a) cos(a)];
    datanew=R*data_init;
    % change the property values of the bar and hinge objects
    set(bar,'xdata',[0 datanew(1)],'ydata',[0 datanew(2)]);
    set(mass,'xdata',datanew(1),'ydata',datanew(2));
    set(hinge,'xdata',0,'ydata',0)
    set(gca,'dataaspectratio',[1 1 1])
    drawnow;
    % get the frame and write to the file
    frame = getframe(gcf);
    writeVideo(v,frame);
end
% close the video writer
close(v);
end

[Alan Stevens]

It's not the same! You should only need to run ode45 once, not every time within the while loop!

[KLETECH MOTORSPORTS]

I realise that there was a line missing at the bottom,

    writeVideo(v,frame);

I was wondering if it was possible to have two pendulums oscillating there instead of one?

One pendulum would follow the equation defined by

    function dXdt = odefnL(~,X)

        theta = X(1);
        thetadot = X(2);
        dXdt = [thetadot;
               -5*(theta)];    %theta instead of sin(theta)
end

and the other would follow

function dYdt = odefnN(~,Y)

        theta = Y(1);
        thetadot = Y(2);
        dYdt = [thetadot;
               -5*sin(theta)];
end

I tried running the code you provided with the extra function, but i only got a single pendulum instead of two in the avi file.

[Alan Stevens]

This should do it (I've assumed the two start in the same place, but this is easily altered):

% A script to animate the motion of the simple pendulum
theta0=deg2rad(50);
tfinal=10;
thetadot0=0; % initial angular speed
tspan = [0 tfinal];
IC = [theta0 thetadot0 theta0 thetadot0];

[t, Y] = ode45(@odefn, tspan,IC);
theta = Y(:,1);
thetadot = Y(:,2);
theta2 = Y(:,3);
theta2dot = Y(:,4);

data_init=[0 0; -4 0]; %pendulum length
R=[cos(theta(1)) -sin(theta(1));sin(theta(1)) cos(theta(1))];
data=R*data_init; % initial pendulum position
R2=[cos(theta2(1)) -sin(theta2(1));sin(theta2(1)) cos(theta2(1))];
data2=R2*data_init; % initial pendulum position

bar=line('xdata',[0 data(1)],'ydata',...
[0 data(2)]','linewidth',3);
mass=line('xdata',data(1),'ydata',data(2),'marker','o',...
'markersize',15,'MarkerFaceColor','b');
hinge=line('xdata',0,'ydata',0,'marker','o',...
'markersize',7);

bar2=line('xdata',[0 data(1)],'ydata',...
[0 data(2)]','linewidth',3);
mass2=line('xdata',data(1),'ydata',data(2),'marker','o',...
'markersize',15,'MarkerFaceColor','g');


axis([-5 5 -5 5])
grid % comment this out if you do like the grid
set(gca,'Fontsize',14)
set(gca,'dataaspectratio',[1 1 1])
box on

for i=2:numel(t)
    R=[cos(theta(i)) -sin(theta(i));sin(theta(i)) cos(theta(i))];
    datanew=R*data_init;
    set(bar,'xdata',[0 datanew(1)],'ydata',[0 datanew(2)]);
    set(mass,'xdata',datanew(1),'ydata',datanew(2));
    set(hinge,'xdata',0,'ydata',0)
    set(gca,'dataaspectratio',[1 1 1])
    
    R2=[cos(theta2(i)) -sin(theta2(i));sin(theta2(i)) cos(theta2(i))];
    datanew2=R2*data_init;
    set(bar2,'xdata',[0 datanew2(1)],'ydata',[0 datanew2(2)]);
    set(mass2,'xdata',datanew2(1),'ydata',datanew2(2));
    set(hinge,'xdata',0,'ydata',0)
    set(gca,'dataaspectratio',[1 1 1])
   
    drawnow;
    frame = getframe(gcf);
   
    
end

function dYdt = odefn(~,Y)

        theta = Y(1);
        thetadot = Y(2);
        theta2 = Y(3);
        theta2dot = Y(4);
        dYdt = [thetadot;
               -5*sin(theta);
               theta2dot;
               -5*theta2];
end

[KLETECH MOTORSPORTS]

Hi, i tried it, and it worked!

1) Any idea what i can do to change the time step in tspan? The pendulums are oscillating way too fast to understand the difference in motion.

2) Also, about this bit of the code:

bar=line('xdata',[0 data(1)],'ydata',...
[0 data(2)]','linewidth',3);
mass=line('xdata',data(1),'ydata',data(2),'marker','o',...
'markersize',15,'MarkerFaceColor','b');
hinge=line('xdata',0,'ydata',0,'marker','o',...
'markersize',7);

are 'bar', 'mass', and 'hinge' inbuilt variables/functions? or just variables?

3) And here:

for i=2:numel(t)
R=[cos(theta(i)) -sin(theta(i));sin(theta(i)) cos(theta(i))];

Why do we take theta as a function of time, 't', when we have already taken theta(1) and theta2(1) in the two cases respectively?

data_init=[0 0; -4 0]; %pendulum length
R=[cos(theta(1)) -sin(theta(1));sin(theta(1)) cos(theta(1))];
data=R*data_init; % initial pendulum position
R2=[cos(theta2(1)) -sin(theta2(1));sin(theta2(1)) cos(theta2(1))];
data2=R2*data_init; % initial pendulum position

[Alan Stevens]

1) Add

pause(0.05)

immediately after

frame = getframe(gcf);

2) "bar", "mass" and "hinge" are variables you have defined in terms of the inbuilt function "line". The fact that you ask the question suggests you didn't write the code, but have picked up coding written by someone else.

3) The angles theta and theta2 are changing all the time, so you need the new angles to reposition the masses.

[KLETECH MOTORSPORTS]

No, i have not written the code. I got it from Rudrapratap in part, and someone here was kind enough to help me modify it a bit. Is that a problem?

[Alan Stevens]

Not a problem!

[KLETECH MOTORSPORTS]

Thank you for your help. It really did clear up some concepts for me. I've only been using matlab for about a month now, so i really appreciate your help. Thanks Alan.

[KLETECH MOTORSPORTS]

Another doubt i had was, what are the values of data(1) and data(2) and where do they come from?

bar=line('xdata',[0 data(1)],'ydata',...
[0 data(2)]','linewidth',3);
mass=line('xdata',data(1),'ydata',data(2),'marker','o',...
'markersize',15,'MarkerFaceColor','b');
hinge=line('xdata',0,'ydata',0,'marker','o',...
'markersize',7);

Why are datanew1 and datanew2 defined?

[Alan Stevens]

The matrix data_init specifies the positions of the pivot point (0,0) and the pendulum mass at rest (0, -4). The R value is a rotation matrix that rotates the line between the pivot and the mass. Initially, R uses theta(1), so the pendulum immediatelyrotates to 50 degrees from the vertical. The value of theta then changes every timestep, so R changes and the rotates the pendulum to its new position (the new positions are ste in datanew).