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

常見問題：

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

2.一定保證符號全英文

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

4.代碼所用方法為歐拉法，精度并不高，個人思考難免有發(fā)現(xiàn)不了的錯誤，還請多多包涵。更好的方法歡迎評論區(qū)討論

所用軟件：MATLAB

代碼：

%b站多一份問候

%iypt2021第七題，求解方法為歐拉法。

%方程求解代碼，有錯誤請聯(lián)系作者。

clear all

R=0.2; % m質(zhì)量

w=80; % rad/s大環(huán)的角速度

T=1000;%計算時間

Q=zeros(T,1); % angle角度

Q(1)=0.01;%角度初始值，為零的話珠子不會抬起

u=0.4; % friction動摩擦系數(shù)

Q1=zeros(T,1); % angle derivative角速度

Q2=zeros(T,1); % angular acceleration角加速度

g=9.8;%重力加速度

t=0.001;%計算時間間隔

for i=1:T

Q2(i+1)=(R*sin(Q(i))*w^2-(R*cos(Q(i))*Q1(i)^2+g+u*2*Q1(i)*R*cos(Q(i))*w)*(sin(Q(i))+u*cos(Q(i)))/(cos(Q(i))-u*sin(Q(i)))-u*2*Q1(i)*R*cos(Q(i))*w+R*sin(Q(i))*Q1(i)^2)/(R*cos(Q(i))+R*sin(Q(i))*(sin(Q(i))+u*cos(Q(i)))/(cos(Q(i))-u*sin(Q(i)))) ;%方程

Q1(i+1)=Q1(i)+Q2(i+1)*t;%角速度

Q(i+1)=Q(i)+Q1(i+1)*t;%速度

end

r=R*sin(Q);%水平位移

h=R-R*cos(Q);%垂直偏離

hold on

plot(Q)%繪圖

珠子動力學(xué)文獻(xiàn)合集：

文獻(xiàn)鏈接: https://pan.baidu.com/s/1_zlXjFJ1Ry4Lv273ltsuOg 提取碼: grg9 復(fù)制這段內(nèi)容后打開百度網(wǎng)盤手機(jī)App，操作更方便哦