數(shù)值分析上機(jī)作業(yè)1 1解析

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

2、，且每個根都是單重的（也稱為簡單的）。現(xiàn)考慮該多項式方程的一個擾動 19?x)?0xp( (E1-2) 19?x的系數(shù)作一個小的擾動。我們希望是一個非常小的數(shù)。這相當(dāng)于是對（E1-1其中）中 ）的解對擾動的敏感性。E1-2）根的差別，從而分析方程（E1-1比較（E1-1）和（ ：實驗內(nèi)容 ”，輸入函數(shù)roots”和“poly為了實現(xiàn)方便，我們先介紹兩個 Matlab函數(shù)：“ ）（aurootsnua1?n的元素依次a維的向量，則該函數(shù)的輸出其中若變量維的向量。設(shè)存儲為一個a,.,aa u，則輸出為的各分量是多項式方程1?2n1nn?1?.axx?axa?0a 1n12n?的全部根，而函數(shù) b=

3、poly(v) 的輸出b是一個n1維變量，它是以n維變量v的各分量為根的多項式的系數(shù)?？梢姟皉oots”和“Poly”是兩個互逆的運算函數(shù). ve=zeros(1,21); ve(2)=ess; roots(poly(1:20)+ve) ?。中的） ess”即是（E1-2程序中的（上述簡單的Matlab程序便得到E1-2）的全部根，“實驗要求： （1）選擇充分小的ess，反復(fù)進(jìn)行上述實驗，記錄結(jié)果的變化并分析它們。如果擾動項?很小，我們自然感覺（E1-1）和(E1-2)的系數(shù)的解應(yīng)當(dāng)相差很小。計算中你有什么出乎意料的發(fā)現(xiàn)？表明有些解關(guān)于如此的擾動敏感性如何？ 18?x E1-22（）將方程（）

4、中的擾動項改成或其他形式，實驗中又有怎樣的現(xiàn)象出現(xiàn)？1 實驗步驟： （1）程序 t_charpt1_1function clcresult=inputdlg(請輸入擾動項:在0 20之間的整數(shù):,charpt );191_1,1, Numb=str2num(char(result);if(Numb20)|(Numb0)errordlg(請輸入正確的擾動項:0 20之間的整 end;!);return數(shù) result=inputdlg(請輸入(0 1)之間的擾動常數(shù):,charpt );0.000011_1,1, ess=str2num(char(result); ve=zeros(1,21);

5、 ve(21-Numb)=ess; root=roots(poly(1:20)+ve); x0=real(root); y0=imag(root); );*plot(x0,y0, disp(對擾動項 ,num2str(Numb),加擾動,num2str(ess),得到的全部根 );:為 disp(num2str(root);二、實驗結(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+

6、3.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.770

7、7+0i 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的圖像如下： 2 當(dāng)?shù)慕庀嗖詈苄?，E.1.1和方程E.1.2 從實驗的圖形中可以看出，當(dāng)ess充分小時，方程 逐漸增大時，方程的解就出現(xiàn)了病態(tài)解，這些解都呈現(xiàn)復(fù)共軛性質(zhì)。ess18時方程的解都比較準(zhǔn)確，沒有出現(xiàn)復(fù)共軛現(xiàn)象。xess=1e-009上后，（2） 將擾動項加到191819小一個數(shù)量級。對上比加到時誤差

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

9、一個等距劃分，分點為考慮區(qū)間i2n,2,.,1,i?0,x?1? in 則拉格朗日插值多項式為n1?)(x)?xlL( in2x?2510i?i)xl(n,.,1,20i?, n其中的是次拉格朗日插值基函數(shù)。，i 實驗要求：)(Lx)xf2n?,3,.(在及插值多項式函數(shù)，畫出原函數(shù)）（l選擇不斷增大的分點數(shù)目n1,?1 上的圖像，比較并分析實驗結(jié)果。 5-52（）選擇其他的函數(shù)，例如定義在區(qū)間，上的函數(shù)3 x,g(x)?arctanh(x)?x 4x1? 重復(fù)上述的實驗看其結(jié)果如何。a,b上切比雪夫點的定義為 （3）區(qū)間?b?abk?1)?a(2?,k?1,2,.,?xn?cos1 ? k2

10、)2(2(n?1?x,x,.,x為插值節(jié)點構(gòu)造上述各函數(shù)的拉格朗日插值多項式，比較其結(jié)果。 以12n?1實驗步驟： （1） 試驗程序： y=Lagrange(x0, y0, x);function % Lagrange插值 n= length(x0); m=length(x); i=1:mfor z=x(i); s=0.0; k=1:n for p=1.0; j=1:n for (j = k) if p = p*(z - x0(j)/(x0(k) - x0(j); end end s = s + p*y0(k); end y(i) = s;end t_charpt2functionpromps

11、 = 請選擇實驗函數(shù)，若選f(x),請輸入f,若選h(x),請輸入h,若選g(x),請輸 ;入g: ;titles = charpt_2 );,1,fresult = inputdlg(promps,charpt 2 Nb_f = char(result);if(Nb_f = f & Nb_f = h & Nb_f = g)errordlg(實驗函數(shù)選擇錯誤！ end;);return );10,N:charpt_2,1,result = inputdlg(請輸入插值結(jié)點數(shù) Nd = str2num(char(result); end;結(jié)點輸入錯誤！);returnif(Nd 1)errord

12、lg( Nb_fswitch4 f case ); a = -1;b = 1;1./(1+25*x.2) f=inline( h case ); a = -5; b = 5;x./(1+x.4) f=inline( g case ); a = -5; b= 5;atan(x) f=inline( end x0 = linspace(a, b, Nd+1); y0 = feval(f, x0); x = a:0.1:b; y = Lagrange(x0, y0, x); );co fplot(f, a b, ;on hold );b- plot(x, y, 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，4)，g(x)=arctan xh(x)=x/(1+x5區(qū)間上，。 5 h(x), n=8 h(x

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

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

16、有spline工具箱，你可以找到極豐富的樣條工具，包括B-樣條。 實驗要求： （1）隨節(jié)點個數(shù)增加，比較被逼近函數(shù)和樣條插值函數(shù)誤差變化情況。分析所得結(jié)果并與拉格朗目多項式插值比較。 （2）樣條插值的思想最早產(chǎn)生于工業(yè)部門。作為工業(yè)應(yīng)用的例子，考慮如下問題：某汽車制造商用三次樣條插值設(shè)計車門的曲線，其中一段的數(shù)據(jù)如下： x k0 1 2 3 4 5 6 7 8 9 10 y k0.0 0.79 1.53 2.19 2.71 3.03 3.27 2.89 3.06 3.19 3.29 y k0.8 0.2 實驗步驟： （1）程序： t_charpt2functionpromps = 請選擇實驗函

17、數(shù)，若選f(x),請輸入f,若選h(x),請輸入h,若選g(x),請輸 ;g:入 ;charpt_2titles = );fcharpt 2,1,result = inputdlg(promps, Nb_f = char(result);if(Nb_f = f & Nb_f = h & Nb_f = g)errordlg(實驗函數(shù)選擇錯誤！ endreturn;); );10,charpt_2,1,N:result = inputdlg(請輸入插值結(jié)點數(shù) Nd = str2num(char(result); end;);returnif(Nd 1)errordlg(結(jié)點輸入錯誤！ Nb_fsw

18、itch f case ); a = -1;b = 1;1./(1+25*x.2) f=inline( hcase ); a = -5; b = 5;x./(1+x.4) f=inline( g case7 ); a = -5; b= 5; f=inline(atan(x) end x0 = linspace(a, b, Nd+1); y0 = feval(f, x0); x = a:0.1:b; cs = spline(x0, y0); y = ppval(cs, x); ); plot(x, y, k-o); hold on plot(x0, y0, );); ylabel(y = f(x) o and y = Spline(x)- xlabel(x實驗結(jié)果： 如圖所示。 f(x), n=5 n=10 n=20 h(x), n=5 h(x), n=10 n=20 g(x), n=5 n=10 n=20 圖中可以看出