《矩陣與數(shù)值分析》課程數(shù)值實(shí)驗(yàn)大作業(yè)

1、優(yōu)質(zhì)文本2011級工科碩士研究生?矩陣與數(shù)值分析?課程數(shù)值實(shí)驗(yàn)班 級： 學(xué) 號(hào)： 姓 名： 任課教師： 大連理工大學(xué)2011年12月20日優(yōu)質(zhì)文本一、 對于數(shù)列，有如下兩種生成方式1、首項(xiàng)為，遞推公式為；2、前兩項(xiàng)為，遞推公式為；給出利用上述兩種遞推公式生成的序列的第50項(xiàng)?！景吹谝环N遞推公式】clearclca=1;for i=1:50-1 a=a a(i)/3;enddisp('數(shù)列第50項(xiàng)小數(shù)表達(dá)為：')format longdisp(a(50)disp('分?jǐn)?shù)表達(dá)為：')format ratdisp(a(50)format short運(yùn)行結(jié)果數(shù)列第50項(xiàng)

2、小數(shù)表達(dá)為： 4.178866707295615e-024分?jǐn)?shù)表達(dá)為： 1/239299329230617530000000【按第二種遞推公式】clearclca=1 1/3;for i=2:50-1 a=a 10/3*a(i)-a(i-1);endformat ratdisp('數(shù)列第50項(xiàng)為：')disp(a(50)format short運(yùn)行結(jié)果數(shù)列第50項(xiàng)為：2060436 【分析】第一種算法數(shù)值穩(wěn)定，計(jì)算過程舍入誤差被嚴(yán)格控制，且按1/3的公差不斷縮小。但第二種算法數(shù)值不穩(wěn)定。另外，在第二種算法中，假設(shè)將遞推公式“a=a 10/3*a(i)-a(i-1)中的分母移動(dòng)位

3、置，改寫成“a=a 10*a(i)/3-a(i-1)，那么程序運(yùn)行結(jié)果為-4966040，可以舍入誤差被放大的十分嚴(yán)重。二、 利用迭代格式及Aitken加速后的新迭代格式求方程在內(nèi)的根【未經(jīng)加速的代碼】clceps=1e-15;i=1;x0=1;format longwhile i<100 x1=sqrt(10/(x0+4); if abs(x1-x0)<=eps break end x0=x1; i=i+1;enddisp('方程的解精度10(-15)')disp(x1)disp('未經(jīng)加速的迭代次數(shù)')disp(i)運(yùn)行結(jié)果方程的解精度10(-1

4、5) 1.36523001341410未經(jīng)加速的迭代次數(shù) 18【經(jīng)Aitken加速的代碼】clceps=1e-15;i=1;x0=1;format longwhile i<100 x1=sqrt(10/(x0+4); y=sqrt(10/(x1+4); z=sqrt(10/(y+4); x1=z-(z-y)2/(z-2*y+x1); if abs(x1-x0)<=eps break end x0=x1; i=i+1;enddisp('方程的解精度10(-15)')disp(x1)disp('未經(jīng)加速的迭代次數(shù)')disp(i)運(yùn)行結(jié)果方程的解精度10

5、(-15) 1.36523001341410未經(jīng)加速的迭代次數(shù) 3【分析】Aitken加速能對數(shù)列xk起明顯的加速作用，在要求相同方程解精度的條件下，它能將迭代次數(shù)顯著降低。實(shí)際上，Aitken有時(shí)甚至能將發(fā)散的數(shù)列加速后變?yōu)槭諗?。三、解線性方程組 1分別Jacobi迭代法和Gauss-Seidel迭代法求解線性方程組迭代法計(jì)算停止的條件為： 2. 用Gauss列主元消去法、QR方法求解如下方程組：【1. Jacobi方法】clci=1;eps=1e-6;A= 6 2 1 -2; 2 5 0 -2; -2 0 8 5; 1 3 2 7;b=4 7 -1 0'x0=zeros(4,1);

6、D=diag(diag(A);L=-tril(A,-1);U=-triu(A,1);B=inv(D)*(L+U);f=inv(D)*b;while i<100 x1=B*x0+f; if norm(x1-x0)<=eps break end x0=x1; i=i+1;enddisp('方程的解精度10(-6)')disp(x1)disp('迭代次數(shù)')disp(i)運(yùn)行結(jié)果方程的解精度10(-6) 0.05204951386229 1.15094332647528 0.24463239101199 -0.57059207123432迭代次數(shù) 16【1

7、. Gauss-Seidel方法】clci=1;eps=1e-6;A= 6 2 1 -2; 2 5 0 -2; -2 0 8 5; 1 3 2 7;b=4 7 -1 0'x0=zeros(4,1);D=diag(diag(A);L=-tril(A,-1);U=-triu(A,1);B=inv(D-L)*(U);f=inv(D-L)*b;while i<100 x1=B*x0+f; if norm(x1-x0)<=eps break end x0=x1; i=i+1;enddisp('方程的解精度10(-6)')disp(x1)disp('迭代次數(shù)&#

8、39;)disp(i)運(yùn)行結(jié)果方程的解精度10(-6) 0.05204946628111 1.15094338455986 0.24463241176981 -0.57059206335719迭代次數(shù) 8【2. Gauss列主元消去法】clcA= 2 2 1 2; 4 1 3 -1; -4 -2 0 1; 2 3 2 3;b=1 2 1 0'Ab=A b;N,N=size(A);x=zeros(N,1);for p=1:N-1 max1,j=max(abs(A(p:N,p); temp=Ab(p,:); Ab(p,:)=Ab(j+p-1,:); Ab(j+p-1,:)=temp; if

9、 Ab(p,p)=0 disp('方程無解') break end for k=p+1:N mult=Ab(k,p)/Ab(p,p); Ab(k,p:N+1)=Ab(k,p:N+1)-mult*Ab(p,p:N+1); endendb=Ab(:,N+1);xx=zeros(N,1);for k=N:-1:1 xx(k)=b(k)/Ab(k,k); i=(1:k-1)' b(i)=b(i)-xx(k)*Ab(i,k);endx=xx運(yùn)行結(jié)果x = 1.54166666666667 -2.75000000000000 0.08333333333333 1.666666666

10、66667【2. QR分解法】clcfor i=1:3 Ai=zeros(5-i); Qi=eye(4);endx=zeros(4,1);y=zeros(4,1);R=zeros(4);A1=2 2 1 2; 4 1 3 -1; -4 -2 0 1; 2 3 2 3;b=1,2,1,0'for i=1:3 I=eye(size(Ai); a=Ai(:,1); w=a-norm(a)*I(:,1); hw=I-2/(w'*w)*(w*w'); Qi(i:4,i:4)=hw; if i=3 break end c=hw*Ai; Ai+1=c(2:max(size(Ai),2

11、:max(size(Ai);endQz=(Q3*Q2*Q1)'R=Q3*Q2*Q1*A1;c=Qz'*b;x(4)=c(4)/R(4,4);for i=3:-1:1 x(i)=(c(i)-R(i,i+1:4)*x(i+1:4)/R(i,i);endx運(yùn)行結(jié)果x = 1.54166666666666 -2.75000000000000 0.08333333333333 1.66666666666666【分析】Jacobi迭代法和Gauss-Seidel迭代法都可用來求解線性方程組。在同等精度下，求解本道題Jacobi迭代法迭代了16次，而Gauss-Seidel僅為8次，是前者的

12、一半。所以可以認(rèn)為，在某些情況下，Gauss-Seidel是比擬好的解法，但不總?cè)绱?。Gauss消去法可能發(fā)生小主元做除數(shù)，從而導(dǎo)致計(jì)算結(jié)果嚴(yán)重偏離真實(shí)值，所以在消元過程中，每一步都按列選主元，也就是Gauss列主元消去法，它可以有效防止過于嚴(yán)重的舍入誤差。正交矩陣是性態(tài)最好的矩陣，它不改變矩陣的條件數(shù)，通過QR分解來計(jì)算線性方程組，也可以防止誤差的放大，保證計(jì)算過程具有數(shù)值穩(wěn)定性。四、一組數(shù)據(jù)點(diǎn)，編寫一程序求解三次樣條插值函數(shù)滿足 并針對下面一組具體實(shí)驗(yàn)數(shù)據(jù)0.250.30.390.450.530.50000.54770.62450.67080.7280求解，其中邊界條件為.【程序代碼，文件

13、名 Spline.m】function s=Spline(x,y,f0,fn) % 輸入實(shí)驗(yàn)數(shù)據(jù)x,y；邊界二階導(dǎo)數(shù)f0,fnclcfigure(1)plot(x,y,'*r')hold onN=max(size(x);syms t s;for k=1:N-1; h(k)=x(k+1)-x(k);endg(1)=3*(y(2)-y(1)/h(1)-h(1)/2*f0;g(N)=3*(y(N)-y(N-1)/h(N-1)+h(N-1)/2*fn;for k=2:N-1 lamda(k)=h(k)/(h(k)+h(k-1); miu(k)=h(k-1)/(h(k)+h(k-1);g

14、(k)=3*(miu(k)*(y(k+1)-y(k)/h(k)+.lamda(k)*(y(k)-y(k-1)/h(k-1);endc=1,miu(2:N-1);b=2*ones(1,N);a=lamda(2:N-1),1;A=diag(c,1)+diag(b,0)+diag(a,-1);d=c;a=0,a;u(1)=b(1);for i=2:N l(i)=a(i)/u(i-1); u(i)=b(i)-l(i)*c(i-1);endL=eye(N)+diag(l(2:N),-1);U=diag(u)+diag(d,1);z(1)=g(1);for i=2:N z(i)=g(i)-l(i)*z(i

15、-1);endm(N)=z(N)/u(N);for i=N-1:-1:1 m(i)=(z(i)-c(i)*m(i+1)/u(i);endfor i=1:N-1 s(i)=(h(i)+2*(t-x(i)*(t-x(i+1)2*y(i)/(h(i)3+.(h(i)-2*(t-x(i+1)*(t-x(i)2*y(i+1)/(h(i)3+.(t-x(i)*(t-x(i+1)2*m(i)/(h(i)2+.(t-x(i+1)*(t-x(i)2*m(i+1)/(h(i)2;Enddigits(8)s=vpa(expand(s) % 輸出分段表達(dá)式for i=1:N-1 % 繪圖 v=x(i):0.005:x

16、(i+1); y=subs(s(i),v); plot(v,y); hold on;endgrid on命令窗口輸入x=0.25 0.3 0.39 0.45 0.53;y=0.5000,0.5477,0.6245,0.6708,0.7280;f0=0;fn=0;Spline(x,y,f0,fn)運(yùn)行結(jié)果ans =4.6988737*t2-.20505552*t+.35547747-6.2651650*t3,-2.6329843*t2+1.9945019*t+.13552173+1.8813439*t3,.10638708*t2+.92614706*t+.27440786-.45999912*t

17、3,-3.4093028*t2+2.5082075*t+.37098798e-1+2.1442156*t3【分析】運(yùn)行結(jié)果是一個(gè)四個(gè)元素的矩陣，各個(gè)元素依次代表四個(gè)分段區(qū)間上的三次樣條曲線。例如在第一個(gè)區(qū)間0.25 0.3上，擬合得到的三次樣條曲線方程4.6988737*t2-0.20505552*t+0.35547747-6.2651650*t3。由圖像可知，三次樣條曲線是很光滑的曲線，擬合效果很好。五、編寫程序構(gòu)造區(qū)間上的以等分結(jié)點(diǎn)為插值結(jié)點(diǎn)的Newton插值公式，假設(shè)結(jié)點(diǎn)數(shù)為包括兩個(gè)端點(diǎn)，給定相應(yīng)的函數(shù)值，插值區(qū)間和等分的份數(shù)，該程序能快速計(jì)算出相應(yīng)的插值公式。以，為例計(jì)算其對應(yīng)的插值公

18、式，分別取不同的值并畫出原函數(shù)的圖像以及插值函數(shù)的圖像，觀察當(dāng)增大時(shí)的逼近效果.【程序代碼，文件名 Newton.m】function f1=Newton(n)clcsyms x t fa=-1;b=1;f=1./(1+25*x.2);y=zeros(1,n);x=linspace(a,b,n); h=-1:0.01:1;m=n;c(1:m) = 0.0;yh=subs(f,h);y=subs(f,x); yk=y;f1=y(1); y1=0;k=1;for i=1:m-1 for j=i+1:m y1(j)=(y(j)-y(i)/(x(j)-x(i); end c(i)=y1(i+1); k=k*(t-x(i); f1=f1+c(i)*k; y=y1;endf1=simplify(f1) yfh=subs(f1,h);figure(1)plot(h,yfh,'-k',x,yk,'*r')grid onhold onx=-1:0.01:1;y=1./(1+25*x.2);plot(x,yh,'b')legend('插值曲線','插值節(jié)點(diǎn)','被插曲線');運(yùn)行結(jié)果N