數(shù)值分析報告實驗報告材料

1、實用標準文檔圖數(shù)值分析實驗報告姓名：張獻鵬學號：173511038專業(yè)：冶金工程班級：重冶二班1 拉格朗日插值 11.1 問題背景 11.2 數(shù)學模型 11.3 計算方法 11.4 數(shù)值分析 22復化辛普森求積公式 21.1 問題背景 21.2 數(shù)學模型 31.3 計算方法 31.4 數(shù)值分析 53 矩陣的LU分解 53.1 問題背景 53.2 數(shù)學模型 63.2.1 理論基礎 63.2.2 實例 63.3 計算方法 63.4 小組元的誤差 84 二分法求方程的根 94.1 問題背景 94.2 數(shù)學模型 94.3 計算方法 94.4 二分法的收斂性 105雅可比迭代求解方程組 101.1 問題

2、背景 101.2 數(shù)學模型 111.2.1 理論基礎 111.2.2 實例 111.3 計算方法 111.4 收斂性分析 126 Romber球積法 136.1 問題背景 136.2 數(shù)學模型： 136.2.1 理論基礎 136.2.2 實例 146.3 計算方法 146.4 誤差分析 157 幕法 157.1 問題背景 157.2 數(shù)學模型 157.2.1 理論基礎 157.2.2 實例 167.3 計算方法 167.4 誤差分析 178 改進歐拉法 178.1 問題背景 178.2 數(shù)學模型 178.2.1 理論基礎 178.2.2 實例 188.3 數(shù)學模型 188.4 誤差分析 19文

3、案大全1拉格朗日插值1.1 問題背景1f (x) 2對于函數(shù) 1 x , 5 X 5求拉格朗日插值。n 10,把f(x)和插值多項式的曲線畫在同一張圖上進行比較，觀察數(shù)值積分中的Lagrange插值。1.2 數(shù)學模型取等距差值節(jié)點？?=-5+10?n, ?=0, 1, . , n,構造n次lagrange插值多項式:?= Z2?=01?“(？)1 + ?2?(? ?街?+1(?當n=10時，十次插值多項式L0(x)以及函數(shù)f(x)的圖像可以由Matlab畫出1.3 計算方法f.m :function f= f( x ) f=1./(1+x.A2);endLagrange.mfunction y

4、=Lagrange(x0,y0,x);n=length(x0);m=length(x);for i=1:mz=x(i);s=0.0;for k=1:np=1.0;for j=1:n if j=k p=p*(z-x0(j)/(x0(k)-x0Q);endends=p*y0(k)+s;endy(i)=s；End拉格朗日插值的曲線:x=-5:1:5;y=1./(1+x.A2);xO=-5:0.001:5;y0=Lagrange(x,y,x0);y1=1./(1+x0A2);plot(x0,y0, hold on'b')plot(x0,y1, 'r')運行這個文件可以得

5、到如下拉格朗日圖形曲線:1.4 數(shù)值分析Lio(x)的斷誤差Rio(x)= f(x)-Lio(x)在區(qū)間卜5 , 5的兩端非常大。例如， Lio(4.8)=1.8O438，而f (4.8)=0.04160。這種現(xiàn)象稱之為龍格現(xiàn)象。不管n取多大，Runge 現(xiàn)象依然存在。因此，對函數(shù)作插值多項式時，必須小心處理，不能認為差值節(jié)點取得越多，差值 余項就越小。此外，當節(jié)點增多時，舍入誤差的影響不能低估。為了克服高次插值的不 足，應采用分段低次插值。2復化辛普森求積公式2.1 問題背景用復化Simpson公式計算定積分？ = /?2?我近似值，要求誤差限己=1/2X10- 7,利用其余項對算法做出步長

6、的事前估計；并將計算結果與精確值進行比較。2.2 數(shù)學模型將積分區(qū)間a,b分為n等分，h=(b-a)/ n, Xk=a+ k h , k=0, 1,門。在每個子區(qū) 問Xk, Xk+i上用Simpson公式可得：?-i?-i?(?+i)? ?x)dx = E ?(?)? E ?(?) + 4?7?j) +?=06 ?=02其中 Xk+1/2=Xk+1/2h。?-1?-1?-1Sn? = yE ? + 4?(?,?+1) + ?+1) =8? + 4 匯？7?g) + 2 匯？?鉤 + ?(?) ?=0?=0?=1此式即為復化Simpson公式。設f(x) C4a,b,由Simpson公式的誤差有

7、一 1? 5?= ? ?= U?=1o- 90(7) ?夕)(?為，?,?+1o則復化Simpson公式的余項為：一 一?-? 一 4 一?=- 288880-? ?夕)(?九? ? ?復化Simpson公式四階收斂。2.3 計算方法程序1(求f(X)的n階導數(shù))：syms xf=X*eXp(X)%定義函數(shù) f (x)n=input('輸入所求導數(shù)階數(shù):')f2=diff(f,x,n)%t f(x)的 n階導數(shù)程序2:clcclearsyms x%!義自變量xf=inline( 'x*exp(x)' , 'x' )%!義函數(shù)f(x)=x*exp(

8、x),換函數(shù)時只需換該函數(shù)表達式即可f2=inline( '(4*exp(x) + x*exp(x)' , 'x' )”定義f(x)的四階導數(shù)，輸入程序 1 里求出的f2即可f3= '-(4*exp(x) + x*exp(x)'%Ufminbnd ()函數(shù)求的是表達式的最小值，且要求表達式帶引號，故取負號，一邊求最大值e=5*10A(-8)嘛度要求值a=1瞅分下限b=2瞅分上限x1=fminbnd(f3,1,2)麻負的四階導數(shù)的最小值點，也就是求四階導數(shù)的最大值點對應的x值for n=2:1000000Rn=-(b-a)/180*(b-a)/(2

9、*n)A4*f2(x1)府算余項endif abs(Rn)<e%用余項進行判斷endh=(b-a)/nSn1=0Sn2=0break%符合要求時結束%R hfor k=0:n-1%求兩組連加和xk=a+k*hxk1=xk+h/2Sn1=Sn1+f(xk1)Sn2=Sn2+f(xk)endSn=h/6*(f(a)+4*Sn1+2*(Sn2-f(a)+f(b)%H Sn2力口了 k=0時的值，故減去f（a）z=exp(2)R=Sn-z%R已知值與計算值的差fprintf（'用Simpson公式計算的結果Sn='）disp(Sn)fprintf('等分數(shù) n='

10、)disp(n)fprintf（'已知值與計算值的誤差 R='）disp(R)運行結果為:E =2. 7284= e-0®用三inpwon公式l+算的結果£n=7. 3S91等分藪“24已知值與計算值的誤差R= 2.7281e-08A»2.4數(shù)值分析誤差分析：在上述計算中，若采用復化梯形公式，則可以知，其與精確值的誤差為 2.8300 Xe8,等分數(shù)為n=7019;而采用復化Simpson公式知與精確值的誤差為 2.7284 Xe-8,等分數(shù)為n=24o故與復化梯形公式相比，復化 Simpson公式誤差相對較小。收斂性分析：若lim Zj?=d?鈔

11、=/：?復化Simpson公式的余項是： ?° ,?-?4?=-痂？ ？夕)(?冽，？?？? ?可以看出誤差是 h4階，實際上若f(x) CC(a,b) , lim ?=4?(?)?？此復化 ?今Simpson公式是收斂的。穩(wěn)定性分析：由于求積公式中 A>0(i=0, 1,.，n)則求積公式是穩(wěn)定的。3矩陣的LU分解3.1 問題背景矩陣的LU分解主要用來求解線性方程組或者計算行列式。在使用初等行變換法求2-110經過E12（-3）、-1-2表示將i行元素與j行元1解線性方程組的過程中，系數(shù)矩陣的變化情況如下：A= 3-112-1日3(1)、E3(1/5)可得到0-53。00-1

12、2/5由上可知：E23(1/5) E13(1) E12(-3)A=U其中U就是上面矩陣A經過行變換后的上三角矩陣，Eij素互換的初等矩陣；Eij (k)表示將i行元素的k倍加到j行上 因此：A=E12(3) E13(-1) &(-1/5)12-110A = 310 = 31-1-1-2-1-1/50 12-10 0 -53 =LU1 00-12/5如果方陣A可以分解成單位下三角矩陣 L與上三角矩陣U的乘積，則式A=LU稱為A的LU分解或三角分解。3.2 數(shù)學模型3.2.1 理論基礎矩陣的LU分解在求解線性方程組時將十分簡便。 如對線性方程組Ax=b,設人=15 其LU分解。我們先求解方

13、程組Ly=b0由于L是下三角矩陣，則解向量y可以通過依次 求出其分量yb y2, , yn而求出，再求解方程組Ux=y。解向量x可以通過該方程組 依次求出分量Xn, Xn-1, , X2, X1而快速得出。于是由兩個方程組 Ux=y, Ly=b的求解 而給出LUx=Ly=b=A的解。?若矩陣A非奇異，則A能分解為LU的充分必要條件是A的順序主子行列式不為00則存在惟一的主對角線上元素全為1的下三角陣L與惟一的上三角陣U,使得A=LU10將矩陣20302045803.2.2 實例3080 進彳亍LU分解。1713.3 計算方法程序：clear allclcA=input('請輸入一個方陣

14、);輸入一個n階方陣n,n=size(A);L=zeros(n,n);U=zeros(n,n);for i=1:n % 將L的主對角線元素賦值1L(i,i)=1;endfor j=1:n 曬矩陣U勺第一行元素U戶A(1,j)；endfor k=2:n 減矩陣L的第一列元素L(k,1)=A(k,1)/U(1,1);endfor i=2:n 減L、U矩陣元素for j=i:ns=0;for t=1:i-1s=s+L(i,t)*U(t,j);endU(i,j)=A(i,j)-s；endfor k=i+1:nr=0;for t=1:i-1r=r+L(k,t)*U(t,i);endL(k,i)=(A(k

15、,i)-r)/U(i,i);endend喻出矢!陣L、ULU輸入一個方陣，輸出結果如下:畬令行葡匚請輸入一個方曬10 20 30.20 45 80.3(1 盯 171L 二00210341L020300S200D13.4 小組元的誤差例如線性方程組Ax = b ,其中：A= 0 1, b=1,可得理論解x=-1 。 1101如果系數(shù)矩陣被擾動成？= 10：01,可手算知：？N -20 0, ?=11101 10-20101 - 10-20 0若上述過程在計算機中進行，由浮點數(shù)運算規(guī)則可知，兩數(shù)相加時，大數(shù)吃掉小數(shù)， 則計算機中產生的矩陣為：?77 = ?7 ?"? = r10 201

16、1'一 .'0 0-10 -20 這時會發(fā)現(xiàn)？ = 101200,且據？ ？ x=b解出的理論解x =0,明顯不再等于前面的理論解。這說明LU分解是穩(wěn)定的，但是將LU分解用到解線性方程組上是不穩(wěn)定的。究其原 因，是因為？中的第一個主元10-20太小，導致第二個主元中的1與值10-20相差懸殊，出 現(xiàn)大數(shù)吃小數(shù)。為了避免上述危害，引入一種選主元手段，即在消去的過程中，通過適當?shù)倪x主元， 避免放大數(shù)據誤差。常用的選主元技術就是列選主元法 （除此之外還有全選主元法、對 角選主元法和隨機選主元法等）：對“< n階矩陣A,在確定第k個主元？?）?照>k）時，先從該列自主元位置

17、（k,jk）至列尾的所有元素中選擇絕對值最大的元素，與 蜜?袞換，然后將？岸？?"?？?雋為 零。繼續(xù)這個過程，直至將矩陣 A化為行階梯形。4二分法求方程的根4.1 問題背景在科學研究與工程計算中，常遇到方程(組)求根問題。若干個世紀以來，工程師 和數(shù)學家花了大量時用于探索求解方程 (組)，研究各種各樣的方程求解方法。對于方程 f(x)=0,當f(x)為線性函數(shù)時，稱f(x)=0為線性方程；當f(x)為非線性函數(shù)時，稱式 f(x)=0為非線性方程。對于線性方程(組)的求解，理論與數(shù)值求法的成果豐富；對于 非線性方程的求解，由于f(x)的多樣性，尚無一般的解析解法。當 f(x)為非線性

18、函數(shù) 時，若f(x)=0無解析解，但如果對任意的精度要求，設計迭代方程，數(shù)值計算出方程 的近似解，則可以認為求根的計算問題已經解決，至少能夠滿足實際要求。4.2 數(shù)學模型使用二分法求方程x”+x-1=0在0,1內的近似根(誤差10A-5)。二分法：二分法是最簡單的求根方法，它是利用連續(xù)函數(shù)的零點定理，將含根區(qū)間 逐次減半縮小，取區(qū)間的中點構造收斂點列 xk來逼近根x。4.3 計算方法二分法程序代碼：function y=erfen1(m,n,er) syms x xka=m;b=n;k=0;ff=xA3+x-1;while b-a>erxk=(a+b)/2;fx=subs(ff,x,xk

19、);fa=subs(ff,x,a);k=k+1;if fx=0y(k)=xk;break;elseif fa*fx<0b=xk;elsea=xk;endy(k)=xk;endplot(y, '.-');grid on在命令窗口下執(zhí)行:>> ab=ezfenl 1, Le-5): vpKabnS)實驗結果如下：可以得到迭代區(qū)間中點數(shù)列分布及圖像，數(shù)值如下：ans =0.5 , 0.75, 0.625, 0.6875, 0.65625, 0.671875, 0.6796875。0.68359375, 0.68164062, 0.68261719, 0.682128

20、91, 0.68237305, 0.68225098, 0.68231201, 0.68234253, 0.68232727, 0.6823349依據題目要求的精度，則需彳十七次，由實驗數(shù)據知x=0.6823349即為所求的根。4.4 二分法的收斂性二分法算法簡單，容易理解，且總是收斂的；但是二分法收斂速度太慢，浪費時間， 并且不能求復根跟偶數(shù)重根。5雅可比迭代求解方程組5.1 問題背景在自然科學和工程技術中有很多問題的解決常常歸結為解線性方程組，例如電學中的網絡問題，化學中的配平方程式模型問題，船體數(shù)學放樣中建立三次樣條函數(shù)問題，用最小二乘法求實驗數(shù)據的曲線擬合問題，非線性方程組求解問題，用

21、差分法或者有限 元法解常微分方程、偏微分方程邊值問題等都導致求解線性方程組。在實踐中，通常采 用數(shù)值方法來討論線性方程組 Ax=b的解，其中迭代法是一種重要方法。5.2 數(shù)學模型5.2.1 理論基礎雅可比迭代法：首先將方程組中的系數(shù)矩陣 A分解成三部分，即：A = L+D+U其中D為對角陣，L為下三角矩陣，U為上三角矩陣。之后確定迭代格式，XA(k+1) = B*XA(k)+f,(這里A表示的是上標，括號內數(shù)字即迭代次數(shù))，其中B稱為迭代矩陣，雅克比迭 代法中一般記為J (k= 0,1 ,)，再選取初始迭代向量 XA(0),開始逐次迭代。設Ax= b,其中A=D+L+UJ非奇異矩陣，且對角陣

22、D也非奇異，則當?shù)仃嘕的譜半徑p (J)<1時，雅克比迭代法收斂0?2?3?0?2？?+0？我?10 對于Ax=b,A=?01? ?+ ?1 ?2?1?2 ?3 0=L+D+U因為？?w 0(Jacob 假設)所以D可逆。故有：(L+D+U)x=bDx=-(L+U)x+bx=-D-1 (L+U)+D-1b5.2.2 實例用雅可比迭代法解方程組：430?2435-1 ?= 30 0 -14? -245.3 計算方法雅可比迭代法程序：function x,k,index=Jacobi(A,b,ep,it_max)if nargin<4it_max=100;end if nargin

23、<3 ep=1e-5;endn=length(A);k=0;x=zeros(n,1);y=zeros(n,1);index=1;while k<=it_max for i=1:n if abs(A(i,i)<1e-10index=0;return ; endy(i)=(b(i)-A(i,1:n)*x(1:n)+A(i,i)*x(i)/A(i,i); end if norm(y-x,inf)<epbreak; endk=k+1;x=y; end在命令窗口輸入的命令和結果如下圖：>> A=H 35 -I;0 -1 4 ;b=24 33 -24J ;ep=le-5

24、;it_wiaz=l00;xT kj indss x = J ac ob i (Aj b, 0p,'x 二4. 2000 2.4000 -b. 40M39index =15.4收斂性分析由上面運行的結果可知方程組的解為4.2000 , 2.4000 , -5.4000,迭代次數(shù)為39, 由index=1表示計算成功。雅克比迭代法的優(yōu)點明顯，計算公式簡單，每迭代一次只需計算一次矩陣和向量的乘法，且計算過程中原始矩陣 A始終不變，比較容易并行計算。然而這種迭代方式收斂 速度較慢，而且占據的存儲空間較大，所以工程中一般不直接用雅克比迭代法，而用其 改進方法。與逐次超松弛迭代法相比，雅可比迭代

25、法收斂速度相對較慢，逐次超松弛迭代法收 斂速度較快。由逐次超松弛迭代法求出的方程組的數(shù)值解與該方程組的精確解十分接近。 并且離散化后線性方程組的逐次超松弛迭代法的精確性較高，故相對于雅可比迭代法， 逐次超松弛迭代法更加廣泛地應用于實際， 可以用逐次超松弛迭代法求解高階稀疏線性 方程組。6 Romberg求積法6.1 問題背景復化求積方法對于提高精度是行之有效的方法，但復化公式的一個主要缺點在于要 事先估計出部長。若步長過大，則精度難于保證；若步長過小，則計算量又不會太大。 而用復化公式的截斷誤差來估計步長， 其結果是步長往往過小，而且f (刈和f(刈在 區(qū)間a，b上的上界M事估計是較為困難的。

26、在實際計算中通常采用變步長的方法，即 把步長逐次分半(也就是把步長二等分)，直到達到某種精度為止，這種方法就是Romberg 積分法的思想。6.2 數(shù)學模型：6.2.1 理論基礎記：TN0) TN ,將區(qū)間N等分的梯形值。T* SN ,將區(qū)間N等分的Simpson；TN2) CN ,將區(qū)間N等分的Coteso TN3)RN ,將區(qū)間N等分的Romberg(k)由其可構造一個序列'N ',次序列稱為Romberg字列，并滿足如下遞推關系:TiMm 3Nb af(ak 12N(2 k以上遞推公式就是4kT2(N1)TN4k 1k 1)一,k1,2,LRomberg積分遞推公式6.2

27、.2 實例利用Romberg積分法的思想，用龍貝格公式數(shù)值積分求函數(shù)L' ?密積分-1.26.3 計算方法龍貝格公式積分程序：function I,step=romberg(f,a,b,eps)%romberg.m為用龍貝格求積分%f為被積函數(shù)%a.b為積分區(qū)間的上下限%eps為積分結果精度%I為積分結果；step為積分的子區(qū)間數(shù)if (nargin=3)eps=1.0e-4;end;M=1;tol=10;k=0;T=zeros(1,1);h=b-a;T(1,1)=(h/2)*(subs(sym(f),findsym(sym(f),a)+subs(sym(f),findsym(sym(

28、f),b);while tol>epsk=k+1;h=h/2;Q=0;for i=1:Mx=a+h*(2*i-1);Q=Q+subs(sym,findsym(sym),x);endT(k+1,1)=T(k,1)/2+h*Q;M=2*M;for j=1:kT(k+1,j+1)=T(k+1,j)+(T(k+1,j)-T(k,j)/(4Aj-1);endtol=abs(T(k+1,j+1)-T(k,j);endI=T(k+1,k+1);step=k;輸入：» clea.r all;a=-L, 2：b=l. 2：fps=le-6:s±ep=rcinberg 1曰 jb,epm

29、)可以得到結果：I =0. 995326000000000st ep 二6.4 誤差分析Romberg求積公式是具有八階精度的算法，收斂且穩(wěn)定，比梯形公式，Simpson公式以及Cotes公式收斂的快。龍貝格公式的余項為:?,?(?= / ?(?)?/,?=- ?+22(?+1 )( ?+2?)2? !(? ?2?+3 ?夕?+2 )(?)Romberg積分公式高速有效，易于編程，適合于計算機計算。但他有一個主要的缺點是，每當把積分區(qū)間對分后，就要對被積函數(shù) ？?(?野算它在新分點處的值，而這些函 數(shù)值的個數(shù)是成倍增加的7曷法7.1 問題背景工程及物理中的許多實際問題需要計算矩陣的特征值及特征

30、向量，對于給定的n階實矩陣A,當n很小時用傳統(tǒng)的方法是可以的，但當 n稍大時，計算工作量將以驚人的 速度增大，并且由于計算存在誤差，用方程(入I-A)x=0求解十分困難。在實際應用中， 有的時候不一定需要求出矩陣的全部特征值和特征向量，例如只需要求出按模最大的或最小的特征值。求這類特征值的方法，通常采用迭代法，其中兩種是幕法和反幕法。7.2 數(shù)學模型7.2.1 理論基礎設An有n個線性相關的特征向量v1, v2,，vn,對應的特征值1, 2,n, 滿足| i| > | 2| n| ,因為v1 , v2,，vn為Cn的一組基，所以任給x(0) 0,x(0) aivii 1,所以有:n /

31、k (0)«k /、A x A (aivi)i 1nkai iVki 1nkai A vii 1n：覘1()kaivii 21若a10,則因1 知，當k充分大時Ak)x(0)1&V1 = CV1屬入1的特征向量另一方面,記max(x)=x,其中|Xi| 二 | x| ，則當k充分大時,max(Akx(0)max( 1ka1v1)1k max(a1v1)Ak 1.(0)kT"k 11max(A x ) max( 1 a1v1)1 max(a1v1)若a0,則因舍入誤差的影響，會有某次迭代向量在v1方向上的分量不為0,迭代 下去可求得入1及對應特征向量的近似值。7.2.

32、2 實例用幕法計算矩陣3-10?= -132023絕對值最大的特征值及相應的特征向量，取精度要求為e=10-4 07.3 計算方法幕法Matlab程序：function m,u,index=pow(A,ep,it_max) if nargin<3 it_max=100; end if nargin<2 ep=1e-5; end n=length(A);u=ones(n,1);index=0;k=0;m1=0;while k<=it_maxv=A*u;,i=max(abs(v);m=v(i);u=v/m;if abs(m-m1)<ep index=1; break ; e

33、ndm1=m;k=k+1;end在命令窗口輸入和輸出如下圖所示:» -l 0 -1 3 2:0 2 31:Mj uf ind«xj -po¥ (A, 11)U1L =5,-0,4136 LOCCO 0,9113index 7.4 誤差分析幕法程序可以用來計算矩陣絕對值最大的特征值及相應的特征向量。 幕法的缺點是 開始的時候并不知道矩陣是否有單一的主特征值，也不知道如何選擇Vo以保證它關于矩 陣特征向量的表達中包含一個與主特征值相關的非零特征向量。采用幕法和反幕法，求矩陣的最大和最小特征值，從原理上看，這兩種方法都是迭 代法，因此迭代初始向量的選擇對計算結果會產生一定影響，主要表現(xiàn)在收斂速度上。同時，原點位移m的選取也影響收斂的速度，但原點位移mo的適當選取依賴于對矩 陣A的大致了解。8改