ode45 problem with plotting
5 views (last 30 days)
Show older comments
Hi all,
I am trying to run this code, aparently it has no problem. But it should give me some plots which it does not. I would be happy if you could help me out.
Best,
Argu
2 Comments
Accepted Answer
KSSV
on 14 Feb 2019
Edited: KSSV
on 14 Feb 2019
You should change view to 1. And there is small change in the name of variable..hope the below one works:
function [t,th,r,J,om]=lab1(eps,view)
% user input:
% (1) eps: this is the variable for the pendulum length: r(t)=r0*(1+eps*t); see lab sheet
% (2) view: if user enters view > 0, a window will pop up on screen to show the
% pendulum motion after the simulation is finished. It is recommended to
% set view=0 since the processing time required to show graphics on screen.
%
% program output:
% (3) t: the time array of 500,000 elements, 0, 0.001, 0.002, ... , 500
% (4) th: a 500000x2 array where each row of this array holds the record
% of the pendulum angle, theta, (d/dt)theta at time t. In particular,
% th(:,1) contains 500,000 elements of theta value at the time instants
% defined in the t array; see (3), and th(:,2) contains the 500,000 elements
% of the (d/dt)theta at the time instants defined by the t array.
% (5) r: the pendulum length as a function of time defined in t array.
% (6) J: the approximation of the total energy of the pendulum mass; see lab sheet
% (7) om: the instantaneous frequency of oscillation of the pendulum; see lab sheet
% assume pendulum mass m=1
r0=1; % pendulum length at t = 0
dt=0.001; % integration time step
tspan=0:dt:500; % all time instants for integration
th0=[0.2;0]; % initial condition: theta(0)=0.2, (d/dt)theta(0)=0
opt=odeset('Reltol',1e-9,'AbsTol',1e-9);
[t,th]=ode45(@pend,tspan,th0,opt,r0,eps);
r = r0*(1+eps*t); % pendulum length over time
J=0.5*(r.^2.*th(:,2).^2)+0.5*9.81*r.*th(:,1).^2;
om=sqrt(9.81./r);
if view
viewmotion(t,th,r);
end
function dthdt=pend(t,th,r0,eps)
r=r0*(1+eps*t); % pendulum length at time t
rdot=r0*eps; % d/dt pendulum length at time t.
dthdt(1,1)=0; % ... enter statement for (d/dt) theta
dthdt(2,1)=(9.81*sin(0.2)-2*rdot*0)/r; % ... enter statement for (d^2/dt^2) theta
end
function viewmotion(t,x,r)
clf;
plot([-1,1],[0,0],'color',[1 1 1]*.5,'linewidth',5);
hold on;
plot([0,.5],[0,.2],'color',[49 79 79]./256,'linewidth',3);
axis([-1.5 1.5 -max(r) 1]);
npt=length(t);
for i=1:100:npt
xi=r(i)*sin(x(i,1));
yi=-r(i)*cos(x(i,1));
g1=plot([0,xi],[0,yi],'color',[49 79 79]./256,'linewidth',3);
g2=plot(xi,yi,'r.','markersize',75);
pause(1e-2);
delete([g1 g2]);
end
end
end
11 Comments
More Answers (0)
See Also
Community Treasure Hunt
Find the treasures in MATLAB Central and discover how the community can help you!
Start Hunting!