上面方程式在M-文件 pendulum.m中定义如下 function dot=pe %状态导数 W=2;L=6,B=0.02;g=9.81;m=W/g; dot=x(2),-B/m*x(2)-W/(m)*sin(x(1) 下面的M-文件对系统进行仿真 0=0: tfinal=s %时间间隔0~5秒 初始条件x1=1,x2=0 tol=0 0001: trace=0 %精度,不打印每一步的计算值 It, x]=ode23( pendulum, to, tfinal, xo, tol, trace); subplot(211),plot(t, x) title( Time response of a rigid pendulum,),xlabel(Time-sed text(3. 2, 3.5, Velocity,), text(3.2, 1.2, Angle-Rad') th=x(,1),w=x(;2) subplot(212), plot(th, w) lane plot of pendulum) xlabel('Position-Rad), ylabel(,Angular velocity) subplot(lll 仿真结果见图1-6 Time response of a rigid pendulum 0.5 Phase plane plot of pendulum 1-0.80.60.40.200.20.40.60.81 Position-Rad 图1-65 1 2 2 2 1 sinx mL W x m B x  x , x    上面方程式在 M-文件 pendulum.m 中定义如下： function xdot=pendulum(t, x); ％状态导数 W=2; L=.6; B=0.02; g=9.81; m=W/g; xdot=[x(2); -B/m*x(2)-W/(m*L)*sin(x(1))]; 下面的 M-文件对系统进行仿真。 t0=0; tfinal=5; ％时间间隔 0～5 秒 x0=[1,0]; ％初始条件 x1=1,x2=0 tol=0.0001; trace=0; ％精度，不打印每一步的计算值 [t, x]=ode23('pendulum',t0,tfinal,x0,tol,trace); subplot(211),plot(t, x) title('Time response of a rigid pendulum'),xlabel('Time-sec') text(3.2,3.5,'Velocity'),text(3.2,1.2,'Angle-Rad') th=x(:,1);w=x(:,2); subplot(212),plot(th, w) title('Phase plane plot of pendulum') xlabel('Position-Rad'),ylabel('Angular velocity') subplot(111) 仿真结果见图 1-6。 0 0.5 1 1.5 2 2.5 3 3.5 4 4.5 5 -4 -2 0 2 4 Time response of a rigid pendulum Time-sec Velocity Angle-Rad -1 -0.8 -0.6 -0.4 -0.2 0 0.2 0.4 0.6 0.8 1 -4 -2 0 2 4 Phase plane plot of pendulum Position-Rad Angular velocity 图 1-6