《數(shù)值分析》06實驗指導(dǎo)書正文

1、精選優(yōu)質(zhì)文檔-傾情為你奉上實驗?zāi)康淖鳛閷嵺`性非常強的課程，安排上機實驗的目的，不僅是為了驗證教材和授課內(nèi)容，更重要的是，要通過實驗深入理解方法的設(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)在掌握

2、數(shù)值計算計算方法的基本原理和思想的同時，注意方法處理的技巧及其與計算機的密切結(jié)合，重視誤差分析、收斂性及穩(wěn)定性的討論。此外，還要注意算法能否在計算機上實現(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)過人工

3、模擬運行檢驗，減少不必要的錯誤，提高上機效率。切忌不編程序、不作人工檢查就進行程序輸入，這只能使上機調(diào)試的難度增加，甚至可能帶來學(xué)習(xí)自信心的下降，影響后續(xù)課程的學(xué)習(xí)。二、上機實驗步驟1、啟動開發(fā)環(huán)境；2、建立源程序文件，輸入源程序；3、編譯產(chǎn)生目標(biāo)程序，連接生成可執(zhí)行程序，運行程序，輸出結(jié)果；4、對數(shù)值計算結(jié)果進行誤差分析，討論數(shù)值算法的收斂性與穩(wěn)定性；5、整理實驗報告。三、實驗報告實驗報告是記錄實驗工作全過程的技術(shù)文檔，實驗報告的撰寫是科學(xué)技術(shù)工作的一個組成部分。數(shù)值分析實驗報告包括下列要求：1、 實驗原理；2、 實驗內(nèi)容和要求；3、 數(shù)值算法描述，包括數(shù)據(jù)輸入、數(shù)據(jù)處理和數(shù)據(jù)輸出；4、算法

4、的實現(xiàn)（1） 給出具體的計算實例，（2） 經(jīng)調(diào)試正確的源程序清單，（3） 對具體的數(shù)值例子給出數(shù)值結(jié)果；5、計算結(jié)果的誤差分析，算法的收斂性與穩(wěn)定性的討論；6、實驗心得。實驗項目實驗一、誤差分析誤差問題是數(shù)值分析的基礎(chǔ)，又是數(shù)值分析中一個困難的課題。在實際計算中，如果選用了不同的算法，由于舍入誤差的影響，將會得到截然不同的結(jié)果。因此，選取算法時注重分析舍入誤差的影響，在實際計算中是十分重要的。同時，由于在數(shù)值求解過程中用有限的過程代替無限的過程會產(chǎn)生截斷誤差，因此算法的好壞會影響到數(shù)值結(jié)果的精度。一、實驗?zāi)康?、 通過上機編程，復(fù)習(xí)鞏固以前所學(xué)程序設(shè)計語言及上機操作指令；2、 通過上機計算，了

5、解誤差、絕對誤差、誤差界、相對誤差界的有關(guān)概念；3、 通過上機計算，了解舍入誤差所引起的數(shù)值不穩(wěn)定性。二、算法實例例1.1 用差商求在處導(dǎo)數(shù)的近似值。取，=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.991.將此程序中改為0.000 1，運行后得yx = 0.385.后者比前者好。再取h = 0.000 000 000 000 001，運行后得yx = 0.

6、006，不如前者好。取h = 0.000 000 000 000 000 1，運行后得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ù)泰勒余項定理可知其截斷誤差

7、為e.如果?。?.1）的前九項，輸入程序>> n=8;s=1;S =1;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.5768e-006 R =2.127因為截斷誤差為所以e的近似值e2.718 28.例1.4 取作為的四舍五入近似值時，求其絕對誤差和相對誤差。

8、解：在MATLAB工作窗口輸入程序>>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*

9、7)y=int(f,x,0,pi/2),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

10、. 所以，的絕對誤差為，故d。的相對誤差為<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,x

11、iangcha,xk= liti112(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ù)程

12、序function y1=fl(x)y1=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），進行比較.將三種算法的對應(yīng)的絕對誤差和相對

13、誤差的值列在一起（見表 1-2），進行比較。表 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.6

14、57 901 652.500 000 0099-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

15、142 860.1666 666 700.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

16、 408 390099InfNaN0.000 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ù)進行真實的計算，所以，最后輸出的結(jié)果與的精確值不符。（2）如果化為，再用MATLAB命令 >

17、;> x=(715)*( (8(-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-x

18、1)運行后輸出結(jié)果y = -4.0923（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位十進制數(shù)值精度計算機上利用軟件M

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

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

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

22、數(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(size (A) - A)命令調(diào)用格式：b=det(a

23、*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ù)l0和l1及其插值多項式的系數(shù)向量P（略）、插值多

24、項式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 (l11), P = l01* Y(1)+ l11*

25、 Y(2), L=poly2sym (P),運行后輸出基函數(shù)l0和l1及其插值多項式的系數(shù)向量P和插值多項式L為l0 = -x+1 l1 = x P = -0.6321 1.0000L =-96761/85248*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ù)和拋物線插值函數(shù).用它們分別計算cos (/6) (取四位有效數(shù)字)，并估計其誤差。解： 輸入程序&

26、gt;> 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 =-22659/40992*x+1l1 =22659/40992*xP = -0.6366 1.0000L =-22659/40992*x+1Y

27、 =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= conv (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),l

28、0=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=8185/0656*x2-08497/42624*x+1l1=-8185/5328*x2+22659/85248*xl2=8185/0656*x2-22659/40992*xP = -0.3357 -0.1092 1.0000L=-80875/*

29、x2-00739/*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

30、(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

31、 (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

32、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=

33、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

34、)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(

35、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;

36、 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 =-08357/70496*x2+07005/42624*x-4691/85248wcgs =1/6*M*(x3-3/4*x2*pi+57799/85248*x-30263/)Cw = 0.1667 -0.3927 0.2970 -0.0718y = 0.6434上述兩種方法計算y的結(jié)

37、果相同.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.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

38、=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.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及

39、其誤差限公式R，三階牛頓插值多項式P及其系數(shù)向量C，差商的矩陣A如下y = 22.3211R =65133/1312*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.9167C =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ù)字)，并估計其誤差.解 首先

40、將名為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,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

41、1.0000 0.3729 -0.3357C1 = -0.3357 1.1640 0P2=-90437/40992*x2+2191/5328*xy2 =0.4345R2 =9.3913e-004A2 = 0 0 0 0 0.5000 0.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=-51963/40992*x3-61127/*x2+51547/70496*x三、實驗任務(wù)1、 已知函數(shù)表 0.56160 0.56280 0.56401 0.56521 0.

42、82741 0.82659 0.82577 0.82495用二次拉格朗日插值多項式求時的函數(shù)近似值。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ù)？實驗三、數(shù)據(jù)擬合法曲線擬合的最小二乘法是計算機數(shù)據(jù)處理的重要內(nèi)容，也是函數(shù)逼近的另一種重要方法，它在工程技術(shù)中有著廣泛的應(yīng)用。對實際問題而言，擬合曲線的選擇是一個

43、極其重要而又比較困難的問題，必要時可由草圖觀察選取幾種不同類型的擬合曲線，再以其偏差小者為優(yōu)，經(jīng)檢驗后再決定最后的取舍。一、實驗?zāi)康?、 理解數(shù)據(jù)擬合的基本概念、基本方法；2、 掌握最小二乘法的基本原理，并會通過計算機解決實際問題；3、了解超定方程組的最小二乘解法。二、算法實例例3.1 給出一組數(shù)據(jù)點列入表31中，試用線性最小二乘法求擬合曲線，并估計其誤差，作出擬合曲線。表31 例3.1的一組數(shù)據(jù)xi-2.5 -1.7 -1.1 -0.8 0 0.1 1.5 2.7 3.6yi-192.9 -85.50 -36.15 -26.52 -9.10 -8.43 -13.12 6.50 68.04解

44、（1）在MATLAB工作窗口輸入程序>> x=-2.5 -1.7 -1.1 -0.8 0 0.1 1.5 2.7 3.6; y=-192.9 -85.50 -36.15 -26.52 -9.10 -8.43 -13.12 6.50 68.04;plot(x,y,'r*'),legend('實驗數(shù)據(jù)(xi,yi)')xlabel('x'), ylabel('y'),title('例3.1的數(shù)據(jù)點(xi,yi)的散點圖')運行后屏幕顯示數(shù)據(jù)的散點圖（略）。 （3）編寫下列MATLAB程序計算在處的函數(shù)值，即

45、輸入程序>> syms a1 a2 a3 a4x=-2.5 -1.7 -1.1 -0.8 0 0.1 1.5 2.7 3.6; fi=a1.*x.3+ a2.*x.2+ a3.*x+ a4運行后屏幕顯示關(guān)于a1,a2, a3和a4的線性方程組fi = -125/8*a1+25/4*a2-5/2*a3+a4, -4913/1000*a1+289/100*a2-17/10*a3+a4, -1331/1000*a1+121/100*a2-11/10*a3+a4, -64/125*a1+16/25*a2-4/5*a3+a4, a4,1/1000*a1+1/100*a2+1/10*a3+a4

46、, 27/8*a1+9/4*a2+3/2*a3+a4, 19683/1000*a1+729/100*a2+27/10*a3+a4, 5832/125*a1+324/25*a2+18/5*a3+a4編寫構(gòu)造誤差平方和的MATLAB程序>> y=-192.9 -85.50 -36.15 -26.52 -9.10 -8.43 -13.12 6.50 68.04;fi=-125/8*a1+25/4*a2-5/2*a3+a4, -4913/1000*a1+289/100*a2-17/10*a3+a4, -1331/1000*a1+121/100*a2-11/10*a3+a4, -64/125

47、*a1+16/25*a2-4/5*a3+a4,a4, 1/1000*a1+1/100*a2+1/10*a3+a4, 27/8*a1+9/4*a2+3/2*a3+a4, 19683/1000*a1+729/100*a2+27/10*a3+a4, 5832/125*a1+324/25*a2+18/5*a3+a4;fy=fi-y; fy2=fy.2; J=sum(fy.2)運行后屏幕顯示誤差平方和如下J=(-125/8*a1+25/4*a2-5/2*a3+a4+1929/10)2+(-4913/1000*a1+289/100*a2-17/10*a3+a4+171/2)2+(-1331/1000*a1

48、+121/100*a2-11/10*a3+a4+723/20)2+(-64/125*a1+16/25*a2-4/5*a3+a4+663/25)2+(a4+91/10)2+(1/1000*a1+1/100*a2+1/10*a3+a4+843/100)2+(27/8*a1+9/4*a2+3/2*a3+a4+328/25)2+(19683/1000*a1+729/100*a2+27/10*a3+a4-13/2)2+(5832/125*a1+324/25*a2+18/5*a3+a4-1701/25)2為求使達(dá)到最小，只需利用極值的必要條件 ，得到關(guān)于的線性方程組，這可以由下面的MATLAB程序完成，即

49、輸入程序 >> syms a1 a2 a3 a4J=(-125/8*a1+25/4*a2-5/2*a3+a4+1929/10)2+(-4913/1000*a1+289/100*a2-17/10*a3+a4.+171/2)2+(-1331/1000*a1+121/100*a2-11/10*a3+a4+723/20)2+(-64/125*a1+16/25*a2-4/5*a3+a4+663/25)2+(a4+91/10)2+(1/1000*a1+1/100*a2+1/10*a3+a4+843/100)2+(27/8*a1+9/4*a2+3/2*a3+a4+328/25)2+(19683/

50、1000*a1+729/100*a2+27/10*a3+a4-13/2)2+(5832/125*a1+324/25*a2+18/5*a3+a4-1701/25)2; Ja1=diff(J,a1); Ja2=diff(J,a2); Ja3=diff(J,a3); Ja4=diff(J,a4);Ja11=simple(Ja1), Ja21=simple(Ja2), Ja31=simple(Ja3), Ja41=simple(Ja4),運行后屏幕顯示J分別對a1, a2 ,a3 ,a4的偏導(dǎo)數(shù)如下Ja11=/10000*a1+/25000*a2+/2500*a3+23667/250*a4-/625J

51、a21=/25000*a1+/2500*a2+23667/250*a3+67*a4+/625Ja31=/2500*a1+23667/250*a2+67*a3+18/5*a4-/125Ja41=23667/250*a1+67*a2+18/5*a3+18*a4+14859/25解線性方程組Ja11 =0，Ja21 =0，Ja31 =0，Ja41 =0，輸入下列程序>>A=/10000, /25000, /2500, 23667/250; /25000, /2500, 23667/250, 67; /2500, 23667/250, 67, 18/5; 23667/250, 67, 18

52、/5, 18;B=/625, -/625, /125, -14859/25; C=B/A, f=poly2sym(C)運行后屏幕顯示擬合函數(shù)f及其系數(shù)C如下C = 5.0911 -14.1905 6.4102 -8.2574f=5759/5328*x3-57579/1312*x2+77693/0656*x-13215/1312 故所求的擬合曲線為.（4）編寫下面的MATLAB程序估計其誤差，并作出擬合曲線和數(shù)據(jù)的圖形.輸入程序>> xi=-2.5 -1.7 -1.1 -0.8 0 0.1 1.5 2.7 3.6; y=-192.9 -85.50 -36.15 -26.52 -9.10 -8.43 -13.12 6.50 68.04;n=length(xi); f=5.0911.*xi.3-14.1905.*xi.2+6.4102.*xi -8.2574;x=-2.5:0.01: 3.6; F=5.0911.*x.3-14.1905.*x.2+6.4102.*x -8.2574;fy=abs(f-y); fy2=fy.2; Ew=max(fy), E1=sum(fy)/n, E2=sqrt(sum(fy2)/n)plot(xi,y,'r*'), hold on, plot(x,F