1、2011級工科碩士研究生矩陣與數值分析課程數值實驗2011級工科碩士研究生矩陣與數值分析課程數值實驗一、 對于數列，有如下兩種生成方式1、首項為，遞推公式為；2、前兩項為，遞推公式為；給出利用上述兩種遞推公式生成的序列的第50項。第一種方法：clc;clear;a=1;for i=2:50 a=1/3*a; i=i+1;enddisp('第50項為:')disp(a)結果如下：第50項為: 4.1789e-024第二種方法：clc;clear;a=1;b=1/3;for i=3:50 c=10/3*b-a; a=b; b=c; i=i+1;enddisp('第50項為:

2、')disp(c)結果如下：第50項為: 2.0604e+006分析：由上述結果可知，第一種遞推公式計算結果是精確的，而第二種遞推公式計算的結果前幾項比較接近，后來計算的誤差就越來較大了。二、 利用迭代格式及Aitken加速后的新迭代格式求方程在內的根Aitken加速后：clc;clear;format long a=1;b=sqrt(10/(a+4);c=sqrt(10/(b+4);x=c-(c-b)2/(a-2*b+c);i=0;while(x3+4*x2-10=0) d=sqrt(10/(c+4); a=b; b=c; c=d; x=c-(c-b)2/(a-2*b+c); i=i

3、+1;enddisp('此方程的根為：')disp(x)disp('迭代次數為：')disp(i)結果如下：此方程的根為： 1.36523001341410迭代次數為： 7未加速：clc;clear;format long x=1; i=0;while(x3+4*x2-10=0) x=sqrt(10/(x+4); i=i+1;enddisp('此方程的根為：')disp(x)disp('迭代次數為：')disp(i)結果如下：此方程的根為： 1.36523001341410迭代次數為： 18三、解線性方程組 1分別Jacobi迭代

4、法和Gauss-Seidel迭代法求解線性方程組迭代法計算停止的條件為：disp('此線性方程組的解為：x=')disp(y')disp('迭代次數為：')disp(k)結果如下：此線性方程組的解為：x= 0.435 1.680 0.24463224043774 -0.575迭代次數為： 15Jacobi迭代法：clc;clear;format long A=6 2 1 -2;2 5 0 -2; -2 0 8 5;1 3 2 7;b=4 7 -1 0;delta=10(-6);%誤差n=length(b);k=1;x=zeros(n,1);y=zeros

5、(n,1);while 1 for i=1:n y(i)=b(i); for j=1:n y(i)=y(i)-A(i,j)*x(j); end y(i)=x(i)+y(i)/A(i,i); end if norm(y-x,inf)<delta break; end x=y; k=k+1;end Gauss-Seidel迭代法：clc;clear;format long A=6 2 1 -2;2 5 0 -2; -2 0 8 5;1 3 2 7;b=4 7 -1 0;delta=10(-6);%誤差n=length(b);k=1;x=zeros(n,1);y=zeros(n,1);whil

6、e 1 for i=1:n y(i)=b(i); for j=1:i-1 y(i)=y(i)-A(i,j)*y(j); end for j=i:n y(i)=y(i)-A(i,j)*x(j); end y(i)=x(i)+y(i)/A(i,i); end if norm(y-x,inf)<delta break; end x=y; k=k+1;enddisp('此線性方程組的解為：x=')disp(y)disp('迭代次數為：')disp(k)結果如下：此線性方程組的解為：x= 0.111 1.986 0.24463241176980 -0.578迭代次數

7、為： 8 分析：Jacobi迭代法稱簡單迭代法，程序編譯比較簡單，但是在迭代過程中，對已經算出來的信息未加充分利用。Gauss-Seidel迭代法占用內存小，收斂快。Jacobi迭代法迭代次數為15，而Gauss-Seidel迭代法次數為8，迭代的次數更少，說明Gauss-Seidel迭代法收斂性更好。2. 用Gauss列主元消去法、QR方法求解如下方程組：for i=n-1:-1:1 sum=0; for j=i+1:n sum=sum+x(j,1)*B(i,j); end x(i,1)=(B(i,n+1)-sum)/B(i,i);enddisp('此方程組的解為：x=')d

8、isp(x)結果如下：此方程組的解為：x= 1.5417 -2.7500 0.0833 1.6667Gauss列主元消去法：clc;clear;A=2 2 1 2;4 1 3 -1; -4 -2 0 1;2 3 2 3;b=1;2;1;0;n=length(b); x=zeros(n,1);B=A b;for k=1:n-1 max=k; for i=k+1:n if B(i,k)>B(max,k) max=i; end end temp=B(k,k:n+1); B(k,k:n+1)=B(max,k:n+1); B(max,k:n+1)=temp; for i=k+1:n B(i,k)=

9、-B(i,k)/B(k,k); B(i,k+1:n+1)=B(i,k+1:n+1)+B(i,k)*B(k,k+1:n+1); endendx(n,1)=B(n,n+1)/B(n,n); QR方法：clc;clear;a=2 2 1 2;4 1 3 -1; -4 -2 0 1;2 3 2 3;b=1;2;1;0;x=qrfenjie(a,b);disp('此方程組的解為：x=')disp(x)單獨m文件：function x=qrfenjie(b,c)h=b;q=zeros(max(size(b),max(size(b)*(max(size(b)-1);g=max(size(b)

10、;for i=1:max(size(b)-1 z=eye(size(b); w=b(:,1)-norm(b(:,1)*z(:,1); q(1:max(size(b),1+(i-1)*g:(i-1)*g+max(size(b)=z-2/(w'*w)*w*w' b=q(1:max(size(b),1+(i-1)*g:(i-1)*g+max(size(b)*b; b=b(2:max(size(b),2:max(size(b);endfor i=1:g-2q(1:g,1:g)=eye(i,i),zeros(i,g-i);zeros(g-i,i),q(1:g-i,i*g+1:i*g+g-

11、i)*q(1:g,1:g);endd=q(1:g,1:g)*h,c;x=zeros(1,g);a1=d(:,1:g);a2=d(:,g+1);x(g)=a2(g)/a1(g,g);for k = g-1:-1:1 x(k)=a2(k); for p=g:-1:k+1 x(k)= x(k)-a1(k,p)*x(p); end x(k)=x(k)/a1(k,k); end結果如下：此方程組的解為：x= 1.5417 -2.7500 0.0833 1.6667四、已知一組數據點，編寫程序求解三次樣條插值函數滿足 并針對下面一組具體實驗數據0.250.30.390.450.530.50000.5477

12、0.62450.67080.7280求解，其中邊界條件為.for i=2:n-1 A(i,i-1)=dx(i-1); A(i,i)=2*(dx(i-1)+dx(i); A(i,i+1)=dx(i); B(i,1)=3*(dy(i)/dx(i)-dy(i-1)/dx(i-1);end%端點一階導數條件 %if flag=1 A(1,1)=2*dx(1);A(1,2)=dx(1); A(n,n-1)=dx(n-1);A(n,n)=2*dx(n-1); B(1,1)=3*(dy(1)/dx(1)-vl); B(n,1)=3*(vr-dy(n-1)/dx(n-1);end% - %端點二階導數條件 %

13、if flag=2 A(1,1)=2;A(n,n)=2; B(1,1)=vl;B(n,1)=vr;end% - %c=AB;for i=1:n-1 d(i)=(c(i+1)-c(i)/(3*dx(i); b(i)=dy(i)/dx(i)-dx(i)*(2*c(i)+c(i+1)/3;endclc;clear;format short g;x1=0.25 0.3 0.39 0.45 0.53'y1=0.5000 0.5477 0.6245 0.6708 0.7280'u1=0;un=0;xx1=x1(1):0.001:x1(end)'yy1 b1 c1 d1=spline

14、3(x1,y1,xx1,2,u1,un);disp('各區間上的三次樣條插值系數為：')fprintf('ttb1tttc1ttttd1n');disp(b1 c1(1:end-1,1) d1);單獨m文件：functionyy b c d=spline3(x,y,xx,flag,vl,vr)%樣條函數為：Si(x)=yi+bi*(x-xi)+ci*(x-xi)2+di*(x-xi)3n=length(x);a=zeros(n-1,1);b=a;d=a;dx=a;dy=a;A=zeros(n);B=zeros(n,1);for i=1:n-1 a(i)=y(i)

15、; dx(i)=x(i+1)-x(i); dy(i)=y(i+1)-y(i);end mm nn=size(xx);yy=zeros(mm,nn);for i=1:mm*nn for ii=1:n-1 if xx(i)>=x(ii) & xx(i)<x(ii+1) j=ii; break; elseif xx(i)=x(n) j=n-1; end end yy(i)=a(j)+b(j)*(xx(i)-x(j)+c(j)*(xx(i)-x(j)2+d(j)*(xx(i)-x(j)3; endend結果如下：各區間上的三次樣條插值系數為：b1c1d1 0.96966 0 -6.

16、2652 0.92267 -0.93977 1.8813 0.79923 -0.43181 -0.46 0.74245 -0.51461 2.1442五、編寫程序構造區間上的以等分結點為插值結點的Newton插值公式，假設結點數為（包括兩個端點），給定相應的函數值，插值區間和等分的份數，該程序能快速計算出相應的插值公式。以，為例計算其對應的插值公式，分別取不同的值并畫出原函數的圖像以及插值函數的圖像，觀察當增大時的逼近效果.Newton插值:clc;clear;syms xa=-1;b=1;n=3;x0=a:2/(n-1):b;y0=1./(1+25.*x0.2);x1=a:0.0001:b;

17、y1=1./(1+25.*x1.2);plot(x1,y1)hold ony=newton(x0,y0);plot(x1,subs(y,x,x1)axis(-1,1,-1,1);單獨m文件：function y=newton(x0,y0)syms x Bn=length(x0);A=zeros(n,n);A(:,1)=y0' for i=2:n for j=i:n A(j,i)=(A(j,i-1)-A(j-1,i-1)/(x0(j)-x0(j-i+1); endendB(1,1)=1;for i=2:n B(i,1)=B(i-1,1)*(x-x0(i-1);endfor i=1:n C

18、(i)=A(i,i).*B(i,1);endp=0;for i=1:n p=p+C(i);endP=expand(p);digits(6);y=vpa(P)n=3時插值結果如下：y = 1.-.961538*x2n=5時插值結果如下：y = 1.-4.27719*x2+3.31565*x4n=7時插值結果如下：y = 1.+.124928e-14*x-8.78409*x2-.284080e-14*x3+20.9574*x4-13.1349*x6n=10時插值結果如下：y= .138809e-14*x-8.26092*x2+.299257e-14*x3+30.7285*x4+.789492e-1

19、5*x7+21.6248*x8-.218329e-14*x5-44.9155*x6+.861538n=20時插值結果如下：y= .216415e-14*x+95604.8*x16-21.6235*x2+121019.*x12-58585.0*x10+.832282e-11*x9-.210106e-12*x3+327.726*x4+.390814e-11*x7+17172.5*x8+.647870e-11*x13-.921334e-12*x5-3055.32*x6+.459534e-11*x17-25671.2*x18-146791.*x14-.292343e-10*x15+.128399e-1

20、0*x11+.992681n=30時插值結果如下：y=-.577993e-8*x11+.242379e-15*x+.999614-.600443e-9*x9-.122518e7*x10-24.5956*x2-.144197e-12*x3+551.257*x4-.387656e-11*x5-9924.48*x6-.820216e-10*x7+131063.*x8+.806398e7*x12-.903419e-8*x13-.375301e8*x14+.468529e-7*x15+.123950e9*x16-.614169e-7*x17-.289894e9*x18+.339302e-7*x19+.474656e9*x20+.134317e-6*x21-.530053e9*x22-.781817e-7*x23+.383368e9*x24-.258139e-7*x25-.161437e9*x26+.411067e-8*x27+.299799e8*x28n=50時插值結果如下：y= .163155e16*x30+.288123e-14*x+.132963e9*x12-.747119e7*x10+.908504e-5*x13+.999999-.212026e-7*x9+357069.*x8+.228985e11*x16+.833425e-4*x17-.570912e1