【数学建模】基于matlab单摆运动仿真【含Matlab源码 3997期】.zip

%% 单摆运动仿真 % 求解：1.单摆周期与角振幅之间关系 % 2.演示单摆振动 % 单摆系统参数 % 已知：单摆的摆锤质量m，角位置theta，摆杆长L % 小角度振动时,单摆周期与角振幅无关，具有等时性，周期为T0；大角度时周期受角振幅影响，为T % 这里简化分析起见，这里计算出T_dot=T/T0，T_dot是约化周期,其中T0=2*pi*sqrt(L/g)，则有小角度振动T_dot≈1. % 通过T_dot=T/T0的计算，可以观察出多大角度下单摆的周期会产生收角振幅的影响 % 假设：摆杆无质量，摆球为质点 % 建立坐标系：摆杆铰接点为原点，水平向右为x轴正向，竖直向上为y轴正向, % 单摆的角位置为摆杆与竖直向下方向的夹角，单摆在右时角位置为正，在左时角位置为负 % 并设单摆逆时针旋转时为旋转正方向 % 单摆振动的周期和运动规律由摆锤的运动方程得出，详细解析见以下网页： Judge = input('是否打开单摆运动规律的网页讲解，1--打开，0--跳过： '); if Judge == 1 url = 'https://wenku.baidu.com/view/29ec28ee81c758f5f61f6738.html'; web(url) end %% 1.单摆的周期与振幅的关系 % 单摆的周期 ( 用内联函数数值积分和符号积分与椭圆积分比较 ) clc; clear all;close all; theta = 1:179; % 角振幅向量（1°-179°）,假设摆杆为轻质杆。 thm = theta*pi/180; % 角振幅由角度化为弧度数 % 用第一类完全椭圆积分计算单摆精确周期T_dot % matlab中第一类完全椭圆积分函数语法为：T_dot = ellipke(M) % 见网页运动规律解析，参数代入时，M = k^2, k = sin(thm/2)，T_dot = T/T0 =ellipke(M)*2/pi. T_dot = ellipke(sin(thm/2).^2)*2/pi; % 用以上网页讲解式中的公式（5.1.5），构造分别内联函数和匿名函数句柄进行积分计算约化周期T_dot T1 = []; % 周期向量置空 T2 = []; % 周期向量置空 n = length(thm); % 角度个数 angle = 1:5:n; % 下标向量 ( 每隔 5度取一个角度 ) for iterNum = angle % 按下标循环 s = ['1./sqrt(cos(x)-cos(' num2str(thm(iterNum)) ') )']; % 式（5.1.5）被积函数字符串表达式 f = inline(s); % 构造被积函数的内联函数 t = quadl(f,0,thm(iterNum))*sqrt(2)/pi; %用数值积分函数quadl按式（5.1.5）求周期T/T0 T1=[T1,t]; % 连接周期，即将各个角振幅对应的单摆周期值依次存入向量T1中 % 创建匿名函数,其中 % 内层函数g = (double(integral(@(x) 1./sqrt(cos(x)-cos(thetam)),0,thm(iterNum)))*sqrt(2)/pi) % 内层函数利用integral积分函数在0-thm(iterNum)之间进行符号积分，积分变量为x,thetam为参数 % 外层函数参数h = @(thetam) (g) 将角振幅参数thetam = thm(iterNum)代入符号积分中，计算约化周期 h = @(thetam) ( (double(integral(@(x) 1./sqrt(cos(x)-cos(thetam)),0,thm(iterNum)))*sqrt(2)/pi) ); t = h(thm(iterNum)); % 用符号积分求周期 T2=[T2,t]; % 连接周期,各个角振幅对应的单摆周期值依次存入向量T2中 end % 结束循环,得到各个角振幅下计算得到的周期向量 fig1 = figure; % 创建图形窗口 plot(theta,T_dot,theta(angle),T1, '.' ,theta(angle),T2, 'o' ) % 画周期曲线 grid on % 加网格 fs=16; % 字体大小 xlabel( '\it角振幅\theta\rm_m/\circ' , 'FontSize' ,fs) % x标签 ylabel( '\it约化周期T/T\rm_0' , 'fontsize' ,fs) % y标签 title( ' 单摆的周期与角振幅的关系 ' , 'FontSize' ,fs) % 标题 legend( ' 椭圆积分 ' , ' 数值积分 ' , ' 符号积分 ') % 图例 text(0,2, '\itT\rm_0=2\pi(\itL/g\rm)^{1/2}' , 'FontSize' ,fs) % 小角单摆周期 % 所得图中可以明显观察出单摆周期与角振幅之间关系 % 图形打印 print(fig1,'单摆的周期与角振幅的关系','-dpng','-r1000') % png格式、分辨率为1000，打印图形 %% 运动仿真 % 解 ： % 利用以上假设，由牛顿第二定律以及法相和切向受力分解（考虑阻尼效应）： % ∑Fn = m*an an=L*w^2 % ∑Fn法向合力，an法向加速度，L单摆长，w为单摆加速度 % T-m*g*cos(θ)=m*L*w^2 得 T=m*L*w^2 +m*g*cos(θ) % T摆杆拉力，g重力加速度，θ单摆角位置， % 单摆处于右侧且在上升时，依据假设：旋转方向为正，角位置为正，重力在切向方向上地分力和阻尼力均为阻力，根据力矩平衡则有 % （-m*g*sin(θ) - c*w）*L = ∑Ft*L = m*L^2*α % 上式中，∑Ft为切向合力，c为阻尼系数，θ_d1代表角速度，α为角加速度，m*L^2为摆球旋转地转动惯量 % 化简得到： α= -g*sin(θ)/L - c/(m*L) * w % 此式就是单摆运动微分方程（当不考虑阻尼时，令c=0即可） % 记r = g/L; k = c/(m*L);微分方程化简为: % α= -r*sin(θ) - k * w % 记向量x = [θ; w]，称为单摆地状态向量，x(1)=θ,x(2)=w % dx(1)/dt = dθ/dt = w 即 dx(1)/dt = x(2) % dx(2)/dt = α= -r*sin(θ) - k * w = -r*sin(x(1)) - k*x(2) % ∴ dx/dt = [x(2); -r*sin(x(1)) - k*x(2)] 该式成为单摆的状态方程 % 建立匿名函数：dynamics_pendulum = @(t,x) [x(2);-k*x(2) - r*sin(x(1))]; % 在规定的仿真时间如0到tSpan秒内，初始状态向量init = [theta0; w0],利用matlab的ode45函数对以上微分方程积分 % 使用ode45函数的语法为：sol = ode45(dynamics_pendulum,[0 tSpan],init); % 积分函数传入了微分方程句柄的名字dynamics_pendulum；积分的始末时间：0和设定的tSpan；初始状态init % 积分之后函数返回 sol 这个结构体变量，结构体变量sol里面存储着[0 tSpan]这个时间段内一定数量的积分结果，为了得到更多的积分点 % 利用t = linspace(0,tSpan,1000);Y = deval(sol,t);这两个命令， % 先在仿真区间[0 tSpan]内创建1000个时间点，Y = deval(sol,t)传入积分结果结构体和时间变量t，返回Y。 % Y是一个2×1000的矩阵，第一行为单摆角位置解向量，第二行为单摆角速度解向量，该解是由deval函数插值得来，以满足仿真时尽量多的时间点 m = 0.3; % (kg) 摆球质量 L = 0.5; % （m）杆长 g = 9.81; % (m/s^2) 重力加速度 judge_damping = input('是否需要考虑阻尼对单摆影响？考虑---1，不考虑---0: '); if judge_damping == 1 c = input('请输入阻尼系数(如c=0.05): '); % （N*s/m）阻尼系数 else c = 0; end % 简化的单摆二阶微分方程系数 r = g/L; k = c/(m*L); % 单摆微分方程(状态空间方程) dynamics_pendulum = @(t,x) [x(2);-k*x(2) - r*sin(x(1))]; th0 = input('请输入初始角振幅(如90度): '); % （度）单摆相对-j轴的初始角度 th0 = th0*pi/180; % 度化为弧度 w0 = 0.0; % (rad/s) 单摆初始角速度 init = [th0; w0]; % 初始状态向量 tSpan = input('请设置仿真时间（s）,如15秒: '); % 仿真时间s % 解微分方程，时间范围设置为0-tSpan s sol = ode45(dynamics_pendulum,[0 tSpan],init); t = linspace(0,tSpan,1000); Y = deval(sol,t); theta = Y(1,:); % 角度向量 omega = Y(2,:); % 角速度向量 nStep = length(t); % 时间向量长度 % 画出解 fig2 = figure; % 新建图窗 plot(t,theta,'r'); grid on; hold on; plot(t,omega,'k'); xlabel('Time (s)') ylabel('Amplitude') title('单摆振动角位移和角速度曲线') legend({'Angle (rad)','Angular speed (rad/s)'}) % 图形打印 if judge_damping == 1 print(fig2,'有阻尼单摆角位移和角速度曲线','-dpng','-r1000'); else print(fig2,'无阻尼单摆角位移和角速度曲线','-dpng','-r1000'); end %% 仿真动画生成 % 建立视频 whether_Generate_video = input('是否需要生成视频（视频生成会捕捉每一帧图片，导致运行速度稍微变慢），是---> 1，否---> 0： '); if whether_Generate_video == 1 % v = VideoWriter(filename) 创建一个 VideoWriter 对象以将视频数据写入采用 Motion JPEG 压缩技术的 AVI 文件。 disp('单摆仿真动画运行中，同时正在生成视频，请稍后...') if judge_damping == 1 v = VideoWriter('有阻尼的单摆运动仿真.avi'); else v = VideoWriter('无阻尼的单摆运动仿真.avi'); end % 视频的播放速率（每秒帧数），指定为正数(默认30)。调用 open 之后，无法更改 FrameRate 值。这里指定为40倍播放 v.FrameRate = 40; open(v); % 打开该对象以进行写入 end % 动画初始化: 画单摆的初始位置，画单摆角位置和角速度的第一个点，并获得各个动画对象的句柄 figure('units','normalized','outerposition',[0 0 1 1]) % 创建2行3列总共6个子图，其中单摆运动子图占据第1、2、4、5子图这四个位置，第3个子图为单摆角位置曲线图，第6个为角速度曲线图。 subplot(2,3,[1:2,4:5]) axis equal % 显示单摆动画图窗的横纵坐标为等刻度 axis([-1.3*L 1.3*L -1.3*L 1.3*L]) % 设置图窗范围 % Title_handle1为单摆动画图窗标题句柄，标题中更新仿真时间 Title_handle1 = title( ['单摆振动仿真: ','共 ',num2str(t(end)),' s','第 ',num2str(t(1)),' s'], 'FontSize' ,12); % 显示标题 hold on % 保持图像 % Suspension_point为悬点句柄 Suspension_point = plot(0,0,'Marker','o','MarkerSize',5,'MarkerFaceColor','k','MarkerEdgeColor','k'); %画悬点 % Dotted_line为竖直虚线句柄 Dotted_line = plot([0;0],[0;-1.1*L], '-.' ); % 画竖虚线 Ball_x = L*sin(theta(1)); % 摆球的起始横坐标, Ball_y = -L*cos(theta(1)); % 摆球的起始纵坐标 % Swing_bar为摆杆句柄 Swing_bar = plot([0 Ball_x],[0 Ball_y] ,'b','