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
clcclearclose 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=50else theta;endtheta=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=10else tfinal;enddata_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 ondt=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))];endfunction xdot = pend_l(t1,x)wsq=3.5; xdot = [x(2); -wsq*x(1)];end
Best Answer