數(shù)值分析上機作業(yè)1-1

1、數(shù)值計算方法上機題目11、實驗1. 病態(tài)問題實驗?zāi)康模核惴ㄓ小皟?yōu)”與“劣”之分，問題也有“好”和“壞”之別。所謂壞問題就是問題本身的解對數(shù)據(jù)變化的比較敏感，反之屬于好問題。希望讀者通過本實驗對此有一個初步的體會。數(shù)值分析的大部分研究課題中，如線性代數(shù)方程組、矩陣特征值問題、非線性方程及方程組等都存在病態(tài)的問題。病態(tài)問題要通過研究和構(gòu)造特殊的算法來解決，當(dāng)然一般要付出一些代價（如耗用更多的機器時間、占用更多的存儲空間等）。問題提出：考慮一個高次的代數(shù)多項式 （E1-1）顯然該多項式的全部根為l，2，20，共計20個，且每個根都是單重的（也稱為簡單的）?，F(xiàn)考慮該多項式方程的一個擾動 (E1-2)其

2、中是一個非常小的數(shù)。這相當(dāng)于是對（E1-1）中的系數(shù)作一個小的擾動。我們希望比較（E1-1）和（E1-2）根的差別，從而分析方程（E1-1）的解對擾動的敏感性。實驗內(nèi)容：為了實現(xiàn)方便，我們先介紹兩個 Matlab函數(shù)：“roots”和“poly”，輸入函數(shù)uroots（a）其中若變量存儲維的向量，則該函數(shù)的輸出為一個維的向量。設(shè)a的元素依次為，則輸出u的各分量是多項式方程的全部根，而函數(shù)b=poly(v)的輸出b是一個n1維變量，它是以n維變量v的各分量為根的多項式的系數(shù)?？梢姟皉oots”和“Poly”是兩個互逆的運算函數(shù). ve=zeros(1,21); ve(2)=ess; roots(

3、poly(1:20)+ve)上述簡單的Matlab程序便得到（E1-2）的全部根，程序中的“ess”即是（E1-2）中的。實驗要求：（1）選擇充分小的ess，反復(fù)進行上述實驗，記錄結(jié)果的變化并分析它們。如果擾動項的系數(shù)很小，我們自然感覺（E1-1）和(E1-2)的解應(yīng)當(dāng)相差很小。計算中你有什么出乎意料的發(fā)現(xiàn)？表明有些解關(guān)于如此的擾動敏感性如何？（2）將方程（E1-2）中的擾動項改成或其他形式，實驗中又有怎樣的現(xiàn)象出現(xiàn)？實驗步驟：（1）程序function t_charpt1_1clcresult=inputdlg('請輸入擾動項:在0 20之間的整數(shù):','charpt

4、 1_1',1,'19');Numb=str2num(char(result);if(Numb>20)|(Numb<0)errordlg('請輸入正確的擾動項:0 20之間的整數(shù)!');return;end result=inputdlg('請輸入(0 1)之間的擾動常數(shù):','charpt 1_1',1,'0.00001');ess=str2num(char(result);ve=zeros(1,21);ve(21-Numb)=ess;root=roots(poly(1:20)+ve);x0=

5、real(root); y0=imag(root);plot(x0',y0', '*');disp('對擾動項 ',num2str(Numb),'加擾動',num2str(ess),'得到的全部根為:');disp(num2str(root);2、 實驗結(jié)果分析 ess分別為1e-6,1e-8.1e-10,1e-12.對擾動項 19加擾動1e-006得到的全部根為:21.3025+1.56717i 21.3025-1.56717i 18.5028+3.6004i 18.5028-3.6004i 15.1651+3.

6、76125i 15.1651-3.76125i 12.4866+2.88278i 12.4866-2.88278i 10.5225+1.71959i 10.5225-1.71959i 9.04485+0.594589i 9.04485-0.594589i 7.9489+0i 7.00247+0i 5.99995+0i 5+0i 4+0i 3+0i 2+0i 1+0i 對擾動項 19加擾動1e-010得到的全部根為:19.9953+0i 19.0323+0i 17.8696+0i 17.2186+0i 15.4988+0.0211828i 15.4988-0.0211828i 13.7707+0

7、i 13.1598+0i 11.9343+0i 11.029+0i 9.99073+0i 9.00247+0i 7.99952+0i 7.00007+0i 5.99999+0i 5+0i 4+0i 3+0i 2+0i 1+0i ess分別為1e-6,1e-8.1e-10,1e-12的圖像如下： 從實驗的圖形中可以看出，當(dāng)ess充分小時，方程E.1.1和方程E.1.2的解相差很小，當(dāng)ess逐漸增大時，方程的解就出現(xiàn)了病態(tài)解，這些解都呈現(xiàn)復(fù)共軛性質(zhì)。（2） 將擾動項加到x18上后，ess=1e-009時方程的解都比較準確，沒有出現(xiàn)復(fù)共軛現(xiàn)象。ess=1e-008時誤差與x19(ess=1e-009

8、)時相當(dāng)，即擾動加到x18上比加到x19小一個數(shù)量級。對x8的擾動ess=1000時沒有出現(xiàn)復(fù)共軛，誤差很小；對x的擾動ess=10e10時沒有出現(xiàn)復(fù)共軛，誤差很小。因此，擾動作用到xn上時，n越小，擾動引起的誤差越小。2、實驗2。多項式插值的振蕩現(xiàn)象，即插值的龍格（Runge）現(xiàn)象問題提出：考慮在一個固定的區(qū)間上用插值逼近一個函數(shù)。顯然，拉格朗日插值中使用的節(jié)點越多，插值多項式的次數(shù)就越高、自然關(guān)心插值多項式的次數(shù)增加時，是否也更加靠近被逼近的函數(shù)。龍格給出的一個例子是極著名并富有啟發(fā)性的。設(shè)區(qū)間上函數(shù)實驗內(nèi)容：考慮區(qū)間的一個等距劃分，分點為則拉格朗日插值多項式為其中的，是n次拉格朗日插值基

9、函數(shù)。實驗要求：（l）選擇不斷增大的分點數(shù)目，畫出原函數(shù)及插值多項式函數(shù)在上的圖像，比較并分析實驗結(jié)果。（2）選擇其他的函數(shù)，例如定義在區(qū)間-5，5上的函數(shù) 重復(fù)上述的實驗看其結(jié)果如何。（3）區(qū)間上切比雪夫點的定義為 以為插值節(jié)點構(gòu)造上述各函數(shù)的拉格朗日插值多項式，比較其結(jié)果。實驗步驟：（1） 試驗程序：function y=Lagrange(x0, y0, x);% Lagrange插值 n= length(x0); m=length(x);for i=1:m z=x(i); s=0.0; for k=1:n p=1.0; for j=1:n if(j = k) p = p*(z - x0(

10、j)/(x0(k) - x0(j); end end s = s + p*y0(k); end y(i) = s;endfunction t_charpt2 promps = '請選擇實驗函數(shù)，若選f(x),請輸入f,若選h(x),請輸入h,若選g(x),請輸入g:'titles = 'charpt_2'result = inputdlg(promps,'charpt 2',1,'f');Nb_f = char(result);if(Nb_f = 'f' & Nb_f = 'h' &

11、 Nb_f = 'g')errordlg('實驗函數(shù)選擇錯誤！');return;end result = inputdlg('請輸入插值結(jié)點數(shù)N:','charpt_2',1,'10');Nd = str2num(char(result);if(Nd <1)errordlg('結(jié)點輸入錯誤！');return;end switch Nb_f case 'f' f=inline('1./(1+25*x.2)'); a = -1;b = 1; case '

12、h' f=inline('x./(1+x.4)'); a = -5; b = 5; case 'g' f=inline('atan(x)'); a = -5; b= 5;endx0 = linspace(a, b, Nd+1); y0 = feval(f, x0); x = a:0.1:b; y = Lagrange(x0, y0, x); fplot(f, a b, 'co'); hold on; plot(x, y, 'b-');xlabel('x'); ylabel('y =

13、f(x) o and y = Ln(x)-'); 增大分點n=2，3，時，拉格朗日插值函數(shù)曲線如圖所示。 n=3 n=6 n=7 n=8 從圖中可以看出，隨著n的增大，拉格朗日插值函數(shù)在x=0附近較好地逼近了原來的函數(shù)f(x)，但是卻在兩端x= -1和x=1處出現(xiàn)了很大的振蕩現(xiàn)象。通過分析圖形，可以看出，當(dāng)n為奇數(shù)時，雖然有振蕩，但振蕩的幅度不算太大，n為偶數(shù)時，其振蕩幅度變得很大。（2） 將原來的f(x)換為其他函數(shù)如h(x)、g(x)，結(jié)果如圖所示。其中h(x), g(x)均定義在-5，5區(qū)間上，h(x)=x/(1+x4)，g(x)=arctan x。h(x), n=7 h(x),

14、 n=8h(x), n=9 h(x), n=10g(x), n=8 g(x), n=9g(x), n=12 g(x), n=13分析兩個函數(shù)的插值圖形，可以看出：隨著n的增大，拉格朗日插值函數(shù)在x=0附近較好地逼近了原來的函數(shù)f(x)，但是卻在兩端x= -5和x=5處出現(xiàn)了很大的振蕩現(xiàn)象。通過圖形可以看出，當(dāng)n為偶數(shù)時，雖然有振蕩，但振蕩的幅度不算太大，n為奇數(shù)時，其振蕩幅度變得很大。原因和上面f(x)的插值類似，h(x)、g(x)本身是奇函數(shù)，如果n為偶數(shù)，那么Lagrange插值函數(shù)Ln(x)的最高次項xn-1是奇次冪，比較符合h(x)、g(x)本身是奇函數(shù)的性質(zhì)；如果n為奇數(shù)，那么Lag

15、range插值函數(shù)Ln(x)的最高次項xn-1是偶次冪，與h(x)、g(x)本身是奇函數(shù)的性質(zhì)相反，因此振蕩可能更劇烈。3、實驗3。 樣條插值的收斂性問題提出：一般的多項式插值不能保證收斂性，即插值的節(jié)點多，效果不一定就好。對樣條函數(shù)插值又如何呢？理論上證明樣條插值的收斂性是比較困難的，也超出了本課程的內(nèi)容。通過本實驗可以驗證這一理論結(jié)果。實驗內(nèi)容：請按一定的規(guī)則分別選擇等距或者非等距的插值節(jié)點，并不斷增加插值節(jié)點的個數(shù)?？紤]實驗2.中的函數(shù)或選擇其它你有興趣的函數(shù)，可以用 Matab的函數(shù) “spline”作此函數(shù)的三次樣條插值。在較新版本的Matlab中，還提供有spline工具箱，你可以

16、找到極豐富的樣條工具，包括B-樣條。實驗要求：（1）隨節(jié)點個數(shù)增加，比較被逼近函數(shù)和樣條插值函數(shù)誤差變化情況。分析所得結(jié)果并與拉格朗目多項式插值比較。（2）樣條插值的思想最早產(chǎn)生于工業(yè)部門。作為工業(yè)應(yīng)用的例子，考慮如下問題：某汽車制造商用三次樣條插值設(shè)計車門的曲線，其中一段的數(shù)據(jù)如下：0123456789100.00.791.532.192.713.033.272.893.063.193.290.80.2實驗步驟：（1）程序：function t_charpt2promps = '請選擇實驗函數(shù)，若選f(x),請輸入f,若選h(x),請輸入h,若選g(x),請輸入g:'titl

17、es = 'charpt_2'result = inputdlg(promps,'charpt 2',1,'f');Nb_f = char(result);if(Nb_f = 'f' & Nb_f = 'h' & Nb_f = 'g')errordlg('實驗函數(shù)選擇錯誤！');return;end result = inputdlg('請輸入插值結(jié)點數(shù)N:','charpt_2',1,'10');Nd = str2nu

18、m(char(result);if(Nd <1)errordlg('結(jié)點輸入錯誤！');return;end switch Nb_f case 'f' f=inline('1./(1+25*x.2)'); a = -1;b = 1; case 'h' f=inline('x./(1+x.4)'); a = -5; b = 5; case 'g' f=inline('atan(x)'); a = -5; b= 5;end x0 = linspace(a, b, Nd+1); y0 = feval(f,