《數(shù)值分析》實驗書2014

1、數(shù)值分析實驗指導(dǎo)書目 錄實驗?zāi)康?實驗基本要求2實驗一、誤差分析3一、實驗?zāi)康?二、算法實例3三、實驗任務(wù)10實驗二、插值法12一、實驗?zāi)康?2二、算法實例12三、實驗任務(wù)19四、思考題20實驗三、解線性方程組的直接法22一、實驗?zāi)康?2二、算法實例22三、實驗任務(wù)25四、思考題25實驗四、解線性方程組的迭代法26一、實驗?zāi)康?6二、算法實例26三、實驗任務(wù)30四、思考題30實驗五、常微分方程初值問題的數(shù)值解法32一、實驗?zāi)康?2二、算法實例32三、實驗任務(wù)42四、思考題4241實驗?zāi)康淖鳛閷嵺`性非常強的課程，安排上機實驗的目的，不僅是為了驗證教材和授課內(nèi)容，更重要的是，要通過實驗深入理解方法

2、的設(shè)計原理與處理問題的技巧，培養(yǎng)自行處理常規(guī)數(shù)值計算問題的能力和綜合運用知識分析、解決問題的能力。1、通過上機實驗加深課堂內(nèi)容的理解。數(shù)值分析的主要任務(wù)就是研究適合于在計算機上使用的數(shù)值計算方法及與此相關(guān)的理論。通過編程上機，就可以加深對方法運行過程的理解，同時在編程中領(lǐng)會和理解數(shù)值計算方法的計算要領(lǐng)和步驟，體會問題的條件和限制范圍，理解一般問題和特殊問題的區(qū)別。2、學(xué)會對數(shù)值計算結(jié)果的分析和處理。數(shù)值分析實驗不只是編寫程序得到一個數(shù)值結(jié)果，我們應(yīng)在掌握數(shù)值計算計算方法的基本原理和思想的同時，注意方法處理的技巧及其與計算機的密切結(jié)合，重視誤差分析、收斂性及穩(wěn)定性的討論。此外，還要注意算法能否在

3、計算機上實現(xiàn)，應(yīng)避免因數(shù)值方法選用不當(dāng)、程序設(shè)計不合理而導(dǎo)致超過計算機的存儲能力，或?qū)е掠嬎憬Y(jié)果精度不高等。3、要能靈活掌握各種數(shù)值計算方法。由于針對同一個問題可以選用不同的數(shù)值計算方法，我們要注意各種方法的使用條件。通過上機，比較各種方法間的異同及優(yōu)缺點，以便更好的使用不同的方法來解決實際問題，使計算機成為我們最好的工具。實驗基本要求一、上機前的準(zhǔn)備工作1、復(fù)習(xí)和掌握與本次實驗有關(guān)的教學(xué)內(nèi)容。2、根據(jù)本次實驗要求，在紙上編寫算法及上機的程序，并經(jīng)過人工模擬運行檢驗，減少不必要的錯誤，提高上機效率。切忌不編程序、不作人工檢查就進(jìn)行程序輸入，這只能使上機調(diào)試的難度增加，甚至可能帶來學(xué)習(xí)自信心的下

4、降，影響后續(xù)課程的學(xué)習(xí)。二、上機實驗步驟1、啟動開發(fā)環(huán)境；2、建立源程序文件，輸入源程序；3、編譯產(chǎn)生目標(biāo)程序，連接生成可執(zhí)行程序，運行程序，輸出結(jié)果；4、對數(shù)值計算結(jié)果進(jìn)行誤差分析，討論數(shù)值算法的收斂性與穩(wěn)定性；5、整理實驗報告。三、實驗報告實驗報告是記錄實驗工作全過程的技術(shù)文檔，實驗報告的撰寫是科學(xué)技術(shù)工作的一個組成部分。數(shù)值分析實驗報告包括下列要求：1、 實驗原理；2、 實驗內(nèi)容和要求；3、 數(shù)值算法描述，包括數(shù)據(jù)輸入、數(shù)據(jù)處理和數(shù)據(jù)輸出；4、算法的實現(xiàn)（1） 給出具體的計算實例，（2） 經(jīng)調(diào)試正確的源程序清單，（3） 對具體的數(shù)值例子給出數(shù)值結(jié)果；5、計算結(jié)果的誤差分析，算法的收斂性與

5、穩(wěn)定性的討論；6、實驗心得。實驗一、誤差分析誤差問題是數(shù)值分析的基礎(chǔ)，又是數(shù)值分析中一個困難的課題。在實際計算中，如果選用了不同的算法，由于舍入誤差的影響，將會得到截然不同的結(jié)果。因此，選取算法時注重分析舍入誤差的影響，在實際計算中是十分重要的。同時，由于在數(shù)值求解過程中用有限的過程代替無限的過程會產(chǎn)生截斷誤差，因此算法的好壞會影響到數(shù)值結(jié)果的精度。一、實驗?zāi)康?、 通過上機編程，復(fù)習(xí)鞏固以前所學(xué)程序設(shè)計語言及上機操作指令；2、 通過上機計算，了解誤差、絕對誤差、誤差界、相對誤差界的有關(guān)概念；3、 通過上機計算，了解舍入誤差所引起的數(shù)值不穩(wěn)定性。二、算法實例例1.1 用差商求在處導(dǎo)數(shù)的近似值。

6、取，=0.000 000 000 000 001和=0.000 000 000 000 000 1分別用MATLAB軟件計算，取十五位數(shù)字計算。解： 在MATLAB工作窗口輸入下面程序>>a=3;h=0.1;y=log(a+h)-log(a);yx=y/h運行后得yx = 0.32789822822991.將此程序中改為0.000 1，運行后得yx = 0.33332777790385.后者比前者好。再取h = 0.000 000 000 000 001，運行后得yx = 0.44408920985006，不如前者好。取h = 0.000 000 000 000 000 1，運行后

7、得yx = 0，算出的結(jié)果反而毫無價值。例1.2 分別求方程組在下列情況時的解，其中.（1）; （2）.解： （1） 首先將方程組化為同解方程，然后在MATLAB工作窗口輸入程序>> b=2,2'A=1,1;1,1.01; X=Ab運行后輸出當(dāng)時，的解為；（2）同理可得，當(dāng)時，的解為.例1.3 計算的近似值。解：泰勒級數(shù)e ,取，得. （1.1）這是一個無限過程，計算機無法求到精確值。只能在（1.1）取有限項時計算，再估計誤差。如果取有限項作為的值必然會有誤差，根據(jù)泰勒余項定理可知其截斷誤差為e.如果?。?.1）的前九項，輸入程序>> n=8;s=1;S =1;

8、fork=1:ns=s*k;S=S+1/s,ends, S,R=3/(s*(n+1)或>>S1=1+1+1/2+1/(1*2*3)+1/(1*2*3*4)+1/(1*2*3*4*5)+1/(1*2*3*4*5*6)+1/(1*2*3*4*5*6*7)+1/(1*2*3*4*5*6*7*8),R1=3/(1*2*3*4*5*6*7*8*9)運行后結(jié)果S =8.267195767195768e-006 R =2.71827876984127因為截斷誤差為所以e的近似值e2.718 28.例1.4 取作為的四舍五入近似值時，求其絕對誤差和相對誤差。解：在MATLAB工作窗口輸入程序>

9、;>juewu=exp(1)-2.71828運行后輸出結(jié)果為juewu = 1.828 459 045 505 326e-006例1.5 計算d 的近似值，并確定其絕對誤差和相對誤差。解 因為被積函數(shù)的原函數(shù)不是初等函數(shù)，故用泰勒級數(shù)求之。 , （1.5）這是一個無限過程，計算機無法求到精確值?？捎茫?.5）的前四項代替被積函數(shù)，得d)d=.根據(jù)泰勒余項定理和交錯級數(shù)收斂性的判別定理，得到絕對誤差= WU，在MATLAB命令窗口輸入計算程序如下：syms xf=1-x2/(1*2*3)+x4/(1*2*3*4*5)-x6/(1*2*3*4*5*6*7)y=int(f,x,0,pi/2),

10、y1=double(y)y11=pi/2-(pi/2)3/(3*3*2)+(pi/2)5/(5*5*4*3*2)-(pi/2)7/(7*7*6*5*4*3*2)inf=int(sin(x)/x,x,0,pi/2) ,infd=double(inf)WU =(pi/2)9/(9*9*8*7*6*5*4*3*2), R =infd-y11因為運行后輸出結(jié)果為： 1.370 762 168 154 49，=1.370 744 664 189 38，1.750 396 510 491 47e-005， WU= 1.782 679 830 970 664e-005. 所以，的絕對誤差為，故d。的相對誤差

11、為<0.007 3%.例1.6 設(shè)計三種算法求方程在的一個正根的近似值，并研究每種算法的誤差傳播情況.解：為解已知方程，我們可以設(shè)計如下三種算法，然后將計算結(jié)果與此方程的精確解比較，觀察誤差的傳播.算法1 將已知方程化為同解方程.取初值，按迭代公式 依次計算，結(jié)果列入表13中。算法2 將已知方程化為同解方程.取初值，按迭代公式 依次計算，結(jié)果列入表11中。算法3 將已知方程化為同解方程.取初值，按迭代公式為 依次計算，結(jié)果列入表11中。我們?yōu)檫@三種算法的計算編寫兩套MATLAB程序如下：（1）MATLAB主程序function k,juecha,xiangcha,xk= liti112(

12、x0,x1,limax)% 輸入的量-x0是初值, limax是迭代次數(shù)和精確值x;% 輸出的量-每次迭代次數(shù)k和迭代值xk,% -每次迭代的絕對誤差juecha和相對誤差xiangcha，x(1)=x0; for i=1:limax x(i+1)=fl(x(i);%程序中調(diào)用的fl.m juecha = abs(x(i)-x1);xiangcha = juecha /( abs(x(i)+eps);xk=x(i);k=i-1;(i-1),juecha,xiangcha,xkend（2）MATLAB調(diào)用函數(shù)程序及其計算結(jié)果 算法2的MATLAB調(diào)用函數(shù)程序function y1=fl(x)y1

13、=15/(2*x+1); 將MATLAB主程序和調(diào)用函數(shù)程序分別命名liti112.m和fl.m，分別保存為M文件，然后在MATLAB工作窗口輸入命令>> k,juecha,xiangcha,xk= liti112(2,2.5,100) 運行后輸出計算結(jié)果列入表11和表 1-2中。 將算法2的MATLAB調(diào)用函數(shù)程序的函數(shù)分別用y1=15-2*x2和y1=x-(2*x2+x-15)/(4*x+1)代替，得到算法1和算法3的調(diào)用函數(shù)程序，將其保存，運行后將三種算法的前8個迭代值列在一起（見表 1-1），進(jìn)行比較.將三種算法的對應(yīng)的絕對誤差和相對誤差的值列在一起（見表 1-2），進(jìn)行比

14、較。表 1-1 例1.6中三種算法的計算結(jié)果算 法迭代次數(shù)算法1的迭代結(jié)果算法2的迭代結(jié)果算法3的迭代結(jié)果022.000 000 002.000 000 00173.000 000 002.555 555 562-832.142 857 142.500 550 063-13 7632.837 837 842.500 000 064-378 840 3232.246 963 562.500 000 005-2.870 42.246 963 562.500 000 006-1.647 82.321 774 842.500 000 007-5.430 72.657 901 652.500 000 0

15、099-Inf2.500 000 012.500 000 00表 1-2 例1.6中三種算法計算結(jié)果的誤差算法迭代 次數(shù)算法1的誤差算法2的誤差算法3的誤差絕對誤差相對誤差絕對誤差相對誤差絕對誤差相對誤差00.500 000 000.250 000 000.500 000 000.250 000 000.500 000 000.250 000 0014.500 000 000.642 857 140.500 000 000.166 666 670.055 555 600.021 739 13285.500 000 001.030 120 480.357 142 860.1666 666 700

16、.000 550 100.000 219 97313 765.500 000.000 100 020.337 837 840.119 047 620.000 000 060.000 000 024378 840 3261.000 000 010.253 036 440.112 612 620.000 000 000.000 000 0052.870 399 8110.230 287 040.084 345 480061.647 839 0110.178 225 160.076 762 470075.430 746 8010.157 901 650.059 408 390099InfNaN0.0

17、00 000 010.000 000 0000 例1.7 求數(shù)的近似值。解 （1）直接用MATLAB命令 >> x=(715)*(sqrt(1+8(-19)-1)運行后輸出結(jié)果x = 0.問題出現(xiàn)在兩個相近的數(shù)與相減時，計算機運行程序>>sqrt(1+8(-19)-1運行后輸出結(jié)果 ans = 0.由于計算機硬件只支持有限位機器數(shù)的運算，因此在計算中可能引入和傳播舍入誤差.因為有效數(shù)字的嚴(yán)重?fù)p失，導(dǎo)致輸出的結(jié)果為0，計算機不能再與數(shù)繼續(xù)進(jìn)行真實的計算，所以，最后輸出的結(jié)果與的精確值不符。（2）如果化為，再用MATLAB命令 >> x=(715)*( (8(-

18、19)/(sqrt(1+8(-19)+1)運行后輸出結(jié)果 x = 1.6471e-005這是因為化為后，計算機運行程序>> x= (8(-19)/(sqrt(1+8(-19)+1)運行后的結(jié)果為x =3.4694e-018由于有效數(shù)字的損失甚少，所以運算的結(jié)果再與繼續(xù)計算，最后輸出的結(jié)果與的精確值相差無幾。例1.8 求數(shù)的近似值。解 （1）直接用MATLAB程序>> x=30;x1= sqrt(x2-1)運行后輸出結(jié)果x1 = 29.9833輸入MATLAB程序>> x=30; x1=29.9833;y=log(x-x1)運行后輸出結(jié)果y = -4.0923

19、（2）因為中的很大，如果采用倒數(shù)變換法，即 .輸入MATLAB程序>> x=30;y=-log(x+sqrt(x2-1)運行后輸出結(jié)果y = -4.0941（3）輸入MATLAB程序>> x=30; y=log(x-sqrt(x2-1)運行后輸出結(jié)果y = -4.0941可見，（2）計算的近似值比（1）的誤差小。參加計算的數(shù)，有時數(shù)量級相差很大.如果不注意采取相應(yīng)的措施，在它們的加減法運算中，絕對值很小的那個數(shù)經(jīng)常會被絕對值較大的那個數(shù)“吃掉”，不能發(fā)揮其作用，造成計算結(jié)果失真。例1.9 請在16位十進(jìn)制數(shù)值精度計算機上利用軟件MATLAB計算下面的兩個數(shù)和將計算結(jié)果與

20、準(zhǔn)確值比較，解釋計算結(jié)果。解 在MATLAB工作窗口輸入下面程序>> x=111111111111111+0.1+0.3, y=1111111111111111+0.1+0.3運行后輸出結(jié)果 x = 1.111111111111114e+014，y =1.111111111111111e+015從輸出的結(jié)果可以看出，x，而y.為什么僅僅比多一位1，而y呢？這是因為計算機進(jìn)行運算時，首先要把參加運算的數(shù)寫成絕對值小于1而“階碼”相同的數(shù)，這一過程稱為數(shù)的“對階”。在16位十進(jìn)制數(shù)值精度計算機上利用軟件MATLAB計算這兩個數(shù)，把運算的數(shù)寫成浮點規(guī)格化形式為在16位十進(jìn)制數(shù)值精度計算機

21、上，三項的數(shù)都表示為小數(shù)點后面16位數(shù)字的數(shù)與之積，所以，計算機沒有對數(shù)進(jìn)行截斷，而是按原來的三個數(shù)進(jìn)行計算。因此，計算的結(jié)果。而三項的數(shù)都表示寫成絕對值小于1而“階碼”都為的數(shù)以后，第一項的純小數(shù)的小數(shù)點后面有16位數(shù)字.但是，后兩項數(shù)的純小數(shù)的小數(shù)點后面有17位數(shù)字，超過了16位十進(jìn)制數(shù)值精度計算機的存儲量，計算機對后兩項的數(shù)都進(jìn)行截斷最后一位，即后兩項的數(shù)都是16位機上的零，再進(jìn)行計算，所以計算結(jié)果與實際不符。三、實驗任務(wù)對，計算定積分 .算法1：利用遞推公式 , ,取 .算法2：利用遞推公式 .注意到,取 .思考：從計算結(jié)果看，哪個算法是不穩(wěn)定的，哪個算法是穩(wěn)定的。實驗二、插值法插值法

22、是函數(shù)逼近的一種重要方法，它是數(shù)值積分、微分方程數(shù)值解等數(shù)值計算的基礎(chǔ)與工具，其中多項式插值是最常用和最基本的方法。拉格朗日插值多項式的優(yōu)點是表達(dá)式簡單明確，形式對稱，便于記憶，它的缺點是如果想要增加插值節(jié)點，公式必須整個改變，這就增加了計算工作量。而牛頓插值多項式對此做了改進(jìn)，當(dāng)增加一個節(jié)點時只需在原牛頓插值多項式基礎(chǔ)上增加一項，此時原有的項無需改變，從而達(dá)到節(jié)省計算次數(shù)、節(jié)約存儲單元、應(yīng)用較少節(jié)點達(dá)到應(yīng)有精度的目的。一、實驗?zāi)康?、理解插值的基本概念，掌握各種插值方法，包括拉格朗日插值和牛頓插值等，注意其不同特點；2、通過實驗進(jìn)一步理解并掌握各種插值的基本算法。二、算法實例1. 與插值有關(guān)

23、的MATLAB 函數(shù)(一) POLY2SYM 函數(shù)調(diào)用格式一：poly2sym (C)調(diào)用格式二：f1=poly2sym(C,'V') 或 f2=poly2sym(C, sym ('V') ),(二) POLYVAL 函數(shù)調(diào)用格式：Y = polyval(P,X)(三) POLY 函數(shù)調(diào)用格式：Y = poly (V)(四) CONV 函數(shù)調(diào)用格式：C =conv (A, B)(五) DECONV 函數(shù)調(diào)用格式：Q,R =deconv (B,A)(六) roots(poly(1:n)命令調(diào)用格式：roots(poly(1:n) (七) det(a*eye(siz

24、e (A) - A)命令調(diào)用格式：b=det(a*eye(size (A) - A)2.線性插值及其MATLAB程序例2.1 已知函數(shù)在上具有二階連續(xù)導(dǎo)數(shù)，且滿足條件。求線性插值多項式和函數(shù)值，并估計其誤差。解：輸入程序>> X=1,3;Y=1,2; l01= poly(X(2)/( X(1)- X(2), l11= poly(X(1)/( X(2)- X(1), l0=poly2sym (l01),l1=poly2sym (l11), P = l01* Y(1)+ l11* Y(2), L=poly2sym (P),x=1.5; Y = polyval(P,x)運行后輸出基函數(shù)l

25、0和l1及其插值多項式的系數(shù)向量P（略）、插值多項式L和插值Y為l0 =-1/2*x+3/2 l1 =1/2*x-1/2 L = 1/2*x+1/2 Y =1.2500輸入程序>> M=5;R1=M*abs(x-X(1)* (x-X(2)/2運行后輸出誤差限為 R1 = 1.8750例2.2 求函數(shù)e在上線性插值多項式，并估計其誤差。解： 輸入程序>> X=0,1; Y =exp(-X) , l01= poly(X(2)/( X(1)- X(2), l11= poly(X(1)/( X(2)- X(1), l0=poly2sym (l01),l1=poly2sym (l

26、11), P = l01* Y(1)+ l11* Y(2), L=poly2sym (P),運行后輸出基函數(shù)l0和l1及其插值多項式的系數(shù)向量P和插值多項式L為l0 = -x+1 l1 = x P = -0.6321 1.0000L =-1423408956596761/2251799813685248*x+1 輸入程序>> M=1;x=0:0.001:1; R1=M*max(abs(x-X(1).*(x-X(2)./2運行后輸出誤差限為 R1 = 0.1250.3.拋物線插值及其MATLAB程序例2.3 求將區(qū)間 0, /2 分成等份，用產(chǎn)生個節(jié)點，然后分別作線性插值函數(shù)和拋物線

27、插值函數(shù).用它們分別計算cos (/6) (取四位有效數(shù)字)，并估計其誤差。解： 輸入程序>> X=0,pi/2; Y =cos(X) ,l01= poly(X(2)/( X(1)- X(2), l11= poly(X(1)/( X(2)- X(1), l0=poly2sym (l01),l1=poly2sym (l11), P = l01* Y(1)+ l11* Y(2), L=poly2sym (P),x=pi/6; Y = polyval(P,x)運行后輸出基函數(shù)l0和l1及其插值多項式的系數(shù)向量P、插值多項式和插值為l0 =-5734161139222659/9007199

28、254740992*x+1l1 =5734161139222659/9007199254740992*xP = -0.6366 1.0000L =-5734161139222659/9007199254740992*x+1Y =0.6667輸入程序>> M=1;x=pi/6; R1=M*abs(x-X(1)*(x-X(2)/2運行后輸出誤差限為R1 = 0.2742.（2） 輸入程序>> X=0:pi/4:pi/2; Y =cos(X) ,l01= conv (poly(X(2),poly(X(3)/( X(1)- X(2)* ( X(1)- X(3), l11= co

29、nv (poly(X(1), poly(X(3)/( X(2)- X(1)* ( X(2)- X(3),l21= conv (poly(X(1), poly(X(2)/( X(3)- X(1)* ( X(3)- X(2),l0=poly2sym (l01),l1=poly2sym (l11),l2=poly2sym (l21),P = l01* Y(1)+ l11* Y(2) + l21* Y(3), L=poly2sym (P),x=pi/6; Y = polyval(P,x)運行后輸出基函數(shù)l01、l11和l21及其插值多項式的系數(shù)向量P、插值多項式L和插值Y為l0=22815502244

30、8185/281474976710656*x2-2150310427208497/1125899906842624*x+1l1=-228155022448185/140737488355328*x2+5734161139222659/2251799813685248*xl2=228155022448185/281474976710656*x2-5734161139222659/9007199254740992*xP = -0.3357 -0.1092 1.0000L=-6048313895780875/18014398509481984*x2-7870612110600739/72057594

31、037927936*x+1Y = 0.8508輸入程序>> M=1;x=pi/6; R2=M*abs(x-X(1)*(x-X(2) *(x-X(3)/6運行后輸出誤差限為R2 =0.0239.4.次拉格朗日（Lagrange）插值及其MATLAB程序例2.4 給出節(jié)點數(shù)據(jù)，作三次拉格朗日插值多項式計算，并估計其誤差。解：輸入程序>> X=-2,0,1,2; Y =17,1,2,17;p1=poly(X(1); p2=poly(X(2);p3=poly(X(3); p4=poly(X(4); l01= conv ( conv (p2, p3), p4)/( X(1)- X

32、(2)* ( X(1)- X(3) * ( X(1)- X(4), l11= conv ( conv (p1, p3), p4)/( X(2)- X(1)* ( X(2)- X(3) * ( X(2)- X(4),l21= conv ( conv (p1, p2), p4)/( X(3)- X(1)* ( X(3)- X(2) * ( X(3)- X(4),l31= conv ( conv (p1, p2), p3)/( X(4)- X(1)* ( X(4)- X(2) * ( X(4)- X(3),l0=poly2sym (l01),l1=poly2sym (l11),l2=poly2sym

33、 (l21), l3=poly2sym (l31),P = l01* Y(1)+ l11* Y(2) + l21* Y(3) + l31* Y(4),運行后輸出基函數(shù)l0，l1，l2和l3及其插值多項式的系數(shù)向量P（略）為l0 =-1/24*x3+1/8*x2-1/12*x，l1 =1/4*x3-1/4*x2-x+1l2 =-1/3*x3+4/3*x，l3 =1/8*x3+1/8*x2-1/4*x輸入程序>> L=poly2sym (P),x=0.6; Y = polyval(P,x)運行后輸出插值多項式和插值為L = x3+4*x2-4*x+1 Y = 0.2560.輸入程序&g

34、t;> syms M; x=0.6; R3=M*abs(x-X(1)*(x-X(2) *(x-X(3) *(x-X(4)/24運行后輸出誤差限為R3 =91/2500*M即 R3 ， .5. 拉格朗日插值及其誤差估計的MATLAB程序拉格朗日插值及其誤差估計的MATLAB主程序function y,R=lagranzi(X,Y,x,M)n=length(X); m=length(x);for i=1:m z=x(i);s=0.0; for k=1:n p=1.0; q1=1.0; c1=1.0;for j=1:n if j=kp=p*(z-X(j)/(X(k)-X(j); end q1=

35、abs(q1*(z-X(j);c1=c1*j; end s=p*Y(k)+s; end y(i)=s;endR=M*q1/c1;例 2.5 已知，用拉格朗日插值及其誤差估計的MATLAB主程序求的近似值，并估計其誤差。解：在MATLAB工作窗口輸入程序>> x=2*pi/9; M=1; X=pi/6 ,pi/4, pi/3;Y=0.5,0.7071,0.8660; y,R=lagranzi(X,Y,x,M)運行后輸出插值y及其誤差限R為y =0.6434 R =8.8610e-004.6.牛頓插值及其誤差估計的MATLAB主程序function y,R= newcz(X,Y,x,M

36、)n=length(X); m=length(x);for t=1:m z=x(t); A=zeros(n,n);A(:,1)=Y' s=0.0; p=1.0; q1=1.0; c1=1.0; for j=2:n for i=j:n A(i,j)=(A(i,j-1)-A(i-1,j-1)/(X(i)-X(i-j+1); end q1=abs(q1*(z-X(j-1);c1=c1*j; end C=A(n,n);q1=abs(q1*(z-X(n);for k=(n-1):-1:1C=conv(C,poly(X(k);d=length(C); C(d)=C(d)+A(k,k);end y(

37、k)= polyval(C, z);endR=M*q1/c1;例 2.6 已知，用牛頓插值法求的近似值，估計其誤差，并與例 6.2.6的計算結(jié)果比較.解 方法1（牛頓插值及其誤差估計的MATLAB主程序）輸入MATLAB程序>> x=2*pi/9;M=1; X=pi/6 ,pi/4, pi/3; Y=0.5,0.7071,0.8660; y,R= newcz(X,Y,x,M)運行后輸出y = R =0.6434 8.8610e-004方法2 （求牛頓插值多項式和差商的MATLAB主程序）輸入MATLAB程序>> x=2*pi/9; X=pi/6 ,pi/4, pi/3;

38、 Y=0.5,0.7071,0.8660; M=1; A,C,L,wcgs,Cw= newploy(X,Y), y=polyval(C,x)運行后輸出結(jié)果A = 0.5000 0 0 0.7071 0.7911 0 0.8660 0.6070 -0.3516C = -0.3516 1.2513 -0.0588L =-1583578379808357/4503599627370496*x2+1408883391907005/1125899906842624*x-132405829044691/2251799813685248wcgs =1/6*M*(x3-3/4*x2*pi+4012734077

39、357799/2251799813685248*x-7757769783530263/18014398509481984)Cw = 0.1667 -0.3927 0.2970 -0.0718y = 0.6434上述兩種方法計算y的結(jié)果相同.7. 牛頓插值法的MATLAB綜合程序求牛頓插值多項式、差商、插值及其誤差估計的MATLAB主程序function y,R,A,C,L=newdscg(X,Y,x,M)n=length(X); m=length(x);for t=1:m z=x(t); A=zeros(n,n);A(:,1)=Y's=0.0; p=1.0; q1=1.0; c1=1.

40、0; for j=2:n for i=j:n A(i,j)=(A(i,j-1)- A(i-1,j-1)/(X(i)-X(i-j+1); end q1=abs(q1*(z-X(j-1);c1=c1*j; end C=A(n,n);q1=abs(q1*(z-X(n);for k=(n-1):-1:1C=conv(C,poly(X(k);d=length(C);C(d)=C(d)+A(k,k);end y(k)= polyval(C, z);endR=M*q1/c1;L(k,:)=poly2sym(C);例 2.7 給出節(jié)點數(shù)據(jù)，作三階牛頓插值多項式，計算，并估計其誤差。解 首先將名為newdscg

41、.m的程序保存為M文件，然后在MATLAB工作窗口輸入程序>> syms M,X=-4,0,1,2; Y =27,1,2,17; x=-2.345; y,R,A,C,P=newdscg(X,Y,x,M)運行后輸出插值y及其誤差限公式R，三階牛頓插值多項式P及其系數(shù)向量C，差商的矩陣A如下y = 22.3211R =1323077530165133/562949953421312*M（即R =2.3503*M）A= 27.0000 0 0 0 1.0000 -6.5000 0 0 2.0000 1.0000 1.5000 0 17.0000 15.0000 7.0000 0.9167

42、C =0.9167 4.2500 -4.1667 1.0000P =11/12*x3+17/4*x2-25/6*x+1例 2.8 將區(qū)間 0,/2 分成等份，用產(chǎn)生個節(jié)點，求二階和三階牛頓插值多項式和.用它們分別計算sin (/7) (取四位有效數(shù)字)，并估計其誤差.解 首先將名為newdscg.m的程序保存為M文件，然后在MATLAB工作窗口輸入程序>> X1=0:pi/4:pi/2; Y1 =sin(X1); M=1; x=pi/7; X2=0:pi/6:pi/2; Y2 =sin(X2); y1,R1,A1,C1,P2=newdscg(X1,Y1,x,M), y2,R2,A2

43、,C2,P3=newdscg(X2,Y2,x,M)運行后輸出插值y1=和y2=及其誤差限R1和R2，二階和三階牛頓插值多項式P2和P3及其系數(shù)向量C1和C2，差商的矩陣A1和A2如下y1 =0.4548R1 =0.0282A1 = 0 0 0 0.7071 0.9003 0 1.0000 0.3729 -0.3357C1 = -0.3357 1.1640 0P2=-3024156947890437/9007199254740992*x2+163820246322191/140737488355328*xy2 =0.4345R2 =9.3913e-004A2 = 0 0 0 0 0.5000 0

44、.9549 0 0 0.8660 0.6991 -0.2443 0 1.0000 0.2559 -0.4232 -0.1139C2 =-0.1139 -0.0655 1.0204 0P3=-1025666884451963/9007199254740992*x3-4717668559161127/72057594037927936*x2+4595602396951547/4503599627370496*x三、實驗任務(wù)1、 已知函數(shù)表 0.56160 0.56280 0.56401 0.56521 0.82741 0.82659 0.82577 0.82495用二次拉格朗日插值多項式求時的函數(shù)

45、近似值。2、 已知函數(shù)表 0.4 0.55 0.65 0.8 0.9 0.41075 0.57815 0.69675 0.88811 1.02652用牛頓插值多項式求和。四、思考題1、 插值多項式的存在唯一性有何意義？2、 多項式插值是如何構(gòu)造的？截斷誤差如何表示？如何估計？3、 在插值區(qū)間上，隨著節(jié)點的增多，插值多項式是否越來越接近被插函數(shù)？實驗三、解線性方程組的直接法解線性方程組的直接法是指經(jīng)過有限步運算后能求得方程組精確解的方法。但由于實際計算中舍入誤差是客觀存在的，因而使用這類方法也只能得到近似解。目前較實用的直接法是古老的高斯消去法的變形，即主元素消去法及矩陣的三角分解法。引進(jìn)選主元

46、的技巧是為了控制計算過程中舍入誤差的增長，減少舍入誤差的影響。一般說來，列主元消去法及列主元三角分解法是數(shù)值穩(wěn)定的算法，它具有精確度較高、計算量不大和算法組織容易等優(yōu)點，是目前計算機上解中、小型稠密矩陣方程組可靠而有效的常用方法。一、實驗?zāi)康?、 了解求線性方程組的直接法的有關(guān)理論和方法；2、 會編制列主元消去法、LU分解法的程序；3、 通過實際計算，進(jìn)一步了解各種方法的優(yōu)缺點，選擇合適的數(shù)值方法。 二、算法實例1. 用高斯消元法解線性方程組的MATLAB程序function RA,RB,n,X=gaus(A,b)B=A b; n=length(b); RA=rank(A); RB=rank(

47、B);zhica=RB-RA;if zhica>0,disp('請注意：因為RA=RB，所以此方程組無解.')returnendif RA=RB if RA=ndisp('請注意：因為RA=RB=n，所以此方程組有唯一解.') X=zeros(n,1); C=zeros(1,n+1); for p= 1:n-1for k=p+1:n m= B(k,p)/ B(p,p); B(k,p:n+1)= B(k,p:n+1)-m* B(p,p:n+1);endend b=B(1:n,n+1);A=B(1:n,1:n); X(n)=b(n)/A(n,n); for q

48、=n-1:-1:1 X(q)=(b(q)-sum(A(q,q+1:n)*X(q+1:n)/A(q,q); endelse disp('請注意：因為RA=RB<n，所以此方程組有無窮多解.')endend例3.1 用高斯消元法和MATLAB程序求解下面的非齊次線性方程組，并且用逆矩陣解方程組的方法驗證。解 在MATLAB工作窗口輸入程序>> A=1 -1 1 -3; 0 -1 -1 1;2 -2 -4 6;1 -2 -4 1; b=1;0; -1;-1; RA,RB,n,X =gaus (A,b)運行后輸出結(jié)果請注意：因為RA=RB=n，所以此方程組有唯一解.X

49、= 0 -0.5000 0.5000 0RA =4RB =4n = 42用列主元消元法解線性方程組的MATLAB程序function RA,RB,n,X=liezhu(A,b)B=A b; n=length(b); RA=rank(A); RB=rank(B);zhica=RB-RA;if zhica>0,disp('請注意：因為RA=RB，所以此方程組無解.')returnendif RA=RB if RA=ndisp('請注意：因為RA=RB=n，所以此方程組有唯一解.') X=zeros(n,1); C=zeros(1,n+1); for p= 1:

50、n-1Y,j=max(abs(B(p:n,p); C=B(p,:);B(p,:)= B(j+p-1,:); B(j+p-1,:)=C;for k=p+1:n m= B(k,p)/ B(p,p); B(k,p:n+1)= B(k,p:n+1)-m* B(p,p:n+1);endend b=B(1:n,n+1);A=B(1:n,1:n); X(n)=b(n)/A(n,n); for q=n-1:-1:1 X(q)=(b(q)-sum(A(q,q+1:n)*X(q+1:n)/A(q,q); endelse disp('請注意：因為RA=RB<n，所以此方程組有無窮多解.')en

51、dend例3.2 用列主元消元法解線性方程組的MATLAB程序解方程組.解 在MATLAB工作窗口輸入程序>> A=0 -1 -1 1;1 -1 1 -3;2 -2 -4 6;1 -2 -4 1; b=0;1;-1;-1; RA,RB,n,X=liezhu(A,b)運行后輸出結(jié)果請注意：因為RA=RB=n，所以此方程組有唯一解.RA = 4，RB = 4，n = 4，X =0 -0.5 0.5 03用LU分解法解線性方程組的MATLAB程序function RA,RB,n,X,Y=LUjfcz(A,b)n n =size(A);B=A b; RA=rank(A); RB=rank(

52、B); for p=1:nh(p)=det(A(1:p, 1:p);endhl=h(1:n);for i=1:nif h(1,i)=0disp('請注意：因為A的r階主子式等于零，所以A不能進(jìn)行LU分解.A的秩RA和各階順序主子式值hl依次如下：') hl;RAreturnendendif h(1,i)=0 disp('請注意：因為A的各階主子式都不等于零，所以A能進(jìn)行LU分解.A的秩RA和各階順序主子式值hl依次如下：')X=zeros(n,1); Y=zeros(n,1); C=zeros(1,n);r=1:1;for p=1:n-1max1,j=max(abs(A(p:n,p); C=A(p,:); A(p,:)= A(j+P