abc.m

来自「something about phisics」· M 代码 · 共 35 行

m=1;
b=[0,0.1,0.1];
p=[0,0,1];
px=[2.1;2.1;2.1];
py=[5.5;4.4;1.8];
strdd{1}='无阻尼';
strdd{2}='阻尼正比于v';
strdd{3}='阻力正比于v^2';
a=['r';'g';'b'];
figure
for i=1:3
[t,y]=ode45('abcfun',[0:0.01:5],[0,5,0,2],[ ],b(i),p(i),m);
H{i}=max(y(:,3))
T{i}=t(find(y(:,3)==H{i})) 
vx0{i}=y(end,2)
subplot(3,1,1)
axis([0 10 -50 3]);
hold on
xlabel('x')
ylabel('y')
comet(y(:,1),y(:,3));
subplot(3,1,2)
axis([0 5 0 6])
hold on
xlabel('t')
ylabel('dx/dt')
text(px(i),py(i),strdd{i});
plot(t,y(:,2),a(i))
subplot(3,1,3)
axis([0 5 -60 2])
hold on
xlabel('t')
ylabel('dy/dt')
plot(t,y(:,4),a(i))
end