ode45 problem with plotting
显示 更早的评论
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
采纳的回答
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 个评论
Argumanh
2019-2-14
KSSV
2019-2-14
There is a copy option......simply copy and paste......I should go for more mouse clicks to attach....:)
Argumanh
2019-2-14
KSSV
2019-2-14
Edited...I guess this what you want.
Argumanh
2019-2-14
Argumanh
2019-2-14
KSSV
2019-2-14
eps = 0.1 ;
view = 0 ;
[t,th,r,J,om]=lab1(eps,view) ;
The above is the way to use the function.
Argumanh
2019-2-14
KSSV
2019-2-14
Already you have th array in hand.......see the outputs of the function.
Argumanh
2019-2-14
Argumanh
2019-2-15
更多回答（0 个）
类别
在 帮助中心 和 File Exchange 中查找有关 Programming 的更多信息
Translated by 
