机械机构计算机辅助分析

1、机械机构计算机辅助分析专业：农业生物环境与能源工程姓名：迟明路v二分法原理 通过每次把f(x)的零点所在小区间收缩一半的方法，使区间的两个端点逐步迫近函数的零点，以求得零点的近似值，这种方法叫做二分法。 二分法也称对分法或区间套法，是逐次逼近一个实根的简单方法。 02xx21.用二分法求解方程二分法程序框图 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(f

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

3、牛顿法求解方程0e-ex牛顿法程序框图 Matlab函数M文件如下：vfunction x=newton(fx,dfx,x0,eps,N)vx=x0; f0=eval(fx); df0=eval(dfx);vn=0;v disp( n xn xn+1 fn+1 );vwhile n=Nv x1=x0-f0/df0;v x=x1;f1=eval(fx);v X=n,x0,x1,f1;v disp(X);v if abs(x0-x1)epsv fprintf(The procedure was successful.)v returnv else v n=n+1;v x0=x1;f0=f1;v e

4、ndvendvif n=N+1v fprintf(the method failed after N iterations. )v endvfx=exp(1)-exp(x);vdfx=-exp(x);vnewton(fx,dfx,0.5,10(-5),100) vThe procedure was successful.vans = 1.0000 梯形法原理 给定步长，确定梯形的高之后，求出所有梯形的面积之和，再根据定积分的几何意义可知曲顶高梯形的面积就是定积分的值。复化梯形公式为：niiixfxfhfT112)()()(nabh/ )( 其中,n是区间的等分数。ba,3、用梯形法求解10si

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

6、求：微分方程的解y(0)=1,h=0.2,区间0,1。yx2y y欧拉法程序框图 Matlab函数M文件如下：vfunction x,y=euler(fun,x0,y0,xfine,h) v N=(xfine-x0)/h;v x(1)=x0;y(1)=y0;v for n=1:Nv x(n+1)=x(n)+h;v y(n+1)=y(n)+h*feval(fun,x(n),y(n);vendvfun=inline(y-2*x/y,x,y);vx,y=euler(fun,0,1,1,0.2)vx = 0 0.2000 0.4000 0.6000 0.8000 1.0000vy =1.0000 1.

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

8、(fun,0,1,1,0.2)vx = 0 0.2000 0.4000 0.6000 0.8000 1.0000vy = 1.0000 1.1867 1.3483 1.4937 1.6279 1.7542v黄金分割法程序框图输入区间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);yynv 5、用黄金分割法求解方程f(x

9、)=3x4-16x3+30 x2-24x+8在区间-10，10，=0.01 v Matlab函数M文件如下：vfunction y=gold(f,a,b,eps)v while abs(b-a)epsv x1=a+0.618*(b-a);v x2=a+0.382*(b-a);v f1=3*x14-16*x13+30*x12-24*x1+8 ;v f2=3*x24-16*x23+30*x22-24*x2+8 ;v if f1f2v b=x1;x1=x2;v x2=a+0.382*(b-a);v elsev a=x2;x2=x1;v x1=a+0.618*(b-a); v endv endvx=(

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

11、1f2v b=x1;x1=x2;v x2=a+F(i-2)/F(i)*(b-a);v elsev a=x2;x2=x1;v x1=a+F(i-1)/F(i)*(b-a);v endv endv endv x=(x1+x2)/2vfibonacci(0.2,1,0.001,50)vx =v 0.99897、用进退法求f(x)= 3x4-16x3+30 x2-24x+8 ，a0=0.5，h=0.5，=0.001 Matlab函数M文件如下：vfunction y=jintui(a0,h,eps)v x1=a0;v x2=x1+h; v i=1;j=1; vwhile(abs(x2-x1)eps)v

12、 f1=3*x14-16*x13+30*x12-24*x1+8;v f2=3*x24-16*x23+30*x22-24*x2+8;v if(f1f2)v x1=x2;v x2=x1+2i*h;v i=i+1;v elsev x2=x1;v x1=x2-h/4j;v j=j+1;v endv endvx=(x1+x2)/2vf= 3*x4-16*x3+30*x2-24*x+8v运行结果如下：vjintui(0.5,0.5,0.001)vx =v 1.9998vf =v 3.5751e-007v8、质点在振动分离筛上的分析v(1)推导出k1,k2,k3。v(2)分析=-25，-10，0，10，25

13、，推导出=030的变化范围。左移右移跃起的k的变化范围。绘制=030时，k的变化规律。 v（1）筛面上任一点位移： R为OA的长度。v速度方程：v加速度方程：v物料在筛面上受到的牵连惯性力为： 当曲柄在左半轴周（2,3象限)时，其惯性力的方向指向x的负方向，摩擦力的方向指向D。当物料向下移动且没有跃起时，则： 为物料相对筛面的加速度。 整理得：v由于-1 coswt 0,可得：0tansin)cos(cos0cos)sin(cosF22NNrFamgawtmRwmaamgawtmRwwtRwxcos2 wtRwxsinwtRw cosS2wtRcosxra0)sin(g)cos(cos)sin

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

15、ace(0.0001,30* *pi/180,100); pi/180,100); vd d=30=30* *(pi/180); (pi/180); vk1=sin(k1=sin(d d- -b b)./cos()./cos(c c- -b b+ +d d); ); vk2=sin(k2=sin(d d+ +b b)./cos()./cos(c c- -b b- -d d); );vk3=cos(k3=cos(b b)./abs(sin()./abs(sin(c c- -b b););0FNa322)sin(cos0cos1cos)sin(cosRkaagRwwtagawtwc1=c*180/

16、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); v9（1）对于空间Vx,Vy,Vz的表达式。v（2）联合收割谷物清选（小石头，饱满谷物，不成熟谷物，短秸杆）.其中V0=5,=10m/s,=10清凉室高度是0.3m,u=0-30m/s，速度分别是25 m/s， 16 m/s ，14 m/s的位移曲线。v（1）气流的速度记为u，VL记为质点的飘浮速度，则质点相对于气流的相对速度v的表达式为：v转换成三维微分方程：uVv

17、FgmdtVdmgVuVVVVgdtdVVuVVVVgdtdVuVuVVVVgdtdVzxzyLzyxzyLyxxzyLx222222222222)()()()(v质点初始速度v0在x、y、z方向上的分量：v将时间分为若干小等分，每等分为t，求出t时间后的速度，再依据已知速度求出下一个t后的速度，以此类推，可以求所有点的速度。v 令：sinsincoscoscos000000vvvvvvzyxgVuVVVVgVVVfVuVVVVgVVVfuVuVVVVgVVVfzxzyLzyxyxzyLzyxxxzyLzyx222232222222221)(),()(),()()(),(),(),(2),()

18、,(2),(),(2t111111111111332211iiiiiiiiiiiiiiiiiiiiiiiizyxzyxzzzyxzyxyyzyxzyxxxVVVfVVVftVVVVVfVVVftVVVVVfVVVfVVv在微小的t时间内,其速度可看作初始速度和终结速度的平均值,物料的位移量为：v质点的位移：v（2）气流分别为25m/s，16m/s，14m/s的位移曲线iiiiiiiiizzzyyyxxx1112/)(2/)(2/)(111tVVztVVytVVxiiiiiizziyyixxifor j=1:1:4 u=0; if j=1 xL=24; elseif j=2 xL=10; els

19、eif 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=-v0*sin(a); u=u+h; x1=0; y1=0; while x1=0.3 & y1=1 f1=-g/xL2*sqrt(vx+u*sin(b)2+(vy-u*cos(b)2)*(vx+u*sin(b)+g; f2=-g/xL2*sqrt(vx+u*sin(b)2+(vy-u*cos(b)2)*(vy-u

20、*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/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

21、; elseif j=2 v2(i)=y2; elseif j=3 v3(i)=y2; else v4(i)=y2; end endendwhile 1 fengsu=input(请输入风速请输入风速:); if fengsu0 break; else num=input(请输入输出参数个数请输入输出参数个数:); n=round(fengsu/h); if n=0n=1; end out(1)=v1(n); out(2)=v2(n); out(3)=v3(n); out(4)=v4(n); if num=1 location=out(1) end if num=2 location=out(

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