子程序ex713f 函数程序应另存成一个文件ex713f.m function zprime=ex713f(t,z) global vt vm zprime-=[0;0];%给出t0之前zprime初值 zprime(1)=-vt-vm*z(1)/sqrt(z(1)2+z(2)^2); zprime(2)=-vm*z(2)/sqrt(z(1)^2+z(2)^2): %上面两句可换成一个矩阵语句： zprime=-vt*[1;0]-vm*z/sqrt(Z(1)^2+z(2)^2)
子程序ex713f 函数程序应另存成一个文件ex713f.m function zprime=ex713f(t,z) global vt vm zprime=[0;0]; % 给出t0之前zprime初值 zprime(1)=-vt-vm*z(1)/sqrt(z(1)^2+z(2)^2); zprime(2)=-vm*z(2)/sqrt(z(1)^2+z(2)^2); %上面两句可换成一个矩阵语句： zprime=-vt*[1;0]-vm*z/sqrt(z(1)^2+z(2)^2);
程序en713运行结果 输入以下参数： t=500;Vm=1000 4000 [x0,y0]=[3000;4000] vt=300m/s 3000 tspan=[to,tfinal]=[0,4.5] 得到轨迹如右图。在给定tfinal2ooo vm=800m/s 时，必须使它小于遭遇点的 m=1000m/s 值，否则积分会进入死循环1000 而得不出结果。 将vm换成800，并相应地把 1000 0 100020003000 tfinal换成6，得到的轨迹位 于图中原轨迹的左上方。 图7-1-3导弹跟踪目标时的相对轨迹
程序exn713运行结果 输入以下参数： vt=500;vm=1000 [x0;y0]=[3000;4000] tspan=[t0,tfinal]=[0,4.5] 得到轨迹如右图。在给定tfinal 时，必须使它小于遭遇点的 值，否则积分会进入死循环 而得不出结果。 将vm换成800，并相应地把 tfinal换成6，得到的轨迹位 于图中原轨迹的左上方。 -1000 0 1000 2000 3000 0 1000 2000 3000 4000 vt=300m/s vm=1000m/s vm=800m/s 图7-1-3 导弹跟踪目标时的相对轨迹
例7-1-4 四连杆运动分析 如图，输入杆L1的转01=w1*t, (w1=100)求输出杆L3转角 3随时间的变化规律，并求 其角速度和角加速度。 解：◆建模： 四连杆的运动方程由X和Y方向的长度关系确定： L1c0s(Θ1)+L2c0s(θ2)-L3c0s(Θ3)-L0=0(1) L1sin(01)+L2sin(02)-L3sin(03)=0 (2) 在上述两个方程中消去2，便可消成一个方程。 给定日1，可求出满足此方程的3
例 7-1-4 四连杆运动分析 如图，输入杆L1的转θ1=w1*t, (w1=100)求输出杆L3转角 θ3随时间的变化规律，并求 其角速度和角加速度。 解：◆建模： 四连杆的运动方程由X和Y方向的长度关系确定： L1cos(θ1)+L2cos(θ2)-L3cos(θ3)-L0 = 0（1） L1sin(θ1)+L2sin(θ2)-L3sin(θ3) = 0 （2） 在上述两个方程中消去θ2，便可消成一个方程。 给定θ1，可求出满足此方程的θ3
计算和编程的思路 由(2) sin(Θ2)=(L3sin(Θ3)-L1sin(e1)/L2 将(1)式中的cos(2)代以sqrt(1-sin(2)2),得出方程： f=L1cos(01)+L2sqrt(1-(L3sin(03)-L1sin(01))2/L22)- L3c0s(03)-L0=0; 给定01，求能使f(03)=0的03值，求出3后，2就可由 sin(e2)=（L3sin(3)-L1sin(日1)/L2求得. 为了求能使f=0的3值，可调用MATLAB中的fzro函数。为此， 要把f=f(3)单独定义为一个MATLAB函数ex714f.m,在主 程序中调用它。为了把长度参数传给子程序，在主程序和 子程序中都加了全局变量语句(global),全局变量容易造成 程序的混乱，要特别小心，在复杂的程序中要尽量避免。 求得01,2和03后，就不难根据杆1的角速度求出杆3的角速 度
计算和编程的思路 由（2） sin(θ2)＝(L3sin(θ3) - L1sin(θ1))/L2 将(1)式中的cos(θ2)代以sqrt(1- sin(θ2)2),得出方程： f=L1cos(θ1)+L2sqrt(1-(L3sin(θ3)-L1sin(θ1))2/L22)- L3cos(θ3)-L0=0; 给定θ1，求能使ff(θ3)=0的θ3值,求出θ3后, θ2就可由 sin(θ2)= ( L3sin(θ3)- L1sin(θ1))/L2求得. 为了求能使f=0的θ3值,可调用MATLAB中的fzero函数。为此， 要把f =f(θ3)单独定义为一个MATLAB函数ex714f.m,在主 程序中调用它。为了把长度参数传给子程序，在主程序和 子程序中都加了全局变量语句(global),全局变量容易造成 程序的混乱，要特别小心，在复杂的程序中要尽量避免。 求得θ1,θ2和θ3后,就不难根据杆1的角速度求出杆3的角速 度
角速度计算的两种方法： 1.求瞬时速度，这是通常理论力学的解法，其依据就是杆2两端 点a和b的速度沿杆长方向的分量相等，通过三角关系，有 L1w1c0s(T/2-01+02)=L3w3c0s(θ3-T/2-02) 从而w3=L1w1c0s(T/2-61+62)/(L3c0s(Θ3-T/2-62) 由杆2两端点a和b的速度沿杆长垂直方向的分量之差，可以求 出杆2的角速度. w2=(-(L3sin(Θ3-T/2-02)-L1w1sin(π/2-61+62)/L2 2.求运动全过程的角位置，角速度，角加速度曲线，这只有借助 于计算工具才能做到，因为用手工算一个点就不胜其烦，算 几十个点是很难想象的.而由MATLAB编程调用fzero函数 时，要求给出一个近似猜测值，若连续算几十点前一个解就 可作为后一个解的猜测值，所以反而带来了方便 这样，本书将提供两个程序ex714a.m和ex714b.m来表述这两 种方法，它们所要调用的函数程序命名为ex714f.m
角速度计算的两种方法: 1. 求瞬时速度,这是通常理论力学的解法,其依据就是杆2两端 点a和b的速度沿杆长方向的分量相等,通过三角关系,有 L1w1cos(π/2-θ1+θ2)=L3w3cos(θ3-π/2-θ2) 从而 w3 = L1w1cos(π/2-θ1+θ2)/ (L3cos(θ3-π/2-θ2)) 由杆2两端点a和b的速度沿杆长垂直方向的分量之差,可以求 出杆2的角速度. w2 = (-(L3sin(θ3-π/2-θ2))- L1w1sin(π/2-θ1+θ2))/L2 2. 求运动全过程的角位置,角速度,角加速度曲线,这只有借助 于计算工具才能做到,因为用手工算一个点就不胜其烦, 算 几十个点是很难想象的.而由MATLAB编程调用fzero函数 时,要求给出一个近似猜测值,若连续算几十点,前一个解就 可作为后一个解的猜测值,所以反而带来了方便. 这样,本书将提供两个程序ex714a.m和ex714b.m来表述这两 种方法,它们所要调用的函数程序命名为ex714f.m