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

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

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

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

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

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

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

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

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

t=198.462

(m=

6.185(s)Mmathematic仿真圖形1

2

3

4

5煙花火箭機械運動過程速度的數(shù)值仿真6204060801

2

3

4

5

6煙花火箭機械運動過程高度上升數(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