含阻尼的单摆仿真

闲的蛋疼ing

公式的推导，用欧拉法，代码如下，虽有瑕疵但能运行，就这么着吧

clc;clear;close;
figure ('name','含阻尼的单摆');
xlabel('时间t/s');
hold on
dt=0.0001;%步长
T=20;%仿真时长
k=1.5;m=8;
t=0:dt:T;
g=9.8;L=1.5;th0=1.5;v0=0;
v=zeros(size(t));
th=zeros(size(t));
v(1)=v0;
th(1)=th0;
for i=1:length(t)
  v(i+1)=v(i)+(g.*sin(th(i))-k/m.*v(i)).*dt;
  th(i+1)=th(i)-(1/L).*v(i).*dt;
end
[x,y1,y2]=plotyy(t,v(1:length(t)),t,th(1:length(t)),'plot');
set(y1,'LineStyle','--');
set(y2,'LineStyle','.');
hold on
legend(y1,'线速度v',2);
legend(y2,'角速度ω',1)
set(get(x(1),'Ylabel'),'string','线速度v');
set(get(x(2),'Ylabel'),'string','角速度v')；

执行结果