數值分析實踐-3-前言當今時代，數學知識已經成為科學研究和工程技術領域基礎知識，計算問題是數學研究過程中一個重大難題，絕大數的問題都得不到準確的解或者計算量很大，需要借助計算機進行數值求解。數值分析就是主要介紹現代科學計算中常用的數值計算方法及其基本原理，研究并解決數值問題的近似解，是數學理論與計算機和實際問題的有機結合。數值分析課程需要和計算機緊密結合，各種數學運行軟件也是學習數值分析的必備課程。Matlab軟件就是其中一個重要的軟件，熟練掌握運用Matlab也成為學習數值分析的必備技能。

本文就是針對日常生活中的若干數值分析知識和數學模型問題以及這些相關問題在Matlab軟件上的運行進行討論。并且為體現出Matlab運行的優越性，與最后設計了C語言相比較同時，本文也在文章末尾總結了對于本次實習過程中所獲得的實驗心得。為以后的學習研究提供借鑒。摘要數值分析本文數值分析是研究分析用計算機求解數學計算問題的數值計算方法及其理論的學科，是數學的一個分支，它以數字計算機求解數學問題的理論和方法為研究對象。Matlab軟件是計算機數學經常使用的程序操作軟件。本文主要研究以下三個問題在Matlab軟件中實際操作問題：一是二分法和牛頓法的算法及實現；二是列主元Gauss消去法解方程組；三是用直接三角分解法解方程組。通過以上題目的練習掌握數值分析的相關知識以及Matlab軟件、C語言軟件的操作。

【關鍵詞】二分法；牛頓法；Gauss消去法；直接三角分解法

目錄前言………………….1摘要………………….2第一章、二分法和牛頓法的算法及現………………41.1二分法的算法及實現……………………….41.1.1問題的提出…………….41.1.2理論知識…………………41.1.3算法…………………………41.1.4程序…………………………41.1.5運行結果…………………51.2牛頓法的算法及實現…………51.2.1問題的提出………………51.2.2理論知識…………………61.2.3算法…………………………61.2.4程序…………………………71.2.5運行結果…………………81.3附注……………...9第二章、列主元高斯消去法的算法及實現………92.1問題的提出………………………92.2理論知識………………………….92.3算法…………………102.4程序…………………102.5結果…………………132.6附注…………………15第三章、矩陣直角三角法的算法及實現…………153.1問題的提出………………………153.2理論與算法………………………153.3程序…………………163.4結果…………………18總結…………………19參考文獻……………………………20

二分法和牛頓法的算法及實現1.1二分法的算法及實現1.1.1問題的提出用二分法計算方程在[1，2]內的根。()1.1.2理論知識對于二分法，其實質就是說對于給定的待求解的方程f(x)，其在[a,b]上連續，f(a)f(b)<0，且f(x)在[a,b]內僅有一個實根x*，取區間中點c，若，則c恰為其根，否則根據f(a)f(c)<0是否成立判斷根在區間[a,c]和[c,b]中的哪一個，從而得出新區間，仍稱為[a,b]。重復運行計算，直至滿足精度為止。這就是二分法的計算思想。1.1.3算法給定區間[a,b]，并設f(a)與f(b)符號相反，取ε為根的容許誤差，δ為|f(x)|的容許誤差。令c=(a+b)/2如果（c-a）<ε或|f(c)|<δ，則輸出c,結束；否則執行③；如果f(a)f(c)>0，則令b=c，重復①②③。1.1.4程序Matlab的M文件編寫：clear%%%給定求解區間b=1.5;a=0;%%%誤差R=1;k=0;%迭代次數初值while(R>5e-6);c=(a+b)/2;iff12(a)*f12(c)>0;a=c;elseb=c;endR=b-a;%求出誤差k=k+1;endx=c%給出解1.1.5運行結果x=1.40441513061523；f(x)=-3.797205105904311e-007；k=18；由f(x)知結果滿足要求，但迭代次數比較多，方法收斂速度比較慢。1.2牛頓法的算法及實現1.2.1問題的提出用Newton法求解下列方程：x0=0.51.2.2理論知識Newton法通常預先要給出一個猜測初值x0，然后根據其迭代公式產生逼近解x*的迭代數列{xk}，這就是Newton法的思想。當x0接近x*時收斂很快，但是當x0選擇不好時，可能會發散，因此初值的選取很重要。另外，若將該迭代公式改進為其中r為要求的方程的根的重數，這就是改進的Newton法，當求解已知重數的方程的根時，在同種條件下其收斂速度要比Newton法快的多。1.2.3算法給定初始值x0如果f'(x)計算。若|x1-x0|<ε或|f(令x01.2.4程序Matlab的M文件編寫：clear%%%%輸入函數f=input('請輸入需要求解函數>>','s')%%%求解f(x)的導數df=diff(f);%%%改進常數或重根數miu=2;%%%初始值x0x0=input('inputinitialvaluex0>>');k=0;%迭代次數max=100;%最大迭代次數R=eval(subs(f,'x0','x'));%求解f(x0)，以確定初值x0時否就是解while(abs(R)>1e-8)x1=x0-miu*eval(subs(f,'x0','x'))/eval(subs(df,'x0','x'));R=x1-x0;x0=x1;k=k+1;if(eval(subs(f,'x0','x'))<1e-10);breakendifk>max;%如果迭代次數大于給定值，認為迭代不收斂，重新輸入初值ss=input('mayberesultiserror,chooseanewx0,y/n?>>','s');ifstrcmp(ss,'y')x0=input('inputinitialvaluex0>>');k=0;elsebreakendendendk;%給出迭代次數x=x0；%給出解1.2.5結果x=0.56714329040978；f(x)=2.220446049250313e-016；k=4；由f(x)知結果滿足要求，而且又迭代次數只有4次看出收斂速度很快。1.3附注本實驗采用Matlab的M文件編寫。其中待求解的方程寫成function的方式，如下functiony=f(x);y=-x*x-sin(x);對于二分法，只要能夠保證在給定的區間內有根，使能夠收斂的，當時收斂的速度和給定的區間有關，二且總體上來說速度比較慢。Newton法，收斂速度要比二分法快，但是最終其收斂的結果與初值的選取有關，初值不同，收斂的結果也可能不一樣，也就是結果可能不時預期需要得結果。第二章、列主元高斯消去法的算法及實現2.1問題的提出求解方程。其中為一小數，當時，分別采用列主元和不列主元的Gauss消去法求解，并比較結果。2.2理論知識由于一般線性方程在使用Gauss消去法求解時，從求解的過程中可以看到，若=0，則必須進行行交換，才能使消去過程進行下去。有的時候即使0，但是其絕對值非常小，由于機器舍入誤差的影響，消去過程也會出現不穩定得現象，導致結果不正確。因此有必要進行列主元技術，以最大可能的消除這種現象。這一技術要尋找行r，使得并將第r行和第k行的元素進行交換，以使得當前的的數值比0要大的多。2.3算法消元過程對k=1,2,…,n-1,進行如下步驟。選主元，記若很小，這說明方程的系數矩陣嚴重病態，給出警告，提示結果可能不對。交換增廣陣A的r，k兩行的元素。(j=k,…,n+1)計算消元(i=k+1,…,n;j=k+1,……,n+1)回代過程對k=n,n-1,…,1,進行如下計算至此，完成了整個方程組的求解。2.4程序Matlab的M文件編寫：cleara=input('輸入系數陣：>>\n')b=input('輸入列陣b：>>\n')n=length(b);A=[ab]x=zeros(n,1);%%%函數主體fork=1:n-1;%%%是否進行主元選取ifabs(A(k,k))<yipusilong;%事先給定的認為有必要選主元的小數yzhuyuan=1;elseyzhuyuan=0;endifyzhuyuan;%%%%選主元t=A(k,k);forr=k+1:n;ifabs(A(r,k))>abs(t)p=r;elsep=k;endend%%%交換元素ifp~=k;forq=k:n+1;s=A(k,q);A(k,q)=A(p,q);A(p,q)=s;endendend%%%判斷系數矩陣是否奇異或病態非常嚴重ifabs(A(k,k))<yipusilongdisp(‘矩陣奇異，解可能不正確’)end%%%%計算消元,得三角陣forr=k+1:n;m=A(r,k)/A(k,k);forq=k:n+1;A(r,q)=A(r,q)-A(k,q)*m;endendend%%%%求解xx(n)=A(n,n+1)/A(n,n);fork=n-1:-1:1;s=0;forr=k+1:n;s=s+A(k,r)*x(r);endt=(A(k,n+1)-s)x(k)=(A(k,n+1)-s)/A(k,k)end2.5結果記Emax為求出的解代入方程后的最大誤差，按要求，計算結果如下：當時，不選主元和選主元的計算結果如下，其中前一列為不選主元結果，后一列為選主元結果，下同。0.999999347683910.999999347826512.000002174219722.000002173911632.999997608594512.99999760869721Emax=9.301857062382624e-010，0此時，由于不是很小，機器誤差就不是很大，由Emax可以看出不選主元的計算結果精度還可以，因此此時可以考慮不選主元以減少計算量。當時，不選主元和選主元的計算結果如下1.000017846308770.999999999993481.999980097208072.000000000021743.000006634247312.99999999997609Emax=2.036758973744668e-005，0此時由Emax可以看出不選主元的計算精度就不好了，誤差開始增大。當時，不選主元和選主元的計算結果如下1.421085471520201.000000000000001.666666666666662.000000000000003.11111111111111300000000000000Emax=0.70770085900503，0此時由Emax可以看出，不選主元的結果應該可以說是不正確了，這是由機器誤差引起的。當時，不選主元和選主元的計算結果如下NaN1NaN2NaN3Emax=NaN,0不選主元時，程序報錯：Warning:Dividebyzero.。這是因為機器計算的最小精度為10-15，所以此時的就認為是0，故出現了錯誤現象。而選主元時則沒有這種現象，而且由Emax可以看出選主元時的結果應該是精確解。2.6附注采用Gauss消去法時，如果在消元時對角線上的元素始終較大（假如大于10-5），那么本方法不需要進行列主元計算，計算結果一般就可以達到要求，否則必須進行列主元這一步，以減少機器誤差帶來的影響，使方法得出的結果正確。第三章、矩陣直角三角法的算法及實現3.1問題的提出用直角三角分解求下列線性方程組的解：13.2理論與算法將方程組Ax=b中的A分解為A=LU，其中L為單位下三角矩陣，U為上三角矩陣，則方程Ax=b化為解兩個方程組Ly=b，Ux=y，具體算法如下：對j=1,2,3,…,n計算u對i=2,3,…,n計算l②對k=2,3,…,n：對j=k,k+1,…,n計算u對i=k+1,k+2,…,n計算l③y1=byxn=yx注：由于計算u的公式與計算y的公式形式上一樣，故可直接對增廣矩陣A|b施行算法②和③，此時U的第(n+1)列元素即為y。3.3程序C語言程序如下：#include(stdio.h)Voidmain(){Floatx[3];IntI;Floata[3][4]={1/4,1/5,1/6,9,1/3,1/4,1/5,8,1/2,1,2};VoidDirectLU(float*,int,float[]);DirectLU(a[0],3,x);For(i=0;i<=3;i++)printf(“x[%d]=%f\n”,I,x[i]);}VoidDirectLU(float*u,intn,floatx[]){IntI,r,k;for(r=0;r<=n-1;r++){for(i=r;i<=n;i++)for(k=0;k<=r-1;k++)*(u+r*(n+1)+i)-=*(u+r*(n+1)+k)*(*(u+k*(n+1)+i));for(i=r+1;i<=n-1;i++){for(k=0;k<=r-1;k++)*(u+i*(n+1)+r)-=*(u+i*(n+1)+k)*(*(u+k*(n+1)+r));*(u+i*(n+1)+r)/=*(u+r*(n+1)+r);}}for(i=n-1