function xdot=pendulumstateeq(t, x, flag, par)
% 考慮空氣阻力的單擺系統(tǒng)的狀態(tài)方程函數(shù)
% 輸入: t 當(dāng)前計(jì)算時(shí)刻,  flag 此處不用
% x為2#1矩陣, x(1)為當(dāng)前擺錘的切向速度; x(2)為當(dāng)前角位移
% par為4#1系統(tǒng)參數(shù)矩陣
% par(1) 為重力加速度參數(shù) g
% par(2) 為空氣阻力比例系數(shù) k
% par(3) 為擺錘質(zhì)量 m
% par(4) 為擺線長(zhǎng)度 L
xdot=zeros(2,1);                    % 狀態(tài)變量矩陣初始化
g=par(1);k=par(2);m=par(3);L=par(4);% 系統(tǒng)參數(shù)設(shè)置
xdot(1)=g*sin(x(2))-k/m*x(1);       % 狀態(tài)方程組
xdot(2)=-1/L*x(1);