2022数学建模A动态仿真

%%  2022数学建模A动态仿真

clc;clear;

figure('numbertitle','off','name','2022数学建模A动态仿真');

%赋值

global m M K C W C1 P F

m=4800;

M=2400;

K=80000;

C=10000;

W=1.8;

C1=656;

P=10000;%改变

F=6000;

%ode45求解绘图

t1=[0:0.1:40];H1=[0;10;0;5];

[t,H]=ode45(@dH,t1,H1);

R=6.5;L=18;L1=15;r=1;Rr=1.1;L2=2;A=1.3;

for i=1:1:length(H)

clf;xlim([-15 15]);ylim([-22 22]);hold on%设置坐标轴范围

%外边

plot([-R,R],[L-H(i,1)-8,L-H(i,1)-8],'LineWidth',2,'Color','k');%

plot([R,R],[L-H(i,1)-8,-H(i,1)-8],'LineWidth',2,'Color','k');%右侧

plot([-R,-R],[L-H(i,1)-8,-H(i,1)-8],'LineWidth',2,'Color','k');%左侧

plot([-R,0],[-H(i,1)-8,-H(i,1)-4-8],'LineWidth',2,'Color','k');%底部右侧

plot([0,R],[-H(i,1)-4-8,-H(i,1)-8],'LineWidth',2,'Color','k');%底部左侧

plot([-R,R],[-H(i,1)-8,-H(i,1)-8],'LineWidth',2,'Color','k');%底部中间

plot([-15,-R],[-H(i,1)-8,-H(i,1)-8],'LineWidth',0.5,'Color','b','LineStyle','--');

%里边

plot([-r,r],[L1-H(i,1)-8,L1-H(i,1)-8],'LineWidth',2,'Color','g');

plot([r,r],[L1-H(i,1)-8,-H(i,1)-8],'LineWidth',2,'Color','g');

plot([-r,-r],[L1-H(i,1)-8,-H(i,1)-8],'LineWidth',2,'Color','g');

%里边正方形

plot([-Rr,Rr],[L2-A*H(i,3),L2-A*H(i,3)],'LineWidth',2,'Color','r');

plot([Rr,Rr],[L2-A*H(i,3),-A*H(i,3)],'LineWidth',2,'Color','r');

plot([-Rr,-Rr],[L2-A*H(i,3),-A*H(i,3)],'LineWidth',2,'Color','r');

plot([-Rr,Rr],[-A*H(i,3),-A*H(i,3)],'LineWidth',2,'Color','r');

plot([-15,-Rr],[-A*H(i,3),-A*H(i,3)],'LineWidth',0.5,'Color','b','LineStyle','--');

%水面

plot([-15,-R],[0,0],'LineWidth',2,'Color','c');%液面高度左侧

plot([R,15],[0,0],'LineWidth',2,'Color','c');%页面高度右侧

hold off;pause(0.05);%可以在这里添加输出动图的程序

print(['Frame ' num2str(i)], '-dpng', '-r150');

drawnow;%实时更新坐标图窗口

end

 

GifName = '2022数学建模A动态仿真.gif';

delay = 0.5;    % Delay between frames (s)

for ii = 1:10:length(t);

    [A, ~] = imread(['Frame ' num2str(ii) '.png']);

    [X, map] = rgb2ind(A, 256);

    if ii == 1

        imwrite(X, map, GifName, 'gif', 'LoopCount', inf, 'DelayTime', delay)

    else

        imwrite(X, map, GifName, 'gif', 'WriteMode', 'append', 'DelayTime', delay)

    end

end

 

%动力学方程

function dHdt=dH(t,H)  

global m M K C W C1 P F 

dHdt=zeros(4,1);

dHdt(1)=H(2);

dHdt(2)=(F*cos(W*t)-C1*H(2)-P*H(1)+K*(H(1)-H(3))+C*(H(2)-H(4)))./m;

dHdt(3)=H(4);

dHdt(4)=(K*(H(1)-H(3))+C*(H(2)-H(4)))./M;

end