ch2example3prg1.m

来自「一些有用的MATLAB通信信号处理的源代码程序」· M 代码 · 共 27 行

% ch2example3prg1.m
dt=0.0001;          % 仿真步进
T=15;               % 仿真时间长度
t=0:dt:T;           % 仿真计算时间序列
g=9.8;              % 重力加速度
L=1;                % 摆线长度
m=10;               % 摆锤质量
k=5;                % 空气阻力比例系数
theta0=3.1;         % 初始摆角设置
v0=0;               % 初始摆速设置
v=zeros(size(t));   % 程序存储变量预先初始化,可提高执行速度
theta=zeros(size(t));
v(1)=v0;            % 初始值赋值
theta(1)=theta0; 
for n=1:length(t)   % 仿真求解开始
    v(n+1)=v(n)+(g*sin(theta(n))-k./m.*v(n)).*dt;
    theta(n+1)=theta(n)-1./L.*v(n).*dt;
end
% 使用双坐标系统来作图, 注意作图和图形标注的技巧
[AX,H1,H2] = plotyy(t,v(1:length(t)),t,theta(1:length(t)),'plot');
set(H1,'LineStyle','-');                            % 设置作图线型
set(H2,'LineStyle','-.');
set(get(AX(1),'Ylabel'),'String','线速度 v(t) m/s');% 坐标标注
set(get(AX(2),'Ylabel'),'String','角位移 \theta(t) rad ');
xlabel('时间 t '); 
legend(H1,'线速度 v(t)',2);
legend(H2,'角位移 \theta(t)',1);