resonance_for_phybook

% 高中物理——演示共振现象——数值模拟
% 新版人教高中物理 选择性必修1
% 第二章 第6节 受迫振动 共振 51页 图2.6-3 演示共振现象
% 从数值模拟结果来看此装置不适合作为共振的演示
clc;clear;close all
%% 模拟参数与数值计算
g=10;                          % 重力加速度
L=1;                           % 基准摆长
m=1;                           % 基准质量
h=0.1*L;                       % 横梁高度
M=0;                           % 横梁质量
Li=[1 1.2 1.4 1 .6 .8 1]'*L;   % 各绳长度
mi=[1   1   1 1  1  1 1]'*m;   % 各球质量
theta=1;                       % 大质量小球A的释放角度
T=190;                         % 模拟时间
num=size(mi,1);                % 小球数量
% 各球初始角度
thetai=zeros(num,1);% 各球初始角度为0
thetai(4)=1;% 大质量小球的释放角度为1，约60°
alpha=0;% 横梁的初始角度
pa=0;p_i=zeros(num,1); % 初始静止
u0=[alpha;pa;thetai;p_i];
ti=0:0.01:T;
odeopts=odeset('RelTol',1e-10,'AbsTol',1e-12, ...
    'MaxStep',1e-2,'InitialStep',1e-4);
[~,u] = ode45(@(t,u) gz34(t,u,g,m,h,M,Li,mi,num),ti,u0,odeopts);
%% 几何
% 小球位置(x,y,z)以及对应摆绳结位置(xo,yo,zo)
x= Li'.*sin(u(:,3:num+2))+h*sin(u(:,1));
z=-Li'.*cos(u(:,3:num+2))-h*cos(u(:,1));
xo=zeros(size(x))+h*sin(u(:,1));
zo=zeros(size(z))-h*cos(u(:,1));
yo=zeros(size(xo))+(0.5:num-0.5)*2*L/num-L;
y=yo;
% 横梁及横梁悬线位置
y1=-0.9*L;y2=0.9*L;
x1=h*sin( u(:,1));x2=x1;
z1=-h*cos( u(:,1));z2=z1;
y1=y1+zeros(size(x1));
y2=y2+zeros(size(x1));
%% 绘图
% figure;
dpi=1.25;scr=[1920 1080];res=[960 720];
fig=gcf;fig.Position=[(scr-res)/2,res]/dpi;
hold on
% 控制显示区域
plot3(-L,-0.5*L,-1.8*L,L,0.5*L,0)
axis tight;axis equal
box on
set(gca,'XTick',[],'YTick',[],'ZTick',[]);
ax=gca;ax.Clipping='off';
ax.CameraViewAngle = 10;%视野范围
ax.CameraTarget = [0 0 -0.9];%目标位置
ax.CameraPosition = [0.2 10 -0.8];
% 支架
plot3([0 0 0 0],[L L -L -L],[-1.8*L 0 0 -1.8*L],...
    '-','Color','k','LineWidth',2)
% 横梁绳子
n=1;
line1=line([0 x1(n)],[y1(n) y1(n)],[0 z1(n)],...
    'linestyle','-','Color','k','LineWidth',1);
line2=line([0 x2(n)],[y2(n) y2(n)],[0 z2(n)],...
    'linestyle','-','Color','k','LineWidth',1);
% 横梁
line3=line([x1(n) x1(n)],[y1(n) y2(n)],[z1(n) z1(n)],...
    'linestyle','-','Color','k','LineWidth',3);
% 颜色集
color_set(1,:)=[0.8500 0.3250 0.0980];
color_set(2,:)=[     0 0.4470 0.7410];
color_set(3,:)=[0.4940 0.1840 0.5560];
color_set(4,:)=[     0      0      0];
color_set(5,:)=[0.9290 0.6940 0.1250];
color_set(6,:)=[0.4660 0.6740 0.1880];
color_set(7,:)=[0.8500 0.3250 0.0980];
% 画小球
ms=40*ones(7,1);ms(4)=60;
for ii=1:num
    sp(ii)=plot3(x(n,ii),y(n,ii),z(n,ii),'.', ...
        'Color',color_set(ii,:),'markersize',ms(ii));
end
% 画绳子
for ii=1:num
    lines(ii)=line( ...
        [xo(n,ii) x(n,ii)], ...
        [yo(n,ii) y(n,ii)], ...
        [zo(n,ii) z(n,ii)], ...
        'Color',color_set(ii,:),'LineWidth',2);
end
%% 动画
mv = 0;% 是否录制视频，1录制，其它不录制
if mv == 1
v = VideoWriter('D:\Programs\MATLAB\物理课本共振模型模拟.mp4','MPEG-4');
v.Quality = 100;v.FrameRate = 60;
open(v);
end
for n=1:length(ti)
    for ii=1:num
        sp(ii).XData=x(n,ii);
        sp(ii).YData=y(n,ii);
        sp(ii).ZData=z(n,ii);
        lines(ii).XData=[xo(n,ii) x(n,ii)];
        lines(ii).YData=[yo(n,ii) y(n,ii)];
        lines(ii).ZData=[zo(n,ii) z(n,ii)];
    end
    line1.XData=[0 x1(n)];
    line1.YData=[y1(n) y1(n)];
    line1.ZData=[0 z1(n)];
    line2.XData=[0 x2(n)];
    line2.YData=[y2(n) y2(n)];
    line2.ZData=[0 z2(n)];
    line3.XData=[x1(n) x1(n)];
    line3.YData=[y1(n) y2(n)];
    line3.ZData=[z1(n) z1(n)];
    drawnow limitrate
    if mv==1;F=getframe(fig);writeVideo(v,F);end
end
if mv == 1;close(v);end
%% 微分方程（哈密顿方程）
function du = gz34(~,u,g,m,h,M,Li,mi,num)
a=u(1);pa=u(2);
ti=u(3:num+2);
p_i=u(num+3:2*num+2);
S=sum(mi.*sin(ti-a).^2)/m+M;
da=(pa-sum(p_i.*cos(ti-a)))/S/h;
dti=(p_i-h*da.*cos(ti))./Li;
ais=da.*dti.*sin(ti-a);
dpa=sum(mi.*Li.*ais)-(M+sum(mi))*g*sin(a);
dp_i=-h*ais-g*sin(ti);
du=[da;dpa;dti;dp_i];
end