數(shù)值計算大作業(yè)——劉

1、閨演2行士孝課程設(shè)計課程名稱：數(shù)值分析設(shè)計題目：數(shù)值計算大作業(yè)學(xué)號：S315070064姓名：劉峰完成時間：2015年10月25日題目一、非線性方程求根1 .題目假設(shè)人口隨時間和當(dāng)時人口數(shù)目成比例連續(xù)增長，在此假設(shè)下人口在短期內(nèi)的增長建立數(shù)學(xué)模型。(1)如果令N(t)表示在t時刻的人口數(shù)目，P表示固定的人口出生率，則人口數(shù)目滿足微分方程"9=PN(t),此方程的解為N(t)=N0e°；dt(2)如果允許移民移入且速率為恒定的v,則微分方程變成處史=PN(t)+v,dt此方程的解為N(t)=Noe+-v-(e-1);假設(shè)某地區(qū)初始有1000000人，在第一年有435000人移

2、入，又假設(shè)在第一年年底該地區(qū)人口數(shù)量1564000人，試通過下面的方程確定人口出生率P,精確到10”；且通過這個數(shù)值來預(yù)測第二年年末的人口數(shù)，假設(shè)移民速度v保持不變。:4350001564000=1000000e：(e-1)2 .數(shù)學(xué)原理采用牛頓迭代法，牛頓迭代法的數(shù)學(xué)原理是，對于方程f(x)=0,如果f(x)是線性函數(shù)，則它的求根是很容易的，牛頓迭代法實質(zhì)上是一種線性化方法，其基本思想是將非線性方程f(x)=0逐步歸結(jié)為某種線性方程來求解。設(shè)已知方程f(x)=0有近似根xk(假定f'(x)±0),將函數(shù)f(x)在點(diǎn)Xk進(jìn)行泰勒展開，有f(x):f(xk)f(xk)(x-xk

3、).于是方程f(x)=0可近似地表示為f(x)f(xk)(x-xk)=0這是個線性方程，記其根為xk+,則xk+的計算公式為xk1=xkf(xk)f(xk)k=0,1,2,；這就是牛頓迭代法，簡稱牛頓法3 .程序設(shè)計作出函數(shù)的圖像，大概估計出根的位置fplot('1000*exp(x)+(435*x)*(exp(x)-1)-1564',03);grid大概估計出初始值x=0.5functionp1,err,k,y=newton(f,df,p0,delta,max1)% f是非線性系數(shù)% df是f的微商% p0是初始值% dalta是給定允許誤差% max1是迭代的最大次數(shù)% p

4、1是牛頓法求得的方程近似解% err是p0誤差估計% k是迭代次數(shù)p0,feval('f,p0)fork=1:max1p1=p0-feval('f,p0)/feval('df,p0);err=abs(p1-p0);p0=p1;p1,err,k,y=feval('f,p1)if(err<delta)|(y=0),break,endp1,err,k,y=feval('f,p1)endfunctiony=f(x)y=1000000*exp(x)+435000*(exp(x)-1)/x-1564000;functiony=df(x)y=1000000*ex

5、p(x)+435000*(exp(x)/x-(exp(x)-1)/xA2);4 .結(jié)果分析與討論在MATLAB的commandwindoW俞入newton('f,'df,1.2,10A(-4),10)運(yùn)行后得出結(jié)果p0=0.5000p1=0.1679err=0.3321k=1y=9.2415e+004p1=0.1031err=0.0648k=2y=2.7701e+003p1=0.1010err=0.0021k=3y=2.6953p1=0.1010err=2.0129e-006k=4y=2.5576e-006ans=0.1010運(yùn)算后的結(jié)果為P=0.1010,通過這個數(shù)值來預(yù)測第

6、二年年末的人口數(shù)，Kt)=1000000e°.1010t435222(/1010t7)0.1010t=2時候?qū)τ趂(x)=2187945.865實踐表明，當(dāng)初始值難以確定時，迭代法就不一定收斂了，因此要根據(jù)問題實際背景或者二分法先得一個較好的初始值，然后再進(jìn)行迭代；再者迭代函數(shù)選擇不合適的話，采用不動點(diǎn)迭代法也有可能出現(xiàn)不收斂的情況；因此我采用的是牛頓法。題目二：線性方程組求解1 .題目假設(shè)一個物體可以位于n+1個等距點(diǎn)X0,X1，Xn的任意位置，當(dāng)物體在Xi位置時，它只能等可能的移動到Xi或者Xi+1,而不能直接移動到其他任何位置，概率Pi表示物體從位置Xi開始在到達(dá)右端點(diǎn)Xn之前

7、到達(dá)左端點(diǎn)X0的概率，顯然11P0=1,Pn=0,且有P=2Pi-1+-Pi+1,i=1,2；,n-1既有下面方程組：11一2一P11P2-ll取n=10對方程組進(jìn)行求解（迭代法或者直接法）2 .數(shù)學(xué)原理在解微分方程的邊值問題、熱傳導(dǎo)方程以及船體數(shù)學(xué)放樣中建立的三次樣條函數(shù)等工程技術(shù)問題時，經(jīng)常遇到下面形式的線性方程組："biC1X2d2a2b2an4bnaCnbXnAXndn凡I方程簡記Ax=d，該線性方程稱為三對角線方程組，其系數(shù)矩陣A滿足條bAG>0,h|*ai十"=0,i=2,n-1bn>Cn>0所以為弱對角陣可以采用追趕法進(jìn)行計算，利用三對角矩陣

8、的LU分解建立計算量更少的線性方程組求解公式。將系數(shù)矩陣A進(jìn)行克勞特分解，即A分解為下三角矩陣和單位上三角矩陣的乘積;A=bia2Cib2riPiianJbnicnianbn其中：iPi工，為待定系數(shù),直接利用矩陣乘法公式可得b=iCi=：i-ib=iLi=2,3,n,G=；-i=2,3,n-i,于是推得計算：Bi,*的公式=bi%Pi，i=2,3，n,;Pi=G/cti,i=2,3，n;求解由此計算出L和U中的全部元素，完成了系數(shù)矩陣A的克勞特分解線性方程組Ax=d等價于求解Ly=dftUx=y因而得到解三對角線性方程組的追趕法公式計算bi的遞推公式：bi=q%。=gbi-ah,i=2,3,

9、n-i(2)解Ly=dyiV=di-aiy_b-aibi4,i=2,3,nXn=yn,我們將計算系數(shù)1，一：解Ux=yXi=V-bixi,i=n-i,i2TT二和yTy2TT丫門二稱為追的過程,將計算方程組的解XnTXnTTXi稱為趕的過程。整個過程為追趕法的思想。3 .程序設(shè)計functionx=chase(a,b,c,f)%求解線性方程組Ax=f,其中A是三對角陣%a是矩陣A的下對角線元素a(1)=0%b是矩陣A的對角線元素%c是矩陣A的上對角線元素c(N)=0%f是方程組的右端向量n=length(b);ifn-1=length(a)fori=n-1:-1:1a(i+1)=a(i);en

10、dendc(1)=c(1)/b(1);f(1)=f/b(1);fori=2:n-1b(i)=b(i)-a(i)*c(i-1);c(i)=c(i)/b(i);f(i)=(f(i)-a(i)*f(i-1)/b(i);endf(n)=(f(n)-a(n)*f(n-1)/(b(n)-a(n)*c(n-1);fori=n-1:-1:1f(i)=f(i)-c(i)*f(i+1);endx=f;4 .結(jié)果分析與討論A的系數(shù)矩陣為A=1,-0.5,0,0,0,0,0,0,0,0;-0.5,1,-0.5,0,0,0,0,0,0,0;0,-0.5,1,-0.5,0,0,0,0,0,0;0,0,-0.5,1,-0.

11、5,0,0,0,0,0;.0,0,0,-0.5,1,-0.5,0,0,0,0;0,0,0,0,-0.5,1,-0.5,0,0,0;0,0,0,0,0,-0.5,1,-0.5,0,0;0,0,0,0,0,0,-0.5,1,-0.5,0;.0,0,0,0,0,0,0,-0.5,1,-0.5;0,0,0,0,0,0,0,0,-0.5,1;所以在MATLABt令窗口輸入>>a=-0.5,-0.5,-0.5,-0.5,-0.5,-0.5,-0.5,-0.5,-0.5,0>>b=1,1,1,1,1,1,1,1,1,1>>c=-0.5,-0.5,-0.5,-0.5,-0.

12、5,-0.5,-0.5,-0.5,-0.5,0>>f=0,5,0,0,0,0,0,0,0,0,0得到此題中的a,b,c,f矩陣:a=-0.5000-0.5000-0.5000-0.5000-0.5000-0.5000-0.5000-0.5000-0.5000b=1111111111c=-0.5000-0.5000-0.5000-0.5000-0.5000-0.5000-0.5000-0.5000-0.50000f=0.5000000000000然后在MATLA呻調(diào)用之前保存的迭代法函數(shù)function,在命令窗口中輸入：chase(a,b,c,f)回車得到結(jié)果：>>x=

13、chase(a,b,c,f)x=0.90000.80000.70000.60000.50000.40000.30000.20000.10000追趕法為一種特殊的LU分解法。追趕法是求解三對角矩陣的常用方法，但從整體編程角度分析，其程序編寫較迭代法復(fù)雜，但通用性較好。追趕法求解三對角矩陣不但節(jié)省存儲單元，而且可以減少計算量，是工程技術(shù)中比較常用的數(shù)學(xué)工具。三、數(shù)值積分1、題目衛(wèi)星軌道是一個橢圓，橢圓周長的計算公式是S=4ajjl-(f)2sin2ede,a這里a是橢圓的半長軸,c是地球中心與軌道中心(橢圓中心)的距離，記h為近地點(diǎn)距離，H為遠(yuǎn)地點(diǎn)距離，R=6371公里為地球半徑，則a/R+H+h

14、Hh,某人造衛(wèi)星近地點(diǎn)距離h=536公里，遠(yuǎn)地點(diǎn)距離22H=2483公里，試用Romberg方法求衛(wèi)星軌道的周長，精確到10”。2.數(shù)學(xué)原理龍貝格方法是在梯形公式、辛普森公式和柯特斯公式之間的關(guān)系的基礎(chǔ)上，構(gòu)造出一種加速計算積分的方法。作為一種外推算法，它在不增加計算量的前提下提高了誤差的精度。龍貝格方法的主要過程是將粗糙的梯形公式Tn(f)逐步加工成精度較高的辛普森公式S(f)和科特斯公式Cn(f)的方法稱為龍貝格方法。復(fù)化梯形公式n1n1Tn(f)='1f(Xk)f(Xk.1)l-“；f(Xi)f(Xi1)k=02i=02在復(fù)化梯形公式中，每個內(nèi)節(jié)點(diǎn)。,''&quo

15、t;,既是前一個小區(qū)間的終點(diǎn)，又是后一個小區(qū)間的起點(diǎn)，因此上式可以改寫為hnTn(f)-f(a)2-f(Xk)f(b)2.kd復(fù)化梯形公式余項E="i2ah2f”()a,b復(fù)化梯形公式的遞推公式為b-aTi=<f(a)f(b)2-8-復(fù)化辛普森求積公式_1_hT2nTnf(Xi)2與復(fù)化梯形公式Tn(f)類似，每個內(nèi)節(jié)點(diǎn)Xl，X2；Xn需用兩次，因此有*f)=hgRR)f顯然復(fù)化辛普森公式在n趨于無窮大時，他的收斂速度比復(fù)化梯形公式更快。以To(k)表示二分k次后求得的梯形值，且以Tk)表示序列飛的的m次加速度，理查森外推法的遞推公式可寫成m(k)-.4T(k1)1(k)1m_

16、mJmm.TmJ,k-l,2,4-14-1龍貝格算法的計算過程如下:(1)取k=0,h=b-a,求T0(0)=-If(a)+f(b)2(2)利用變步長梯形公式To(k),其中k為區(qū)間的二分次數(shù)，即2k1-._1f(Xi)j2To(k)(f)=1To(k)(f)22j王2k-1嚶一記小)fa(2j1)22jzo(3)依橫行次序求加速值，逐個求出的第k行其余各元素Tj(J)(j=12，k)(4)當(dāng)相鄰對角元素之差的絕對值小于預(yù)先給定的精度時，終止計算。表3-1龍貝格算法遞推表khT0(k)T1(k)T2(k)T3(k)T4(k)一.0b-aT0(0)1b-a2T0(1)T1(0)2b-a4T。T1

17、T2(0)3b-a8T工T2(1)丁T34ba16T0(4)T1T2T3T4(0)aa-aa:+3.程序設(shè)計functionR=romberg(f,a,b,n)formatlongR=zeros(n+1,n+1);R(0+1,0+1)=(b-a)/2*(feval(f,a)+feval(f,b);fori=1:n,h=(b-a)/2Ai;s=0;fork=1:2A(i-1),s=s+feval(f,a+(2*k-1)*h);endR(i+1,0+1)=R(i-1+1,0+1)/2+h*s;endforj=1:n,fac=1/(4Aj-1);form=j:n,R(m+1,j+1)=R(m+1,j-1+1)+fac*(R(m+1,j-1+1)-R(m-1+1,j-1+1);endend-10-4.結(jié)果分析與討論本題根據(jù)算法原理在matlab中編寫完龍貝格算法的自定義程序后，直接輸入符合格