7.1轨迹举例说明:重力场中有两个物体,其中质量为m2的物体固定,而质量为m1的物体绕m2做平面圆周运动.做圆周运动的m1物体的轨道半径用变量r表示,角度用变量a表示.m1ram2两物体系统卫星绕地球转动时，m2等于地球的质量，m1等于卫星的质量，r为卫星球心与地球球心间的距离。其运动轨迹由下列方程组决定：式中：，其中t是时间变量，p为物体在地球表面做圆周运动的周期。在地球表面，r=6.373x106m。用龙格—库塔法可以实现求解：引入新状态变量：带入前面的微分方程组，可得四个一阶微分方程。建立函数文件orbit.mfunctionxd=orbit(t,x)xd=[x(2);x(1)*x(4)^2-4.0*pi^2/x(1)^2;x(4);-2.0*x(2)*x(4)/x(1)];三组初始条件(t=0)：由初始条件建立执行文件menu71.minitcond=[2001.5;1002*pi;2004];tspan=linspace(0,5,1000);options=odeset('RelTol',1e-6,'AbsTol',[1e-61e-61e-61e-6]);lintype=['-.''-.''-.'];fori=1:3[t,x]=ode45('orbit',tspan,[initcond(i,:)],options);polar(x(:,3),x(:,1),lintype(2*(i-1)+1:2*i));holdonendtext(0.5,-1.2,'椭圆轨迹');text(-1.2,1,'圆轨迹');text(1.75,2,'双曲线轨迹');7.2单自由度系统7.2.1概述一.力学模型弹簧—质量—阻尼系统其中：振体质量为m，弹簧的线性系数为k，非线性系数为a，阻尼系数为c，外力F（t）。二.运动微分方程用x表示系统的位移，则运动微分方程为：式中：固有频率非线性系数阻尼因子引入新变量转化状态空间方程形式：7.2.2线性系统的自由振动一.运动微分方程当时，得到线性振动系统的自由振动方程。二.MATLAB求解对应的函数文件FreeOcillation.mfunctionxdot=FreeOcillation(t,x,dummy,zeta)xdot=[x(2);-2.0*zeta*x(2)-x(1)];三种阻尼系数（1）阻尼系数为0.1时是欠阻尼情况（2）阻尼系数为1时是临界阻尼情况（3）阻尼系数为5时是过阻尼情况由初始条件（位移和速度均为1时）建立执行文件menu72.mzeta=[0.11.05.0];tspan=linspace(0,40,400);%生成0-40的四百个线性点lintype=['-b''--r''---r'];fori=1:3[t,x]=ode45('FreeOcillation',tspan,[11],[],zeta(i));subplot(2,1,1);plot(t,x(:,1),lintype(2*(i-1)+1:2*i));holdonsubplot(2,1,2);plot(x(:,1),x(:,2),lintype(2*(i-1)+1:2*i));holdonendsubplot(2,1,1);xlabel('Time(\tau)');ylabel('Displacementx(\tau)');title('Displacementasafunctionof(\tau)');axis([040-2.02.0]);text(2.7,-1.3,'阻尼系数=0.1');text(3.6,-0.1,'1.0');text(3.6,1.0,'5.0');subplot(2,1,2);xlabel('Displacement');ylabel('Velocity');title('Phaseportrait');axis([-2.02.0-2.02.0]);text(0.7,-1.25,'阻尼系数=0.1');text(0.8,-0.65,'1.0');text(0.8,0.1,'5.0');程序运行结果7.2.3线性系统的强迫振动一.运动微分方程二.MATLAB求解若对应的函数文件ForceOcillation.mfunctionxdot=ForceOcillation(t,x,dummy,zeta,Omega,x0)xdot=[x(2);-2.0*zeta*x(2)-x(1)+x0*cos(Omega*t)];为了获得频谱图,建立函数文件AmplitudeSpectrum.mf