数值分析在机械领域应用

《机械运动的数值仿真》研究学校：河北联合大学研究生院—机械工程学院班级：2011级研究生5班专业：机械设计及自动化第21组主讲：杨超学习的主要内容一、机械运动相关内容二、数值分析算法的研究三、常微分方程初值问题数值求解四、龙格-库塔算法五、龙格-库塔方法的实际应用六、报告总结一、机械运动相关内容1、机械运动的相关简单概念机械运动在物理学中，把一个物体相对于另一个物体位置的变化称作为机械运动，简称运动。参照物要判断一个物体是否在运动，必须选择另一个物体作为标准，这个作为标准的物体叫做参照物。对于同一个物体的运动，选择的参照物不同，得出的结论也有可能是不同的。运动和静止的相对性自然界中一切物体都在运动，因为地球本身在自转，所以绝对静止的物体是不存在的。通常所描述的物体的运动或静止都是相对于某一个参照物而言的。2、机械运动的前沿科学—导弹的制导•（1）导弹制导的一般原理

在大气层内飞行的导弹，可由改变空气动力获得控制，有翼导弹一般用改变空气动力的方法来改变控制力。在大气层中或大气层外飞行的导弹，都可以用改变推力的方法获得控制。无翼导弹主要是用改变推力的办法来改变控制力，因无翼导弹在稀薄大气层内飞行时，弹体产生的空气动力很小。二、数值分析算法的研究1、数值分析方法意义数学是一种工具，用于解决日常生活、工业工程上的相关问题。针对于数值分析中的数学方法，我们小组将主要内容概括分解，

将使用到的方法进行对比分析。2、我们要求数值算法的稳定性排除病态算法的数值稳定性与病态问题：若某算法受初始误差或运算过程中的舍入误差影响较小，则称为数值稳定。若微小的初始误差都会对最终结果产生极大的影响，则称之为病态问题。3、数值分析主要部分。1各类插值方法我们讲过拉格朗日插值、牛顿插值、分段插值、样条插值2函数逼近及拟合3数值积分、欧拉法解常微分方程、龙格-库塔法解常微分方程、方程组。【1】插值对于牛顿插值相对于拉氏插值增加一个节点，所有的插值基本多项式要重新取、重新算.2而牛顿插值，节点增加，次数增加，即高次插值函数计算量大，有剧烈震荡，数值稳定性较差（例如龙格现象）；分段插值在分段点上仅连续（即函数值相等），但是有尖点，不光滑（尖点导数不连续）；样条函数可以解决以上问题：使插值函数既是低次阶分段函数，又是光滑的函数。【2】理解逼近问题与拟合问题：逼近问题：函数f(x)在区间[a,b]具有一阶光滑度，求多项式p(x)是f(x)-p(x)在某衡量标准下最小的问题。拟合问题：从理论上讲y=f(x)是客观存在的，但在实际中，仅仅从一些离散的数据（xi,yi）(i=1,2…)是不可能求出f(x)的准确表达式，只能求出其近似表达式φ（x）。插值问题与逼近问题的特点和区别：相同点：它们都是求某点值的算法。不同点：

A，被插值函数是未知的，而被逼近函数是已知的。B，插值函数在节点处与被插值函数相等。而逼近函数的值只要满足很好的均匀逼近即可。C，求p(x)的方法不同。【3】Romberg（龙贝格）求积法和Gauss求积法的基本思想：(主要研究方法)复化求积公式精度较高，但需要事先确定步长，欠灵活性，在计算过程中将步长逐次减半得到一个新的序列，用此新序列逼近I的算法为Romberg求积法。对插值型求积公式，若能选取适当的xk.Ak使其具有2n+1阶代数精度，则称此类求积公式为Gauss型。【4】Runge-Kutta方法的基本思想：借助于Taylor级数法的思想，将yn+1=yn+hy’(ξ)中的

y’(ξ)（平均斜率）表示为f在若干点处值的线性组合，通过选择适当的系数使公式达到一定的阶。1、导弹轨迹的描述涉及到常微分方程初值问题数值求解问题2、《常微分方程初值问题数值求解》方法有多少？Euler方法、向后Euler方法、梯形方法、改进

Euler方法、龙格-库塔方法三、常微分方程初值问题数值求解3、经典的《常微分方程初值问题数值求解》方法是什么方法？Euler方法4、《常微分方程初值问题数值求解》的优缺点分析Euler方法计算简单但精度差；向后Euler方法与

Euler方法误差相似；梯形方法比Euler方法精度高但算法复杂、计算量很大；改进Euler方法结合了Euler方法和梯形法的优点；5、最好的方法是？龙格-库塔方法四、龙格-库塔算法

y(0)

=1y

y¢=

y

-

2x

(0

<

x

<1)1、龙格-库塔算法应用举例求(3阶R

-K公式)解初值问题(步长h

=0.1)

2

262

3=

f

(x

+h,

y

-hK

+2hK

).K3hh=

f

(x

,

y

),K1=

y

+

h(K

+4K

+K

),yn

n

1

2K2

=

f

xn

+

,

yn

+

K1

,n

nn

1n+1Clear[a,b,x,y]x[0]=0;y[0]=1;h=0.1;x[n_]:=n*h;f[u_,v_]:=v-2u/v;K1[n_]:=f[x[n-1],y[n-1]]K2[n_]:=f[x[n-1]+h/2,y[n-1]+h/2*K1[n]]K3[n_]:=f[x[n-1]+h,y[n-1]-h*K1[n]+2h*K2[n]];y[n_]:=y[n-1]+h/6*(K1[n]+4K2[n]+K3[n]);Table[{x[n],y[n]},{n,0,6}]//N;MatrixForm[%]01.0.11.095440.21.183220.31.264910.41.341650.51.414220.61.48326y[0]

->

1y[0.1]

->

1.09545y[0.2]

->

1.18322y[0.3]

->

1.26491y[0.4]

->

1.34164y[0.5]

->

1.41421y[0.6]

->

1.483240

1.0.1

1.097740.2

1.187570.3

1.271290.4

1.350130.5

1.424990.6

1.496570

1.0.1

1.095440.2

1.183220.3

1.264910.4

1.341650.5

1.414220.6

1.48326精确解改进Euler近似解3阶R-K近似解五、龙格-库塔方法的实际应用

如下图4.1所示，假设有一烟花火箭，其初始条件为零。将其放在地方然后点火，该烟花火箭的初始质量为，其中粉末燃料占。经过实验得知，燃料的持续时间为。燃料所产生的恒定推力为。这也说明燃料的消耗率恒定。空气产生的阻力和烟花火箭的速度的平方成正比：。这里，要求选择一种数值方法对其运动过程进行仿真并且其截断误差为或者更高。要求计算出该烟花火箭的最高高度，同时计算出从燃料消耗到该烟花火箭运动到最高点的时间延迟。该实例要求其截断误差要求大于或等于，这就使得较为简单的欧拉法，中点法不适合本例。龙格-库塔法以其优异的数值特性成为解决本问题的首选，且该数值算法很容易Mmathematic中实现。

很显然，该问题属于变质量的运动学问题，在该运动过程中，其前两秒是在驱动力和阻力的共同作用下加速上升的，而后的时间内，该烟花火箭是在空气的阻力下减速上升的，同时注意到空气的阻力和速度的平方成正比。为了对该运动过程进行数值仿真，那么必须建立相应的微分方程组。分析该运动过程可知，应该将该运动过程分为两部分：加速上升过程和减速上升过程。从而得到相应的微分方程组。加速上升过程：111000T

1000kv2dh

dv--

g 1

=

dt

120

-35t

120

-35t=

v

dt

v1

(0)

=

0,

h1

(0)

=

0减速上升过程：211000kv2

dh2=

v

dt

dv-

g

2

=

-

dt

120

-

35

·

2

h2

(0)

=

h1

(2)

v

(0)

=

v

(2)

2•

式中，为上升的高度，为上升过程的速度，为重力加速度。表示加速上升过程的最终高度，表示加速上升过程的最终速度。使用龙格-库塔法求解如上的微分方程组。该系统的数值仿真结果已在编程体现。如图是该烟花火箭的上升过程高度的数值仿真，如图是其上升过程速度的数值仿真。同时亦可以得到烟花火箭上升的最大高度和问题中所需的时间延迟：delayhmax

t=198.462

(m=

6.185(s)Mmathematic仿真图形1

2

3

4

5烟花火箭机械运动过程速度的数值仿真6204060801

2

3

4

5

6烟花火箭机械运动过程高度上升数值仿真50100150200编程过程：f[x_,y_]:=5200/(120-35t)-0.4x^2/(120-35t)-10;g[x_,y_]:=x;{x,y}={0,0};h=0.1;t=0.1;xx=Table[0,{i,1,70}];yy=Table[0,{i,1,70}];tt=Table[0,{i,1,70}];Do[a=f[x,y];xa=x+h

(a+f[x+h,y+h*a])/2;b=g[x,y];ya=y+h

(b+g[x+h,y+h*b])/2;Print[k,"

",t,"

",xa,"",ya];{t,x,y,xx[[k]],yy[[k]],tt[[k]]}={t+h,xa,ya,xa,ya,t+h},{k,1,20}]tt{0.6,1.1,1.6,2.1,2.6,3.1,3.6,4.1,4.6,5.1,5.6,6.1,6.6,7.1,7.6}xx{17.3174,37.9564,61.2759,85.0712,104.674,113.404,109.842,75.5189,8.73997,-59.2801,-96.8113,-111.485,-117.825,-121.341,-123.752}yy{0.125,8.90869,28.0119,58.7749,101.435,153.898,210.725,265.771,303.655,308.15,278.635,230.354,174.737,115.949,55.4041}0.51.01.52.020406020406080ListPlot[Table[{tt[[i]],xx[[i]]},{i,1,20}],PlotStyle®PointSize[0.02]]0.5

1.0

1.5

2.0ListPlot[Table[{tt[[i]],xx[[i]]},{i,1,20}],PlotJoined®True]802040600.51.01.52.0204060ListPlot[Table[{tt[[i]],yy[[i]]},{i,1,20}],PlotStyle•

PointSize[0.02]]800.5

1.0

1.5

2.0ListPlot[Table[{tt[[i]],yy[[i]]},{i,1,20}],PlotJoined®True]80f[x_,y_]:=-0.4x^2/(120-35*2)-10;g[x_,y_]:=x;{x,y}={85.3977,78.2347};h=0.1;t=2.1;Do[a=f[x,y];xa=x+h

(a+f[x+h,y+h*a])/2;b=g[x,y];ya=y+h

(b+g[x+h,y+h*b])/2;Pri