數(shù)值分析課程設(shè)計

摘要實驗一線性方程組數(shù)值解法中，本文選取列主元Gauss消去法，并選取數(shù)據(jù)于我們課堂學(xué)習(xí)該章節(jié)時的課件中的例2進行實驗。列主元消去法是在Gauss消去法的根底上改良而得到的一種比擬快速和合理的求解線性方程組的方法。它的主要思路是通過對每次消元過程中主元的屢次選取以到達減小誤差和加快求解速度的一種消去法。使用列主元消去法相比于Gauss消去法根本上能控制舍入誤差的影響，并且選主元素較全主元消去法更為方便。實驗二插值法和數(shù)據(jù)擬合中，本文選取最小二乘擬合方法進行實驗，數(shù)據(jù)來源于我們課堂學(xué)習(xí)該章節(jié)時的課件中的多項式擬合例子進行實驗。最小二乘擬合是一種數(shù)學(xué)上的近似和優(yōu)化，利用的數(shù)據(jù)得出一條直線或者曲線，使之在坐標系上與數(shù)據(jù)之間的距離的平方和最小。利用excel的自帶函數(shù)可以較為方便的擬合線性的數(shù)據(jù)分析。實驗三數(shù)值積分中，本文選取復(fù)化Simpon積分方法進行實驗，通過將復(fù)化Simpson公式編譯成MATLAB語言求積分，完成實驗過程的同時，也對復(fù)化Simpon積分章節(jié)的知識進行了穩(wěn)固。實驗四常微分方程數(shù)值解，本文選取Runge-Kutta方法進行實驗，通過實驗了解Runge-Kutta法的收斂性與穩(wěn)定性同時學(xué)會了學(xué)會用Matlab編程實現(xiàn)Runge-Kutta法解常微分方程，并在實驗的過程中意識到盡管我們熟知的四種方法，事實上，在求解微分方程初值問題，四階法是單步長中最優(yōu)秀的方法，通常都是用該方法求解的實際問題，計算效果比擬理想的。最后，在數(shù)值方法實際應(yīng)用中，本文通過非線性擬合方法，并學(xué)會用這些函數(shù)解決實際人口增長模型問題，在1949年到1994年的數(shù)據(jù)根底上，通過構(gòu)建“線性增長模型”和“指數(shù)指數(shù)增長模型”兩個模型來對1999年的中國人口數(shù)據(jù)進行預(yù)測，最后通過中國人口年鑒中的數(shù)據(jù)對兩種方法進行準確性檢驗，并且比照得出中國人口數(shù)據(jù)在1949年到1999年期間更符合線性增長的結(jié)論。關(guān)鍵詞：列主元Gauss消去法最小二乘擬合復(fù)化Simpon積分Runge-Kutta目錄實驗一列主元Gauss消去法 11.1 實驗?zāi)康?11.2根本原理 11.3實驗內(nèi)容 31.4數(shù)據(jù)來源 41.5實驗結(jié)論 5實驗二最小二乘擬合 62.1實驗?zāi)康?62.2根本原理 62.3實驗內(nèi)容 72.4實驗結(jié)論 11試驗三復(fù)化Simpon積分 123.1實驗?zāi)康?123.2根本原理 123.3實驗內(nèi)容 123.4數(shù)據(jù)來源 133.5實驗結(jié)論 14實驗四Runge-Kutta方法 154.1實驗?zāi)康?154.2根本原理 154.3實驗內(nèi)容 164.4數(shù)據(jù)來源 164.5實驗結(jié)論 16實驗五數(shù)值方法實際應(yīng)用 175.1實驗?zāi)康?175.2根本原理 175.3實驗內(nèi)容 185.4數(shù)據(jù)來源 185.5 實驗結(jié)論 21參考文獻 23實驗一列主元Gauss消去法實驗?zāi)康睦斫獠⒄莆樟兄髟狦auss消去法的消去原理，熟悉算法流程；能有效利用列主元消去法解方程組；通過運用Matlab軟件編程實現(xiàn)列主元Gauss消去法；1.2根本原理設(shè)有線性方程組,其中，為非奇異矩陣,方程組的增廣矩陣為首先在的第1列選取絕對值最大的元素作為主元素，即選擇然后交換的第1行與第行〔交換后增廣矩陣為簡單起見仍記為，其元素仍記為〕。經(jīng)過第1次消元計算得到與原方程組等價的方程，其中，上述過程可記為重復(fù)上述計算過程，現(xiàn)假設(shè)已完成第步的選主元素過程，交換兩行并進行消元計算，此時約化為其中的元素仍記為，的元素仍記為.第步選主元素〔在右下角方陣的第1列內(nèi)選〕，即確定，使交換第行與行的元素，再進行消元計算，最后將原線性方程組化為回代可求解得算法描述：輸入數(shù)據(jù)和,置；對于，做如下步驟：按列選主元，即確定使(如果，那么為奇異矩陣，停止計算)；如果，那么交換第行與第行元素；消元計算：置。如果，那么停止計算〔此時〕;否那么進行回執(zhí)，對于置；置；輸出為線性方程組的解，為矩陣的行列式的值。1.3實驗內(nèi)容functionx=gauss(A,b)[n,n]=size(A);x=zeros(n,1);Aug=[A,b];fork=1:n-1[piv,r]=max(abs(Aug(k:n,k)));r=r+k-1;ifr>ktemp=Aug(k,:);Aug(k,:)=Aug(r,:);Aug(r,:)=temp;endifAug(k,k)==0,error('對角元出現(xiàn)0’),endAugforp=k+1:nAug(p,:)=Aug(p,:)-Aug(k,:)*Aug(p,k)/Aug(k,k);endendAugA=Aug(:,1:n);b=Aug(:,n+1);x(n)=b(n)/A(n,n);fork=n-1:-1:1x(k)=b(k);forp=n:-1:k+1x(k)=x(k)-A(k,p)*x(p);endx(k)=x(k)/A(k,k);end1.4數(shù)據(jù)來源為了驗證以上編寫的Matlab程序的正確性，試用列主元Gauss消去法解方程組解：只需輸入>>A=[1,2,-3;2,-1,3;3,-2,2];>>b=[1;5;1];>>gauss(A,b)程序運行得到如下結(jié)果Aug=12-312-1353-221Aug=3-2212-13512-31Aug=3.0000-2.00002.00001.000002.6667-3.66670.666700.33331.66674.3333Aug=3.0000-2.00002.00001.000002.6667-3.66670.6667002.12504.2500ans=132可見方程組的解為1.5實驗結(jié)論在順序Gauss消去法中，必須要求否那么將無法進行計算。這里需要指出的是：即使但其絕對值很小，那么由于舍入誤差的影響，也可能引起很大的誤差，而選取列主元就可以很好的解決這個問題，而且經(jīng)過實例驗證發(fā)現(xiàn)該程序也是正確的。實驗二最小二乘擬合2.1實驗?zāi)康膶W(xué)會并掌握數(shù)值最小二乘擬合的根本原理；自行學(xué)會用Matlab實現(xiàn)算法并將其應(yīng)用于實際問題；了解插值與擬合的異同點；2.2根本原理如果函數(shù)在假設(shè)干點處的值,便可根據(jù)插值原理來建立插值多項式作為的近似。但在科學(xué)實驗和生產(chǎn)實踐中，往往會遇到這樣一種情況，即節(jié)點上的函數(shù)值并不是很精確的，這些函數(shù)值是由實驗或觀測得到的數(shù)據(jù)，不可防止地帶有測量誤差，如果要求所得的近似函數(shù)曲線精確無誤地通過所有的點(xi,yi),就會使曲線保存著一切測試誤差。當(dāng)個別數(shù)據(jù)的誤差較大時,插值效果顯然是不理想的。此外,由實驗或觀測提供的數(shù)據(jù)個數(shù)往往很多,如果用插值法,勢必得到次數(shù)較高的插值多項式，這樣計算起來很煩瑣。換句話說:求一條曲線,使數(shù)據(jù)點均在離此曲線的上方或下方不遠處,所求的曲線稱為擬合曲線,它既能反映數(shù)據(jù)的總體分布,又不至于出現(xiàn)局部較大的波動,更能反映被逼近函數(shù)的特性,使求得的逼近函數(shù)與函數(shù)從總體上來說其偏差按某種方法度量到達最小,這就是曲線擬合。與函數(shù)插值問題不同,曲線擬合不要求曲線通過所有點,而是要求得到的近似函數(shù)能反映數(shù)據(jù)的根本關(guān)系。在某種意義上,曲線擬合更有實用價值。在對給出的實驗(或觀測)數(shù)據(jù)一般希望各實驗(或觀測)數(shù)據(jù)與擬合曲線的偏差的平方和最小,這就是最小二乘原理。函數(shù)插值是插值函數(shù)與被插函數(shù)在節(jié)處函數(shù)值相同,即而曲線擬合函數(shù)不要求嚴格地通過所有數(shù)據(jù)點,也就是說擬合函數(shù)在處的偏差(亦稱殘差〕不都嚴格地等于零。但是,為了使近似曲線能盡量反映所給數(shù)據(jù)點的變化趨勢,要求按某種度量標準最小。假設(shè)記向量,即要求向量某種范數(shù)最小,如的范數(shù)或范數(shù)即，或最小。為了方便于計算、分析與應(yīng)用，通常要求的范數(shù)即為最小。這種要求誤差〔偏差〕平方和最小的擬合稱為曲線擬合的最小二乘法。2.3實驗內(nèi)容給出一組數(shù)據(jù)點,利用最小二乘法擬合曲線，并估計誤差，作出擬合曲線。123456012345521123解：在MATLAB工作窗口中輸入程序>>x=[0,1,2,3,4,5];y=[5,2,1,1,2,3];plot(x,y,'r*'),legend('實驗數(shù)據(jù)(xi,yi)')xlabel('x'),ylabel('y'),title('數(shù)據(jù)點(xi,yi)的散點圖')可以看出這些點接近一條拋物線，因此設(shè)所求的多項式為在MATLAB中輸入編碼為：symsa1a2a0;fi=a2.*x.^2+a1.*x+a0編寫構(gòu)造誤差平方和的MATLAB程序fy=fi-y;J=sum(fy.^2)gJa0=diff(J,a0);Ja1=diff(J,a1);Ja2=diff(J,a2);Ja01=simple(Ja0)Ja11=simple(Ja1)Ja21=simple(Ja2)fi=[a0,a0+a1+a2,a0+2*a1+4*a2,a0+3*a1+9*a2,a0+4*a1+16*a2,a0+5*a1+25*a2]運行后屏幕顯示誤差平方和如下J=(a0-5)^2+(a0+2*a1+4*a2-1)^2+(a0+3*a1+9*a2-1)^2+(a0+4*a1+16*a2-2)^2+(a0+5*a1+25*a2-3)^2+(a0+a1+a2-2)^2運行后屏幕顯示關(guān)于和的線性方程組Ja01=12*a0+30*a1+110*a2-28Ja11=30*a0+110*a1+450*a2-60Ja21=110*a0+450*a1+1958*a2-244>>A=[12,30,110;30,110,450;110,450,1958];>>B=[28,60,244];>>C=B/A,f=poly2sym(C)C=4.7143-2.78570.5000f=(33*x^2)/7-(39*x)/14+1/2即，解一下方程組：解得：所求的多項式為：xi=[0,1,2,3,4,5];y=[5,2,1,1,2,3];n=length(xi)f=4.7143-2.7857*xi+0.5000*xi.^2x=0:1:5;F=4.7143-2.7857*xi+0.5000*xi.^2fy=abs(f-y);fy2=fy.^2;Ew=max(fy),E1=sum(fy)/n,E2=sqrt((sum(fy2))/n)plot(xi,y,'r*'),holdon,plot(x,F,'b-'),holdofflegend('數(shù)據(jù)點(xi,yi)','擬合曲線y=f(x)'),xlabel('x'),ylabel('y'),title('數(shù)據(jù)點(xi,yi)和擬合曲線y=f(x)的圖形')運行后屏幕顯示數(shù)據(jù)與擬合函數(shù)f的最大誤差Ew，平均誤差E1和均方根誤差E2及其數(shù)據(jù)點和擬合曲線y=f(x)的圖形.n=6f=4.71432.42861.14290.85721.57153.2858F=4.71432.42861.14290.85721.57153.2858Ew=0.4286E1=0.2857E2=0.30862.4實驗結(jié)論由圖可得，擬合出的曲線還是比擬合理的。同時可還得出在科學(xué)實驗和生產(chǎn)實踐中，往往會遇到這樣一種情況，即節(jié)點上的函數(shù)值并不是很精確的，這些函數(shù)值是由實驗或觀測得到的數(shù)據(jù)，不可防止地帶有測量誤差，如果要求所得的近似函數(shù)曲線精確無誤地通過所有的點(xi,yi),就會使曲線保存著一切測試誤差。當(dāng)個別數(shù)據(jù)的誤差較大時,插值效果顯然是不理想的。此外,由實驗或觀測提供的數(shù)據(jù)個數(shù)往往很多,如果用插值法,勢必得到次數(shù)較高的插值多項式，這樣計算起來很煩瑣。所以在一些預(yù)測的問題上，擬合還是比擬有優(yōu)勢的。試驗三復(fù)化Simpon積分3.1實驗?zāi)康膶W(xué)會并熟練掌握復(fù)化Simpson求積公式的的編程與應(yīng)用；進一步掌握Matlab數(shù)學(xué)軟件，尤其是提供計算積分的各種函數(shù)的使用方法；3.2根本原理將區(qū)間等分，子區(qū)間的長度為，在每個子區(qū)間上采用Simpson公式，在用Simpson公式時，還需要將子區(qū)間再二等分，因此有個分點，即。經(jīng)推導(dǎo)得到稱為復(fù)化Simpson值，稱如上式子為復(fù)化Simpson公式。類似于復(fù)化梯形公式，復(fù)化Simpson公式的余項為事實上，=。3.3實驗內(nèi)容functionI=S_quad(x,y)%復(fù)化Simpson求積公式，其中，X為向量，被積函數(shù)自變量的等距節(jié)點；%y為向量，被積函數(shù)在節(jié)點處的函數(shù)值；n=length(x);m=length(y);%積分自變量的節(jié)點數(shù)應(yīng)與它的函數(shù)值的個數(shù)相同ifn~=merror('ThelengthsofXandYmustbeequal');return;endifrem(n-1,2)~=0%如果n-1不能被2整除，那么調(diào)用復(fù)化梯形公式I=T_quad(x,y);return;endN=(n-1)/2;h=(x(n)-x(1))/N;a=zeros(1,n);fork=1:Na(2*k-1)=a(2*k-1)+1;a(2*k)=a(2*k)+4;a(2*k+1)=a(2*k+1)+1;endI=h/6*sum(a.*y);3.4數(shù)據(jù)來源根據(jù)以上用Matlab數(shù)學(xué)軟件所編代碼來上機求解實際數(shù)學(xué)問題，用復(fù)化Simpson公式求積分，在積分區(qū)間中點與點之間的間隔取為0.1。我們輸入如下數(shù)據(jù)：>>x=-1:0.1:1;>>y=exp(-x.^2);>>I=S_quad(x,y)運行得到結(jié)果為：I=1.49363.5實驗結(jié)論學(xué)過的常用的求積公式歸結(jié)為以下幾種：梯形公式及復(fù)化梯形公式；辛普森公式及其復(fù)化辛普森公式；科茨公式及其復(fù)化科茨公式；龍貝格求積公式等，在本實驗中重點研究復(fù)化辛普森公式，考察用數(shù)學(xué)軟件實現(xiàn)求積過程，經(jīng)過手算檢驗發(fā)現(xiàn)計算結(jié)果與計算機結(jié)果相一致。實驗四Runge-Kutta方法4.1實驗?zāi)康牧私釸unge-Kutta法的收斂性與穩(wěn)定性；學(xué)會用Matlab編程實現(xiàn)Runge-Kutta法解常微分方程；熟記Runge-Kutta法的各階表達形式；4.2根本原理考慮差分格式，視為在點處的斜率，視為在點處的預(yù)報斜率，假設(shè)參數(shù)及的取值使得上式具有2階精度，那么稱之為二階Runge-Kutta格式。另一種二階Runge-Kutta法的形式為它本質(zhì)上用和中間點的值進行計算，因此也稱為中點公式?？疾焖碾A經(jīng)典Runge-Kutta法是用和的加權(quán)平均值來近似。最常見的四階Runge-Kutta法有以上五個式子稱為標準四階Runge-Kutta公式，相應(yīng)的方法稱為標準四階Runge-Kutta法4.3實驗內(nèi)容function[x,y]=Runge_Kutta4(ydot_fun,x0,y0,h,N)%標準四階Runge_Kutta方法，其中，ydot_fun為一階微分方程的函數(shù)；%x0,y0為初始條件；%h為區(qū)間步長；%N為區(qū)間的個數(shù)；%x為Xn構(gòu)成的向量；%y為Yn構(gòu)成的向量；x=zeros(1,N+1);y=zeros(1,N+1);x(1)=x0;y(1)=y0;forn=1:Nx(n+1)=x(n)+h;k1=h*feval(ydot_fun,x(n),y(n));k2=h*feval(ydot_fun,x(n)+1/2*h,y(n)+1/2*k1);k3=h*feval(ydot_fun,x(n)+1/2*h,y(n)+1/2*k2);k4=h*feval(ydot_fun,x(n)+h,y(n)+k3);y(n+1)=y(n)+1/6*(k1+2*k2+2*k3+k4);end4.4數(shù)據(jù)來源>>ydot_fun=inline('x+y','x','y');>>[x,y]=Runge_Kutta4(ydot_fun,0,1,0.1,1)x=00.1000y=1.00001.1103顯然這里取到達改良方法中的計算結(jié)果〔實際上更精確，只是顯示結(jié)果相同〕，這說明，四階法有很好的計算效果。4.5實驗結(jié)論盡管我們熟知的四種方法，事實上，在求解微分方程初值問題，四階法是單步長中最優(yōu)秀的方法，通常都是用該方法求解的實際問題，計算效果還是很理想的。實驗五數(shù)值方法實際應(yīng)用5.1實驗?zāi)康膶W(xué)會函數(shù)的最正確平方逼近，理解最正確平方逼近的意義；學(xué)會數(shù)據(jù)的最小二乘擬合；學(xué)會Matlab提供的擬合方法，掌握非線性擬合方法，并學(xué)會用這些函數(shù)解決實際人口增長模型問題；5.2根本原理在生產(chǎn)和科學(xué)實驗中常常遇到這樣一類問題：要求依據(jù)一組實驗數(shù)據(jù)確定函數(shù)的近似表達式。一般給定數(shù)據(jù)的數(shù)量較大，且準確度不一定高，甚至個別點有很大的誤差，形象的稱為“噪聲”。假設(shè)用插值法求之，欲使?jié)M足插值條件，勢必將噪聲帶進近似函數(shù)，因而不能較好的描繪。數(shù)據(jù)擬合的最小二乘法是解決上述問題的一種有效、應(yīng)用廣泛的方法。設(shè)給定一組實驗數(shù)據(jù)及各點的權(quán)系數(shù)，要求在函數(shù)類中求函數(shù)滿足式中為中任意函數(shù)，稱這種函數(shù)近似表達式的方法為數(shù)據(jù)〔曲線〕的最小二乘擬合，稱為最小二乘解。在指定的函數(shù)類中求擬合數(shù)據(jù)的最小二乘解的關(guān)鍵在于求系數(shù)使得滿足以上式子，令那么是多元函數(shù)取得極小值的點。由極值點的必要充分條件，那么數(shù)據(jù)的最小二乘問題可轉(zhuǎn)化成求解方程組的問題，我們稱上式為法方程組，當(dāng)線性無關(guān)時，方程組有唯一解，于是有。5.3實驗內(nèi)容多項式擬合的通用程序：functionp=mafit(x,y,m)formatshort;A=zeros(m+1,m+1);fori=0:mforj=0:mA(i+1,j+1)=sum(x.^(i+j));endb(i+1)=sum(x.^i.*y);enda=A\b';p=fliplr(a');5.4數(shù)據(jù)來源例如：1949年到1994年我國人口數(shù)據(jù)資料如下：年份1949195419591964196919741979198419891994人口數(shù)5.46.06.77.08.19.19.810.311.311.8建模分析我國人口增長的規(guī)律,預(yù)報1999年我國人口數(shù)。用以下兩個簡單人口模型模型一：線性增長模型觀測值的模型：

線性模型：擬合的精度:

誤差平方和只需輸入：>>x=[1949195419591964196919741979198419891994];>>y=[5.466.778.19.19.810.311.311.8];>>p=mafit(x,y,1)p=0.1480-283.2320容易算出：=-283.232，=0.148模型：y=-283.2320+0.1480x預(yù)測1999年：>>x=1999;>>y=-283.2320+0.1480*xy=12.6200即1999年人口為12.6200億人。模型二：指數(shù)增長模型〔用簡單的線性最小二乘法〕用Matlab軟件計算得：>>x=[1949195419591964196919741979198419891994];>>y=[5.46.06.77.08.19.19.810.311.311.8];>>p=mafit(x,log(y),1)p=0.0179-33.1170

a=0.0179,b=-33.1170即:預(yù)測1999年人口：>>x=1999;>>y=exp(-33.1170+0.0179*x)y=14.3694即預(yù)測所得1999年人口數(shù)為14.3694億人。同時，執(zhí)行擬合圖形程序：>>x=[19491