机械机构计算机辅助分析

1、 二分法原理 通过每次把f(x)的零点所在小区间收缩 一半的方法，使区间的两个端点逐步迫近函 数的零点，以求得零点的近似值，这种方法 叫做二分法。 二分法也称对分法或区间套法，是逐次 逼近一个实根的简单方法。 02xx 2 1.用二分法求解方程 二分法程序框图 Matlab函数M文件如下： vfunction x=erfen(fx,xa,xb,n,eps) vx=xa;fa=eval(fx); vx=xb;fb=eval(fx); v disp( n xa xb xc fc ); vfor i=1:n xc=(xa+xb)/2;x=xc;fc=eval(fx); v X=i,xa,xb,xc,

2、fc; vdisp(X), if fc*fa0 xb=xc; else xa=xc; end if (xb-xa)/2eps,break,end end fx=x2-x-2; erfen(fx,-1.5,0.5,20,10(-5) ans = -1.0000 fx=x2-x-2; erfen(fx,1.5,2.5,20,10(-5) ans = 2.0000 牛顿法原理 牛顿法是解非线性方程（组）的一种重要迭 代方法。 牛顿法的一般公式： )( / )f(n1nnnxfxxx 牛顿法也常称作切线法，其求单实根的收敛速 度较快，但求重实根的收敛速度较快。 2、用牛顿法求解方程 0e-e x 牛顿

3、法程序框图 Matlab函数M文件如下： function x=newton(fx,dfx,x0,eps,N) x=x0; f0=eval(fx); df0=eval(dfx); n=0; disp( n xn xn+1 fn+1 ); while n=N x1=x0-f0/df0; x=x1;f1=eval(fx); X=n,x0,x1,f1; disp(X); if abs(x0-x1)eps fprintf(The procedure was successful.) return else n=n+1; x0=x1;f0=f1; end end if n=N+1 fprintf(the

4、 method failed after N iterations. ) end fx=exp(1)-exp(x); dfx=-exp(x); newton(fx,dfx,0.5,10(-5),100) The procedure was successful. ans = 1.0000 梯形法原理 给定步长，确定梯形的高之后，求出所有 梯形的面积之和，再根据定积分的几何意义 可知曲顶高梯形的面积就是定积分的值。 复化梯形公式为： n i iixfxf hfT 1 1 2 )()( )( nabh/ )( 其中,n是区间的等分数。 ba, 3、用梯形法求解 1 0 sinxdx 梯形法程序框图

5、 Matlab函数M文件如下： vfunction y=T_quad(a,b,n) v h=(b-a)/n; v x1=a; v S=0; v for i=1:n v T=1/2*h*(sin(x1)+sin(x1+h); v S=S+T; v x1=x1+h; vend v S vT_quad(0,1,100) vS = 0.4597 欧拉法原理 欧拉（Euler）法（亦称切线法）是基于查 分法的最简单的显式单步方法。尽管其精度较低， 但算法简单，起源最早，也是其它较复杂数值解 法的一个很好的向导。 欧拉法公式： ),(1jjjjtyhfyy ) 1,.,2 , 1 , 0(nj 4、欧拉法

6、和改进求欧拉法求：微分方程 的解y(0)=1,h=0.2,区间0,1。 y x2 y y 欧拉法程序框图 Matlab函数M文件如下： function x,y=euler(fun,x0,y0,xfine,h) N=(xfine-x0)/h; x(1)=x0;y(1)=y0; for n=1:N x(n+1)=x(n)+h; y(n+1)=y(n)+h*feval(fun,x(n),y(n); end fun=inline(y-2*x/y,x,y); x,y=euler(fun,0,1,1,0.2) x = 0 0.2000 0.4000 0.6000 0.8000 1.0000 y =1.0

7、000 1.2000 1.3733 1.5315 1.6811 1.8269 改进欧拉法： 改进欧拉法程序框图 Matlab函数M文件如下： function x,y=euler_g(fun,x0,y0,xfine,h) N=(xfine-x0)/h; x(1)=x0;y(1)=y0; for n=1:N x(n+1)=x(n)+h; ym=y(n)+h*feval(fun,x(n),y(n); y(n+1)=y(n)+h/2*(feval(fun,x(n),y(n) +feval(fun,x(n+1),ym); end fun=inline(y-2*x/y,x,y); x,y=euler_g

8、(fun,0,1,1,0.2) x = 0 0.2000 0.4000 0.6000 0.8000 1.0000 y = 1.0000 1.1867 1.3483 1.4937 1.6279 1.7542 黄金分割法程序框图 输入区间a,b，精度 a1=b-0.618*(b-a); y1=f(a1); a2=a+0.618*(b-a); y2=f(a2); a2-a1y1 停a=0.5(a1+a2) b=a2; a2=a1; y2=y1; a1=b-0.618*(b-a) y1=f(a1); a=a1; a1=a2; y1=y2; a2=a+0.618*(b-a) y2=f(a2); y y

9、n 5、用黄金分割法求解方程f(x)=3x4-16x3+30 x2-24x+8在 区间-10，10，=0.01 Matlab函数M文件如下： function y=gold(f,a,b,eps) while abs(b-a)eps x1=a+0.618*(b-a); x2=a+0.382*(b-a); f1=3*x14-16*x13+30*x12-24*x1+8 ; f2=3*x24-16*x23+30*x22-24*x2+8 ; if f1f2 b=x1;x1=x2; x2=a+0.382*(b-a); else a=x2;x2=x1; x1=a+0.618*(b-a); end end x

10、=(x1+x2)/2 运行结果如下： f=3*x4-16*x3+30*x2-24*x+8; gold(f,-10,10,0.01) x = 1.9984 6、斐波那契法求：f(x)=x+20/x在0.2，1 上的极值，=0.001 function y=fibonacci(a,b,eps,n) F(1)=1;F(2)=1; for i=3:n F(i)=F(i-1)+F(i-2); x1=a+F(i-1)/F(i)*(b-a); x2=a+F(i-2)/F(i)*(b-a); while abs(x2-x1)eps f1=x1+20/x1; f2=x2+20/x2; if f1f2 b=x1;

11、x1=x2; x2=a+F(i-2)/F(i)*(b-a); else a=x2;x2=x1; x1=a+F(i-1)/F(i)*(b-a); end end end x=(x1+x2)/2 fibonacci(0.2,1,0.001,50) x = 0.9989 7、用进退法求f(x)= 3x4-16x3+30 x2-24x+8 ，a0=0.5，h=0.5，=0.001 Matlab函数M文件如下： function y=jintui(a0,h,eps) x1=a0; x2=x1+h; i=1;j=1; while(abs(x2-x1)eps) f1=3*x14-16*x13+30*x12-

12、24*x1+8; f2=3*x24-16*x23+30*x22-24*x2+8; if(f1f2) x1=x2; x2=x1+2i*h; i=i+1; else x2=x1; x1=x2-h/4j; j=j+1; end end x=(x1+x2)/2 f= 3*x4-16*x3+30*x2-24*x+8 运行结果如下： jintui(0.5,0.5,0.001) x = 1.9998 f = 3.5751e-007 8、质点在振动分离筛上的分析 (1)推导出k1,k2,k3。 (2)分析=-25，-10，0，10，25，推导出 =030的变化范围。左移右移跃起的k的变化范围。 绘制=030时

13、，k的变化规律。 （1）筛面上任一点位移： R为OA的长度。 速度方程： 加速度方程： 物料在筛面上受到的牵连惯性力为： 当曲柄在左半轴周（2,3象限)时，其惯性力的方向指向x的负方 向，摩擦力的方向指向D。当物料向下移动且没有跃起时，则： 为物料相对筛面的加速度。 整理得： 由于-1 coswt 0,可得： 0tansin)cos(cos 0cos)sin(cosF 2 2 N Nr FamgawtmRwma amgawtmRw wtRwxcos 2 wtRwxsin wtRw cosS 2 wtRcosx r a 0)sin(g)cos(cos )sin(coscos 2 2 aawtRw

14、 awtmRwamgF N 1 2 )cos( )sin( k a a g Rw 当曲柄在右半轴周（1,4象限)时，惯性力的方向指向x的正方向，摩 擦力的方向指向C。当物料向上移动且没有跃起时，则： 由于-1 coswt 0,可得： 物料跃起条件： 当 时，由 式可以看出， 只有coswt 时，只有coswt0时，才有可能 ， 即物料在曲柄转至右半周时跃起。 （2 2）clear; clear; clcclc; ; a=input(a=input(输入角度：输入角度：); ); b b=a=a* *pi/180pi/180; ; c c=linspace(0.0001,30=linspace(

15、0.0001,30* *pi/180,100); pi/180,100); d d=30=30* *(pi/180); (pi/180); k1=sin(k1=sin(d d- -b b)./cos()./cos(c c- -b b+ +d d); ); k2=sin(k2=sin(d d+ +b b)./cos()./cos(c c- -b b- -d d); ); k3=cos(k3=cos(b b)./abs(sin()./abs(sin(c c- -b b);); 0 FNa 3 2 2 )sin( cos 0cos1 cos)sin(cosR k a a g Rw wt agawtw

16、 c1=c*180/pi; plot(c1,k1,-b,c1,k2,- g,c1,k3,-r); xlabel(c); ylabel(k); title(k1,k2,k3曲线图曲线图); legend(k1 ,k2 ,k3 ); axis(0 30 0 8); 9（1）对于空间Vx,Vy,Vz的表达式。 （2）联合收割谷物清选（小石头，饱满谷物， 不成熟谷物，短秸杆）.其中 V0=5,=10m/s,=10清凉室高度是 0.3m,u=0-30m/s，速度分别是25 m/s， 16 m/s ，14 m/s的位移曲线。 （1）气流的速度记为u，VL记为质点的飘浮速度，则 质点相对于气流的相对速度v的

17、表达式为： 转换成三维微分方程： uVv Fgm dt Vd m gVuVVV V g dt dV VuVVV V g dt dV uVuVVV V g dt dV zxzy L z yxzy L y xxzy L x 222 2 222 2 222 2 )( )( )()( 质点初始速度v0在x、y、z方向上的分量： 将时间分为若干小等分，每等分为t，求出t时间后的速度，再依据 已知速度求出下一个t后的速度，以此类推，可以求所有点的速度。 令： sin sincos coscos 00 00 00 vv vv vv z y x gVuVVV V g VVVf VuVVV V g VVVf u

18、VuVVV V g VVVf zxzy L zyx yxzy L zyx xxzy L zyx 222 2 3 222 2 2 222 2 1 )(),( )(),( )()(),( ),(),( 2 ),(),( 2 ),(),( 2 t 1111 1111 1111 33 22 11 iiiiiiii iiiiiiii iiiiiiii zyxzyxzz zyxzyxyy zyxzyxxx VVVfVVVf t VV VVVfVVVf t VV VVVfVVVfVV 在微小的t时间内,其速度可看作初始速度和终结速度的平均值,物料的位移 量为： 质点的位移： （2）气流分别为25m/s，16

19、m/s，14m/s的位移曲线 iii iii iii zzz yyy xxx 1 1 1 2/)( 2/)( 2/)( 1 1 1 tVVz tVVy tVVx ii ii ii zzi yyi xxi for j=1:1:4 u=0; if j=1 xL=24; elseif j=2 xL=10; elseif j= 3 xL=6.5; else xL=3; end clear; clc; w=30; H=0.3; L=1; v0=5; a=10*pi/180; b=a; g=9.8; t=0.001; h=w/2000; for i=1:1:2000 vx=v0*cos(a); vy=-v

20、0*sin(a); u=u+h; x1=0; y1=0; while x1=0.3 f2=-g/xL2*sqrt(vx+u*sin(b)2+(vy-u*cos(b)2)*(vy-u*cos(b); vx1=vx+t*f1; vy1=vy+t*f2; f3=-g/xL2*sqrt(vx1+u*sin(b)2+(vy1-u*cos(b)2)*(vx1+u*sin(b)+g; f4=-g/xL2*sqrt(vx1+u*sin(b)2+(vy1-u*cos(b)2)*(vy1-u*cos(b); vx2=vx+t/2*(f1+f3); vy2=vy+t/2*(f2+f4); dx=(vx2+vx)*t

21、/2; dy=(vy2+vy)*t/2; vx=vx1; vy=vy1; x2=x1+dx; y2=y1+dy; x1=x2; y1=y2; end if j=1 v1(i)=y2; elseif j=2 v2(i)=y2; elseif j=3 v3(i)=y2; else v4(i)=y2; end end end while 1 fengsu=input(请输入风速请输入风速:); if fengsu0 break; else num=input(请输入输出参数个数请输入输出参数个数:); n=round(fengsu/h); if n=0 n=1; end out(1)=v1(n);

22、out(2)=v2(n); out(3)=v3(n); out(4)=v4(n); if num=1 location=out(1) end if num=2 location=out(1:2) end if num=3 location=out(1:3) end if num=4 location=out(1:4) end end end 运行结果 请输入风速:25 请输入输出参数个数:4 location = -0.0311 0.0786 0.3931 1.0029 请输入风速:16 请输入输出参数个数:4 location = -0.0423 0.0008 0.0884 1.0061 请输入输出参数个数:4 location = -0.0442 -0.0111 0.0520 1.0023 10、曲柄滑块机构运动学分析 clear; clc; w=input(角速度(转数):r=); R=input(OA杆长度:R=); L=input(AB杆长度:L=); H=input(高度