代碼允許轉(zhuǎn)載，但是放在附錄時一定要記得把出處和制作人標(biāo)清楚?。。。。。。。。。。。。。。。。。。。。。。。。。?！

常見問題：

1.代碼復(fù)制一定要復(fù)制到編輯器中，而不是命令行?

2.一定保證符號全英文

?3.對代碼修改后尤其是方程部分修改后，求解錯誤是正常現(xiàn)象，請一定要首先檢查方程是否正確?

4.這個代碼求解方程是龍格庫塔法求解的方程，如果不了解數(shù)值計算的可以先了解一下數(shù)值計算！！

軟件：matlab

代碼：

%--------------

%2020iypt下沉的碗仿真動畫

%使用ode45求解動力學(xué)方程后進(jìn)行動畫模擬

%這里的模型使用的是無限長的碗，圖畫的為有限長所以會有沒過邊界后還在下降情況

%作者：多一份問候 時間：2020/1/22 僅在b戰(zhàn)發(fā)布

%引用或轉(zhuǎn)載請標(biāo)明出處

%----------------

global M g p R r hd N N2 L

M=0.13581;%千克

g=9.8;

p=1000;%水的密度

R=0.05;

r=0.003;

hd=0.002;%碗的厚度

N=0.0001;%阻尼系數(shù)小孔的

N2=1;%碗的阻力系數(shù)

L=0.1;%碗的高度

t1=[0,14];

H1=[0;0;0;0];

[t,H] = ode45(@dH,t1,H1);%求解動力學(xué)方程

% hold on

% plot(t,H(:,1:3));%繪制的曲線隨時間的變化圖

% plot(t,H(:,1));%碗的底端高度

% plot(t,H(:,2));%碗的下落速度

% plot(t,H(:,3));%碗內(nèi)液面高度

%動畫繪制

%在上述運行完畢后得到

A=size(H);

loops = A(1,1);

F(loops) = struct('cdata',[],'colormap',[]);

v = VideoWriter('peaks135.31kg0.5_0.003.avi');%創(chuàng)建動畫

open(v);

for j=1:loops

Hz{j}=num2str(H(j));

end

for i=1:A(1,1)

figure(i)

hold on

axis([-0.25 0.25 -0.4 0.1])%設(shè)置坐標(biāo)軸范圍x：-0.25，0.25 y：-0.4 0.1

%碗的圖形繪制

%這里的圖形繪制認(rèn)為是隨著頁面的高度變畫而變化的

plot([R,R],[L-H(i,1),-H(i,1)],'LineWidth',2,'Color','b');%右側(cè)的碗

plot([-R,-R],[L-H(i,1),-H(i,1)],'LineWidth',2,'Color','b');%左側(cè)的碗

plot([-R,-r],[-H(i,1),-H(i,1)],'LineWidth',2,'Color','b');%碗底部的小孔左側(cè)

plot([r,R],[-H(i,1),-H(i,1)],'LineWidth',2,'Color','b');%碗底部的小孔右側(cè)

plot([-R,R],[-H(i,1)+H(i,3),-H(i,1)+H(i,3)],'LineWidth',2,'Color','b');%碗內(nèi)液面

plot([-0.25,R],[-H(i,1)+H(i,3),-H(i,1)+H(i,3)],'LineWidth',0.5,'Color','b','LineStyle','--');%碗內(nèi)液面指示線

text(-0.25,-H(i,1)+H(i,3)-0.006,Hz{1,i});%碗內(nèi)液面指示線數(shù)字

plot([-0.25,R],[L-H(i,1),0.1-H(i,1)],'LineWidth',0.5,'Color','b','LineStyle','--');%碗內(nèi)液面指示線

text(-0.25,L-H(i,1)-0.006,Hz{1,i});%碗內(nèi)液面指示線數(shù)字

%水面高度認(rèn)為是不變的

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

plot([R,0.25],[0,0],'LineWidth',2,'Color','c');%頁面高度右側(cè)

F(i) = getframe(gcf);%保存幀

writeVideo(v,F(i));%寫入動畫

close

end

close(v);

function dHdt=dH(t,H)

global M g p R r hd N N2

dHdt=zeros(4,1);

dHdt(1)=H(2);

dHdt(2)=g-p*g*pi*(H(1)-H(3))*(R^2-r^2)/M-N2*H(2);

dHdt(3)=sqrt(-2*g*r^4/R^4*(H(3)-H(1))-N*H(4)*hd);

dHdt(4)=H(3);

end