數(shù)值分析在機(jī)械領(lǐng)域應(yīng)用

《機(jī)械運(yùn)動(dòng)的數(shù)值仿真》研究學(xué)校：河北聯(lián)合大學(xué)研究生院—機(jī)械工程學(xué)院班級(jí)：2011級(jí)研究生5班專業(yè)：機(jī)械設(shè)計(jì)及自動(dòng)化第21組主講：楊超學(xué)習(xí)的主要內(nèi)容一、機(jī)械運(yùn)動(dòng)相關(guān)內(nèi)容二、數(shù)值分析算法的研究三、常微分方程初值問題數(shù)值求解四、龍格-庫塔算法五、龍格-庫塔方法的實(shí)際應(yīng)用六、報(bào)告總結(jié)一、機(jī)械運(yùn)動(dòng)相關(guān)內(nèi)容1、機(jī)械運(yùn)動(dòng)的相關(guān)簡(jiǎn)單概念機(jī)械運(yùn)動(dòng)在物理學(xué)中，把一個(gè)物體相對(duì)于另一個(gè)物體位置的變化稱作為機(jī)械運(yùn)動(dòng)，簡(jiǎn)稱運(yùn)動(dòng)。參照物要判斷一個(gè)物體是否在運(yùn)動(dòng)，必須選擇另一個(gè)物體作為標(biāo)準(zhǔn)，這個(gè)作為標(biāo)準(zhǔn)的物體叫做參照物。對(duì)于同一個(gè)物體的運(yùn)動(dòng)，選擇的參照物不同，得出的結(jié)論也有可能是不同的。運(yùn)動(dòng)和靜止的相對(duì)性自然界中一切物體都在運(yùn)動(dòng)，因?yàn)榈厍虮旧碓谧赞D(zhuǎn)，所以絕對(duì)靜止的物體是不存在的。通常所描述的物體的運(yùn)動(dòng)或靜止都是相對(duì)于某一個(gè)參照物而言的。2、機(jī)械運(yùn)動(dòng)的前沿科學(xué)—導(dǎo)彈的制導(dǎo)?（1）導(dǎo)彈制導(dǎo)的一般原理

在大氣層內(nèi)飛行的導(dǎo)彈，可由改變空氣動(dòng)力獲得控制，有翼導(dǎo)彈一般用改變空氣動(dòng)力的方法來改變控制力。在大氣層中或大氣層外飛行的導(dǎo)彈，都可以用改變推力的方法獲得控制。無翼導(dǎo)彈主要是用改變推力的辦法來改變控制力，因無翼導(dǎo)彈在稀薄大氣層內(nèi)飛行時(shí)，彈體產(chǎn)生的空氣動(dòng)力很小。二、數(shù)值分析算法的研究1、數(shù)值分析方法意義數(shù)學(xué)是一種工具，用于解決日常生活、工業(yè)工程上的相關(guān)問題。針對(duì)于數(shù)值分析中的數(shù)學(xué)方法，我們小組將主要內(nèi)容概括分解，

將使用到的方法進(jìn)行對(duì)比分析。2、我們要求數(shù)值算法的穩(wěn)定性排除病態(tài)算法的數(shù)值穩(wěn)定性與病態(tài)問題：若某算法受初始誤差或運(yùn)算過程中的舍入誤差影響較小，則稱為數(shù)值穩(wěn)定。若微小的初始誤差都會(huì)對(duì)最終結(jié)果產(chǎn)生極大的影響，則稱之為病態(tài)問題。3、數(shù)值分析主要部分。1各類插值方法我們講過拉格朗日插值、牛頓插值、分段插值、樣條插值2函數(shù)逼近及擬合3數(shù)值積分、歐拉法解常微分方程、龍格-庫塔法解常微分方程、方程組?！?】插值對(duì)于牛頓插值相對(duì)于拉氏插值增加一個(gè)節(jié)點(diǎn)，所有的插值基本多項(xiàng)式要重新取、重新算.2而牛頓插值，節(jié)點(diǎn)增加，次數(shù)增加，即高次插值函數(shù)計(jì)算量大，有劇烈震蕩，數(shù)值穩(wěn)定性較差（例如龍格現(xiàn)象）；分段插值在分段點(diǎn)上僅連續(xù)（即函數(shù)值相等），但是有尖點(diǎn)，不光滑（尖點(diǎn)導(dǎo)數(shù)不連續(xù)）；樣條函數(shù)可以解決以上問題：使插值函數(shù)既是低次階分段函數(shù)，又是光滑的函數(shù)?！?】理解逼近問題與擬合問題：逼近問題：函數(shù)f(x)在區(qū)間[a,b]具有一階光滑度，求多項(xiàng)式p(x)是f(x)-p(x)在某衡量標(biāo)準(zhǔn)下最小的問題。擬合問題：從理論上講y=f(x)是客觀存在的，但在實(shí)際中，僅僅從一些離散的數(shù)據(jù)（xi,yi）(i=1,2…)是不可能求出f(x)的準(zhǔn)確表達(dá)式，只能求出其近似表達(dá)式φ（x）。插值問題與逼近問題的特點(diǎn)和區(qū)別：相同點(diǎn)：它們都是求某點(diǎn)值的算法。不同點(diǎn)：

A，被插值函數(shù)是未知的，而被逼近函數(shù)是已知的。B，插值函數(shù)在節(jié)點(diǎn)處與被插值函數(shù)相等。而逼近函數(shù)的值只要滿足很好的均勻逼近即可。C，求p(x)的方法不同。【3】Romberg（龍貝格）求積法和Gauss求積法的基本思想：(主要研究方法)復(fù)化求積公式精度較高，但需要事先確定步長，欠靈活性，在計(jì)算過程中將步長逐次減半得到一個(gè)新的序列，用此新序列逼近I的算法為Romberg求積法。對(duì)插值型求積公式，若能選取適當(dāng)?shù)膞k.Ak使其具有2n+1階代數(shù)精度，則稱此類求積公式為Gauss型?！?】Runge-Kutta方法的基本思想：借助于Taylor級(jí)數(shù)法的思想，將yn+1=yn+hy’(ξ)中的

y’(ξ)（平均斜率）表示為f在若干點(diǎn)處值的線性組合，通過選擇適當(dāng)?shù)南禂?shù)使公式達(dá)到一定的階。1、導(dǎo)彈軌跡的描述涉及到常微分方程初值問題數(shù)值求解問題2、《常微分方程初值問題數(shù)值求解》方法有多少？Euler方法、向后Euler方法、梯形方法、改進(jìn)

Euler方法、龍格-庫塔方法三、常微分方程初值問題數(shù)值求解3、經(jīng)典的《常微分方程初值問題數(shù)值求解》方法是什么方法？Euler方法4、《常微分方程初值問題數(shù)值求解》的優(yōu)缺點(diǎn)分析Euler方法計(jì)算簡(jiǎn)單但精度差；向后Euler方法與

Euler方法誤差相似；梯形方法比Euler方法精度高但算法復(fù)雜、計(jì)算量很大；改進(jìn)Euler方法結(jié)合了Euler方法和梯形法的優(yōu)點(diǎn)；5、最好的方法是？龍格-庫塔方法四、龍格-庫塔算法

y(0)

=1y

y￠=

y

-

2x

(0

<

x

<1)1、龍格-庫塔算法應(yīng)用舉例求(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精確解改進(jìn)Euler近似解3階R-K近似解五、龍格-庫塔方法的實(shí)際應(yīng)用

如下圖4.1所示，假設(shè)有一煙花火箭，其初始條件為零。將其放在地方然后點(diǎn)火，該煙花火箭的初始質(zhì)量為，其中粉末燃料占。經(jīng)過實(shí)驗(yàn)得知，燃料的持續(xù)時(shí)間為。燃料所產(chǎn)生的恒定推力為。這也說明燃料的消耗率恒定?？諝猱a(chǎn)生的阻力和煙花火箭的速度的平方成正比：。這里，要求選擇一種數(shù)值方法對(duì)其運(yùn)動(dòng)過程進(jìn)行仿真并且其截?cái)嗾`差為或者更高。要求計(jì)算出該煙花火箭的最高高度，同時(shí)計(jì)算出從燃料消耗到該煙花火箭運(yùn)動(dòng)到最高點(diǎn)的時(shí)間延遲。該實(shí)例要求其截?cái)嗾`差要求大于或等于，這就使得較為簡(jiǎn)單的歐拉法，中點(diǎn)法不適合本例。龍格-庫塔法以其優(yōu)異的數(shù)值特性成為解決本問題的首選，且該數(shù)值算法很容易Mmathematic中實(shí)現(xiàn)。

很顯然，該問題屬于變質(zhì)量的運(yùn)動(dòng)學(xué)問題，在該運(yùn)動(dòng)過程中，其前兩秒是在驅(qū)動(dòng)力和阻力的共同作用下加速上升的，而后的時(shí)間內(nèi)，該煙花火箭是在空氣的阻力下減速上升的，同時(shí)注意到空氣的阻力和速度的平方成正比。為了對(duì)該運(yùn)動(dòng)過程進(jìn)行數(shù)值仿真，那么必須建立相應(yīng)的微分方程組。分析該運(yùn)動(dòng)過程可知，應(yīng)該將該運(yùn)動(dòng)過程分為兩部分：加速上升過程和減速上升過程。從而得到相應(yīng)的微分方程組。加速上升過程：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?

式中，為上升的高度，為上升過程的速度，為重力加速度。表示加速上升過程的最終高度，表示加速上升過程的最終速度。使用龍格-庫塔法求解如上的微分方程組。該系統(tǒng)的數(shù)值仿真結(jié)果已在編程體現(xiàn)。如圖是該煙花火箭的上升過程高度的數(shù)值仿真，如圖是其上升過程速度的數(shù)值仿真。同時(shí)亦可以得到煙花火箭上升的最大高度和問題中所需的時(shí)間延遲：delayhmax

t=198.462

(m=

6.185(s)Mmathematic仿真圖形1

2

3

4

5煙花火箭機(jī)械運(yùn)動(dòng)過程速度的數(shù)值仿真6204060801

2

3

4

5

6煙花火箭機(jī)械運(yùn)動(dòng)過程高度上升數(shù)值仿真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