計(jì)算方法 上機(jī)作業(yè)

1、目錄上機(jī)作業(yè)11上機(jī)作業(yè)25上機(jī)作業(yè)315上機(jī)作業(yè)417上機(jī)作業(yè)524上機(jī)作業(yè)627上機(jī)作業(yè)729計(jì)算方法上機(jī)作業(yè)：要求：1）抄題；2）迭代公式（初值）或簡(jiǎn)單原理；3）程序，結(jié)果；（打?。?）結(jié)果分析。上機(jī)作業(yè)11.分別用不動(dòng)點(diǎn)迭代與Newton法求解方程 2x-ex+3=0的正根與負(fù)根。解：1.不動(dòng)迭代法：迭代公式：pn+1=g(pn)，n=0,1,正根：選取的迭代格式為x(k+1)=ln(2x(k)+3)編寫(xiě)程序如下：x0=0;N=2000;Tol=10(-4);k=1;while k<=Nx=log(2*x0+3);if abs(x-x0)<Tol break; endk=k

2、+1;x0=x;enddisp(x);disp(k)結(jié)果：1.9239 10負(fù)根：選取的迭代格式為x(k+1)= exp(x(k)-3/2編寫(xiě)程序如下：x0=0;N=2000;Tol=10(-4);k=1;while k<=Nx=(exp(x0)-3)/2;if abs(x-x0)<Tol break; endk=k+1;x0=x;enddisp(x);disp(k)結(jié)果：-1.3734 72.Newton法：其迭代格式為x（k+1）=x（k）-f(x(k)/f，即x=(1-x)ex-3)/(2-ex)編寫(xiě)程序如下：正根：x0=3;N=2000;Tol=10(-7);for k=1

3、:N x=(1-x0)*exp(x0)-3)/(2-exp(x0); if abs(x-x0)<Tol break; end x0=x; enddisp(x);disp(k)結(jié)果：1.9239 6負(fù)根：x0=0;N=2000;Tol=10(-7);for k=1:N x=(1-x0)*exp(x0)-3)/(2-exp(x0); if abs(x-x0)<Tol break; end x0=x; enddisp(x);disp(k)結(jié)果：-1.3734 5結(jié)果分析：不動(dòng)點(diǎn)迭代法的原理比較簡(jiǎn)單。相比較之下，牛頓迭代法收斂于根的速度更快。2.用Newton法求解方程 x-sinx0 的

4、根。再用Steffensens method加速其收斂。解：Newton法迭代公式：，n=0,1,Steffensen加速原理：程序：Newton法：x0=-100, x1=100;N=2000;Tol=10(-9);for k=1:N x=x1-(x1-sin(x1)*(x1-x0)/(x1-sin(x1)-x0+sin(x0); if abs(x-x1)<Tol break; end x0=x1;x1=x; enddisp(x);disp(k)運(yùn)行結(jié)果：x0 = -100 1.4211e-014 2Steffensens method加速：M文件：function s=steff(f

5、,x0,eps,max)f=inline(f);x(1)=x0;disp('k x y z');for k=1:max y(k)=feval(f,x(k); z(k)=feval(f,y(k); x(k+1)=x(k)-(y(k)-x(k)2/(z(k)-2*y(k)+x(k); if abs(x(k+1)-x(k)<eps breakenddisp(sprintf('%d %f %f %f',k,x(k),y(k),z(k);ends=x(k+1);命令窗口：s=steff('x-sin(x)',1,1.0e-7,200)運(yùn)行結(jié)果如下：k

6、 x y z1 1.000000 0.158529 0.0006632 -0.035793 -0.000008 -0.000000s = 2.0680e-025結(jié)果分析：Steffensen加速極大地改善了原迭代的收斂性質(zhì)，它加速了牛頓法的收斂速度，且結(jié)果更為精確，但有些時(shí)候計(jì)算量可能會(huì)有所增加。上機(jī)作業(yè)2分別用Jacobi、Seidel、SOR(w=1.1,1.2,1.3,1.4,1.5)方法求解方程組Ax=b，這里A=-1211-12 1 1 1 1-1210×10，b=-1-110×1x0=0 . 0；e=10-5解：Jacobi迭代：公式：x(k+1)=(I-D-1

7、A)x(k)+ D-1b程序：A=-12 1 1 1 1 1 1 1 1 1;1 -12 1 1 1 1 1 1 1 1;1 1 -12 1 1 1 1 1 1 1;1 1 1 -12 1 1 1 1 1 1;1 1 1 1 -12 1 1 1 1 1;1 1 1 1 1 -12 1 1 1 1;1 1 1 1 1 1 -12 1 1 1;1 1 1 1 1 1 1 -12 1 1; 1 1 1 1 1 1 1 1 -12 1; 1 1 1 1 1 1 1 1 1 -12;b=-1 -1 -1 -1 -1 -1 -1 -1 -1 -1'N=500;eps=1e-5;n=length(b

8、);x0=zeros(n,1);x=zeros(n,1); k=0;D=zeros(n,n);for i=1:nfor j=1:nif i=jD(i,j)=A(i,j);endendendBJ=eye(n)-inv(D)*A;dJ=inv(D)*bwhile k<Nx=BJ*x0+dJ;if norm(x-x0,inf)<eps,break;endx0=x;k=k+1;endif k=N,warning('已達(dá)到迭代次數(shù)上限');enddisp('k= ',num2str(k),x運(yùn)行結(jié)果如下：dJ = 0.0833 0.0833 0.0833 0.

9、0833 0.0833 0.0833 0.0833 0.0833 0.0833 0.0833k= 32x = 0.3333 0.3333 0.3333 0.3333 0.3333 0.3333 0.3333 0.3333 0.33330.3333故x=0.3333 0.3333 0.3333 0.3333 0.3333 0.3333 0.3333 0.3333 0.3333 0.3333Gauss-Seidel迭代:公式：x(k+1)=-(L+D)-1U x(k)+ (L+D)-1b程序：A=-12 1 1 1 1 1 1 1 1 1;1 -12 1 1 1 1 1 1 1 1;1 1 -12

10、 1 1 1 1 1 1 1;1 1 1 -12 1 1 1 1 1 1;1 1 1 1 -12 1 1 1 1 1;1 1 1 1 1 -12 1 1 1 1;1 1 1 1 1 1 -12 1 1 1;1 1 1 1 1 1 1 -12 1 1; 1 1 1 1 1 1 1 1 -12 1; 1 1 1 1 1 1 1 1 1 -12;b=-1 -1 -1 -1 -1 -1 -1 -1 -1 -1'N=500;ep=1e-5;n=length(b);x0=zeros(n,1);x=zeros(n,1); k=0;while k<N for i=1:n if i=1 x(1)=

11、(b(1)-A(1,2:n)*x0(2:n)/A(1,1); else if i=n x(n)=(b(n)-A(n,1:n-1)*x(1:n-1)/A(n,n); else x(i)=(b(i)-A(i,1:i-1)*x(1:i-1)-A(i,i+1:n)*x0(i+1:n)/A(i,i); end end end if norm(x-x0,inf)<ep, break; end x0=x;k=k+1; end if k=N,warning('已達(dá)到迭代次數(shù)上限');end disp('k= ',num2str(k),x運(yùn)行結(jié)果如下：k= 18x = 0.

12、3333 0.3333 0.3333 0.3333 0.3333 0.3333 0.3333 0.3333 0.33330.3333故x=0.3333 0.3333 0.3333 0.3333 0.3333 0.3333 0.3333 0.3333 0.3333 0.3333SOR迭代：公式：x(k+1)=(D+L)-1(1-)D-U x(k)+ (D+L)-1b程序：w=1.1A=-12 1 1 1 1 1 1 1 1 1;1 -12 1 1 1 1 1 1 1 1;1 1 -12 1 1 1 1 1 1 1;1 1 1 -12 1 1 1 1 1 1;1 1 1 1 -12 1 1 1 1

13、 1;1 1 1 1 1 -12 1 1 1 1;1 1 1 1 1 1 -12 1 1 1;1 1 1 1 1 1 1 -12 1 1; 1 1 1 1 1 1 1 1 -12 1; 1 1 1 1 1 1 1 1 1 -12;b=-1 -1 -1 -1 -1 -1 -1 -1 -1 -1'n=length(b);N=500;eps=1e-5;x0=zeros(n,1);omega=1.1;x=zeros(n,1); k=0;while k<Nfor i=1:nif i=1x1(1)=(b(1)-A(1,2:n)*x0(2:n)/A(1,1);else if i=nx1(n)=

14、(b(n)-A(n,1:n-1)*x(1:n-1)/A(n,n);elsex1(i)=(b(i)-A(i,1:i-1)*x(1:i-1)-A(i,i+1:n)*x0(i+1:n)/A(i,i);endendx(i)=(1-omega)*x0(i)+omega*x1(i);endif norm(x0-x,inf)<eps, break; endk=k+1;x0=x;endif k=N,warning('已達(dá)到迭代次數(shù)上限');enddisp('k= ',num2str(k),x運(yùn)行結(jié)果如下：b = -1 -1 -1 -1 -1 -1 -1 -1 -1 -1k

15、= 15x = 0.3333 0.3333 0.3333 0.3333 0.3333 0.3333 0.3333 0.3333 0.33330.3333故x=0.3333 0.3333 0.3333 0.3333 0.3333 0.3333 0.3333 0.3333 0.3333 0.3333w=1.2A=-12 1 1 1 1 1 1 1 1 1;1 -12 1 1 1 1 1 1 1 1;1 1 -12 1 1 1 1 1 1 1;1 1 1 -12 1 1 1 1 1 1;1 1 1 1 -12 1 1 1 1 1;1 1 1 1 1 -12 1 1 1 1;1 1 1 1 1 1 -

16、12 1 1 1;1 1 1 1 1 1 1 -12 1 1; 1 1 1 1 1 1 1 1 -12 1; 1 1 1 1 1 1 1 1 1 -12;b=-1 -1 -1 -1 -1 -1 -1 -1 -1 -1'n=length(b);N=500;eps=1e-5;x0=zeros(n,1);omega=1.2;x=zeros(n,1); k=0;while k<Nfor i=1:nif i=1x1(1)=(b(1)-A(1,2:n)*x0(2:n)/A(1,1);else if i=nx1(n)=(b(n)-A(n,1:n-1)*x(1:n-1)/A(n,n);elsex

17、1(i)=(b(i)-A(i,1:i-1)*x(1:i-1)-A(i,i+1:n)*x0(i+1:n)/A(i,i);endendx(i)=(1-omega)*x0(i)+omega*x1(i);endif norm(x0-x,inf)<eps, break; endk=k+1;x0=x;endif k=N,warning('已達(dá)到迭代次數(shù)上限');enddisp('k= ',num2str(k),x運(yùn)行結(jié)果如下：b = -1 -1 -1 -1 -1 -1 -1 -1 -1 -1k= 11x = 0.3333 0.3333 0.3333 0.3333 0.

18、3333 0.3333 0.3333 0.3333 0.33330.3333故x=0.3333 0.3333 0.3333 0.3333 0.3333 0.3333 0.3333 0.3333 0.3333 0.3333 w=1.3A=-12 1 1 1 1 1 1 1 1 1;1 -12 1 1 1 1 1 1 1 1;1 1 -12 1 1 1 1 1 1 1;1 1 1 -12 1 1 1 1 1 1;1 1 1 1 -12 1 1 1 1 1;1 1 1 1 1 -12 1 1 1 1;1 1 1 1 1 1 -12 1 1 1;1 1 1 1 1 1 1 -12 1 1; 1 1 1

19、 1 1 1 1 1 -12 1; 1 1 1 1 1 1 1 1 1 -12;b=-1 -1 -1 -1 -1 -1 -1 -1 -1 -1'n=length(b);N=500;eps=1e-5;x0=zeros(n,1);omega=1.3;x=zeros(n,1); k=0;while k<Nfor i=1:nif i=1x1(1)=(b(1)-A(1,2:n)*x0(2:n)/A(1,1);else if i=nx1(n)=(b(n)-A(n,1:n-1)*x(1:n-1)/A(n,n);elsex1(i)=(b(i)-A(i,1:i-1)*x(1:i-1)-A(i,i+

20、1:n)*x0(i+1:n)/A(i,i);endendx(i)=(1-omega)*x0(i)+omega*x1(i);endif norm(x0-x,inf)<eps, break; endk=k+1;x0=x;endif k=N,warning('已達(dá)到迭代次數(shù)上限');enddisp('k= ',num2str(k),x運(yùn)行結(jié)果如下：b = -1 -1 -1 -1 -1 -1 -1 -1 -1 -1k= 9x = 0.3333 0.3333 0.3333 0.3333 0.3333 0.3333 0.3333 0.3333 0.33330.3333

21、故x=0.3333 0.3333 0.3333 0.3333 0.3333 0.3333 0.3333 0.3333 0.3333 0.3333 w=1.4A=-12 1 1 1 1 1 1 1 1 1;1 -12 1 1 1 1 1 1 1 1;1 1 -12 1 1 1 1 1 1 1;1 1 1 -12 1 1 1 1 1 1;1 1 1 1 -12 1 1 1 1 1;1 1 1 1 1 -12 1 1 1 1;1 1 1 1 1 1 -12 1 1 1;1 1 1 1 1 1 1 -12 1 1; 1 1 1 1 1 1 1 1 -12 1; 1 1 1 1 1 1 1 1 1 -1

22、2;b=-1 -1 -1 -1 -1 -1 -1 -1 -1 -1'n=length(b);N=500;eps=1e-5;x0=zeros(n,1);omega=1.4;x=zeros(n,1); k=0;while k<Nfor i=1:nif i=1x1(1)=(b(1)-A(1,2:n)*x0(2:n)/A(1,1);else if i=nx1(n)=(b(n)-A(n,1:n-1)*x(1:n-1)/A(n,n);elsex1(i)=(b(i)-A(i,1:i-1)*x(1:i-1)-A(i,i+1:n)*x0(i+1:n)/A(i,i);endendx(i)=(1-om

23、ega)*x0(i)+omega*x1(i);endif norm(x0-x,inf)<eps, break; endk=k+1;x0=x;endif k=N,warning('已達(dá)到迭代次數(shù)上限');enddisp('k= ',num2str(k),x運(yùn)行結(jié)果如下：b = -1 -1 -1 -1 -1 -1 -1 -1 -1 -1k= 11x = 0.3333 0.3333 0.3333 0.3333 0.3333 0.3333 0.3333 0.3333 0.33330.3333故x=0.3333 0.3333 0.3333 0.3333 0.3333

24、 0.3333 0.3333 0.3333 0.3333 0.3333 w=1.5A=-12 1 1 1 1 1 1 1 1 1;1 -12 1 1 1 1 1 1 1 1;1 1 -12 1 1 1 1 1 1 1;1 1 1 -12 1 1 1 1 1 1;1 1 1 1 -12 1 1 1 1 1;1 1 1 1 1 -12 1 1 1 1;1 1 1 1 1 1 -12 1 1 1;1 1 1 1 1 1 1 -12 1 1; 1 1 1 1 1 1 1 1 -12 1; 1 1 1 1 1 1 1 1 1 -12;b=-1 -1 -1 -1 -1 -1 -1 -1 -1 -1'

25、;n=length(b);N=500;eps=1e-5;x0=zeros(n,1);omega=1.5;x=zeros(n,1); k=0;while k<Nfor i=1:nif i=1x1(1)=(b(1)-A(1,2:n)*x0(2:n)/A(1,1);else if i=nx1(n)=(b(n)-A(n,1:n-1)*x(1:n-1)/A(n,n);elsex1(i)=(b(i)-A(i,1:i-1)*x(1:i-1)-A(i,i+1:n)*x0(i+1:n)/A(i,i);endendx(i)=(1-omega)*x0(i)+omega*x1(i);endif norm(x0-

26、x,inf)<eps, break; endk=k+1;x0=x;endif k=N,warning('已達(dá)到迭代次數(shù)上限');enddisp('k= ',num2str(k),x運(yùn)行結(jié)果如下：b = -1 -1 -1 -1 -1 -1 -1 -1 -1 -1k= 15x = 0.3333 0.3333 0.3333 0.3333 0.3333 0.3333 0.3333 0.3333 0.33330.3333故x=0.3333 0.3333 0.3333 0.3333 0.3333 0.3333 0.3333 0.3333 0.3333 0.3333三種

27、迭代方法對(duì)比如下表：表1：三種迭代方法對(duì)比表迭代方法迭代次數(shù)近似解Jacobi迭代法32x=0.3333 0.3333 0.3333 0.3333 0.3333 0.3333 0.3333 0.3333 0.3333 0.3333Gauss-Seidel迭代法18SOR迭代法w=1.115w=1.211w=1.39w=1.411w=1.515結(jié)果分析：由以上計(jì)算過(guò)程可以看出，這三個(gè)迭代法都是收斂的，但不同的是Jacobi迭代法的收斂速度是最慢的，SOR迭代法收的斂速度是最快，且隨著松弛因子w的變化而不同，從表中可以看出，當(dāng)w=1.3的時(shí)候收斂速度最快。方程的近似解為x=0.3333 0.333

28、3 0.3333 0.3333 0.3333 0.3333 0.3333 0.3333 0.3333 0.3333。上機(jī)作業(yè)3用Newton法與最速下降法求方程組x2+2y2-2=0x2=y在(0.8 , 0.7) 附近的根。解：Newton法：迭代公式：x(k+1)= x(k)-F(x(k) -1 F(x(k)程序：x0=0.8;0.7;for k=1:10z=-2*x0(1),4*x0(2);2*x0(1),-1x0(1)2+2*x0(2)2-2;x0(1)2-x0(2);x=x0+z;if norm(z)<1.0e-6;disp(x);disp(k)return;endx0=x;e

29、nd運(yùn)行結(jié)果： 0.8836 0.7808 4最速下降法：原理：對(duì)于方程組f1x1,x2=0f2x1,x2=0構(gòu)造函數(shù)x=f12x+f22x，xR2，求x的零最小值；即求x*R2，使得x*=0。x*也是方程組的解。M文件：function r,k= FastDown(F,x0,h,eps)format long;if nargin=3eps=1.0e-8;endn=length(x0);x0=transpose(x0);k=1;tol=1;while tol>epsfx=subs(F,findsym(F),x0); J=zeros(n,n); for i=1:n x1=x0; x1(i)

30、=x1(i)+h; J(:,i)=(subs(F,findsym(F),x1)-fx)/h; endlamda=fx/sum(diag(transpose(J)*J); r=x0-J*lamda;fr=subs(F,findsym(F),r);tol=dot(fr,fr); x0=r; k=k+1; if(k>500)disp('迭代步數(shù)太多,可能不收斂'); break; endendformat short;命令窗口：syms x y;f=x2+2*y2-2;x2-y;r,k=FastDown(f,0.8,0.7,1.0e-6)運(yùn)行結(jié)果如下：r = 0.8836 0.

31、7807k =19結(jié)果分析：Newton法具有較快的收斂速度，缺點(diǎn)是對(duì)初值要求過(guò)高，計(jì)算量較大；最速下降法對(duì)任意初值總是收斂的，缺點(diǎn)是收斂速度較慢，通常是越靠近解，逼近速度越慢。上機(jī)作業(yè)41.用冪法與反冪法求矩陣A的按模最大、最小特征值與對(duì)應(yīng)的特征向量。A=4 11 311-11 1-1 1 1 2 0 0 2解：冪法：原理：（1） 任取一非零向量u0=V0Rn，一般可取V0=1,1,1T（2） 計(jì)算Vk=Auk-1（3） mk=maxVk,uk=Vk/mk當(dāng)k足夠大時(shí)，即可得到：1mkukx1max(x1)程序：M文件：functionl,v,s=pmethod(A,x0,eps)if na

32、rgin=2eps=1.0e-5;endv=x0;N=500;m=0;l=0;for k=1:N y=A*v; m=max(y); v=y/m; if(abs(m-l)<eps) l=m; s=k; return; else if(k=N)disp('迭代步數(shù)太多,收斂速度太慢'); l=m; s=N; else l=m; end endend命令窗口：A=4 1 1 1;1 3 -1 1;1 -1 2 0;1 1 0 2x0=1 1 1 1'l,v,s=pmethod(A,x0)運(yùn)行結(jié)果如下：l = 5.2361v = 1.0000 0.6180 0.1180

33、0.5000s = 25故最大特征值為5.2361，對(duì)應(yīng)的特征向量為(1.0000,0.6180,0.1180,0.5000)T反冪法：原理：（1）任取一非零向量u0=V0Rn，一般可取V0=1,1,1T（2）計(jì)算Vk=A-1uk-1求解AVk=uk-1（3）mk=maxVk,uk=Vk/mk當(dāng)k足夠大時(shí)，即可得到：n1/mkukxnmax(xn)程序：M文件：functionl,v=ipmethod(A,x0,eps)if nargin=2eps=1.0e-6;endv=x0;N=500;m=0;l=0;for k=1:N y=Av; m=max(y); v=y/m; if(abs(m-l)

34、<eps) l=1/m; return; else if(k=N)disp('迭代次數(shù)太多,收斂速度太慢!'); l=1/m; else l=m; end endend命令窗口：A=4 1 1 1;1 3 -1 1;1 -1 2 0;1 1 0 2x0=1 1 1 1'l,v=ipmethod(A,x0)運(yùn)行結(jié)果如下：l = 0.7639v = -0.4721 0.7639 1.0000 -0.2361結(jié)果分析：冪法適用求矩陣的按模最大特征值及相應(yīng)的特征向量，算法簡(jiǎn)單，容易編寫(xiě)程序在計(jì)算機(jī)上實(shí)現(xiàn)，但收斂速度慢，其有效性依賴(lài)與矩陣特征值的分布情況。反冪法的適用范圍是

35、求矩陣的按模最小特征值及對(duì)應(yīng)的特征向量。這兩種算法的取舍決定于最大、最小特征值及相應(yīng)特征向量求得的難易程度。2.用Householder變換求矩陣A的QR分解，并用QR 方法做3次迭代A=4 11 311-11 1-1 1 1 2 0 0 2解：原理：ARm×n非奇異構(gòu)造Householder陣HkRm×n(k=1,2,n-1)使Hn-1Hn-2H2H1A=R(上三角陣)A=H1-1H2-1Hn-1-1R=H1H2Hn-1R=QR其中Q=H1H2Hn-1Rm×n為正交陣R=QTA=Hn-1Hn-2H2H1AHk-1=HkT=Hk,Q-1=QT=Hn-1Hn-2H2

36、H1程序：M文件1：function H,rho=householder(x,y)x=x(:);y=y(:);if length(x)=length(y) error(' X 和Y 必須要有相同的維數(shù)!');endrho=-sign(x(1)*norm(x)/norm(y);y=rho*y;v=(x-y)/norm(x-y);I=eye(length(x);H=I-2*v*v'M文件2：function Q,R=qrhs(A)n=size(A,1);R=A;Q=eye(n);for i=1:n-1 x=R(i:n,i); y=1;zeros(n-i,1);Ht=hous

37、eholder(x,y); H=blkdiag(eye(i-1),Ht); Q=Q*H;R=H*R; disp(Q);disp(R);end命令窗口：運(yùn)行結(jié)果如下：第一次QR分解：>>A=4 1 1 1;1 3 -1 1;1 -1 2 0;1 1 0 2;>>qrhs(A) -0.9177 -0.2294 -0.2294 -0.2294 -0.2294 0.9726 -0.0274 -0.0274 -0.2294 -0.0274 0.9726 -0.0274 -0.2294 -0.0274 -0.0274 0.9726 -4.3589 -1.6059 -1.1471 -

38、1.6059 0.0000 2.6882 -1.2569 0.6882 0.0000 -1.3118 1.7431 -0.3118 0.0000 0.6882 -0.2569 1.6882 -0.9177 0.1543 -0.3168 -0.1835 -0.2294 -0.8574 0.3895 -0.2462 -0.2294 0.4458 0.8647 0.0291 -0.2294 -0.2058 0.0132 0.9512 -4.3589 -1.6059 -1.1471 -1.6059 -0.0000 -3.0694 1.9034 -1.1146 0.0000 -0.0000 1.0231

39、 0.0990 0.0000 0.0000 0.1209 1.4727 -0.9177 0.1543 0.3362 -0.1451 -0.2294 -0.8574 -0.3579 -0.2902 -0.2294 0.4458 -0.8622 -0.0725 -0.2294 -0.2058 -0.1247 0.9431(R1)-4.3589 -1.6059 -1.1471 -1.6059 -0.0000 -3.0694 1.9034 -1.1146 -0.0000 0.0000 -1.0303 -0.2711 0.0000 0.0000 0.0000 1.4510(Q1)ans = -0.917

40、7 0.1543 0.3362 -0.1451 -0.2294 -0.8574 -0.3579 -0.2902 -0.2294 0.4458 -0.8622 -0.0725 -0.2294 -0.2058 -0.1247 0.9431第二次QR分解：>>R=-4.3589 -1.6059 -1.1471 -1.6059;0 -3.0694 1.9034 -1.1146;0 0 -1.0303 -0.2711;0 0 0 1.4510;>>Q=-0.9177 0.1543 0.3362 -0.1451; -0.2294 -0.8574 -0.3579 -0.2902; -

41、0.2294 0.4458 -0.8622 -0.0725; -0.2294 -0.2058 -0.1247 0.9431;>>A=R*Q;>>qrhs(A) -0.9907 -0.1037 -0.0591 0.0659 -0.1037 0.9946 -0.0031 0.0034 -0.0591 -0.0031 0.9982 0.0020 0.0659 0.0034 0.0020 0.9978 -5.0472 -0.8989 -0.3204 0.4616 0.0000 3.6356 -0.4358 -0.2571 0.0000 -0.4458 0.9037 -0.157

42、4 0 -0.2515 -0.1604 1.3421 -0.9907 0.1000 -0.0716 0.0589 -0.1037 -0.9850 0.1177 0.0716 -0.0591 0.1244 0.9905 -0.0024 0.0659 0.0652 -0.0018 0.9957 -5.0472 -0.8989 -0.3204 0.4616 -0.0000 -3.6714 0.5303 0.3274 0.0000 -0.0000 0.8448 -0.1930 0.0000 0 -0.1937 1.3220 -0.9907 0.1000 0.0829 0.0415 -0.1037 -0

43、.9850 -0.0987 0.0961 -0.0591 0.1244 -0.9660 0.2190 0.0659 0.0652 0.2243 0.9701(R2)-5.0472 -0.8989 -0.3204 0.4616 -0.0000 -3.6714 0.5303 0.3274 -0.0000 0.0000 -0.8667 0.4836 0.0000 -0.0000 0 1.2454(Q2)ans = -0.9907 0.1000 0.0829 0.0415 -0.1037 -0.9850 -0.0987 0.0961 -0.0591 0.1244 -0.9660 0.2190 0.06

44、59 0.0652 0.2243 0.9701第三次QR分解：>>R=-5.0472 -0.8989 -0.3204 0.4616;0 -3.6714 0.5303 0.3274;0 0 -0.8667 0.4836;0 0 0 1.2454;>>Q=-0.9907 0.1000 0.0829 0.0415;-0.1037 -0.9850 -0.0987 0.0961;-0.0591 0.1244 -0.9660 0.2190;0.0659 0.0652 0.2243 0.9701;>>A=R*Q;>>qrhs(A) -0.9972 -0.071

45、9 -0.0161 -0.0159 -0.0719 0.9974 -0.0006 -0.0006 -0.0161 -0.0006 0.9999 -0.0001 -0.0159 -0.0006 -0.0001 0.9999 -5.1575 -0.6363 -0.0973 -0.1111 -0.0000 3.6674 -0.0830 0.0740 -0.0000 -0.0844 0.9442 0.2778 0 0.0732 0.2779 1.2066 -0.9972 0.0718 -0.0178 -0.0145 -0.0719 -0.9969 0.0224 -0.0205 -0.0161 0.02

46、36 0.9996 0.0001 -0.0159 -0.0194 0.0001 0.9997 -5.1575 -0.6363 -0.0973 -0.1111 0.0000 -3.6691 0.0991 -0.0916 -0.0000 0.0000 0.9422 0.2797 0.0000 0 0.2797 1.2050 -0.9972 0.0718 0.0212 -0.0088 -0.0719 -0.9969 -0.0156 -0.0260 -0.0161 0.0236 -0.9583 -0.2844 -0.0159 -0.0194 -0.2846 0.9583(R3)-5.1575 -0.6

47、363 -0.0973 -0.1111 0.0000 -3.6691 0.0991 -0.0916 0.0000 -0.0000 -0.9828 -0.6111 0.0000 -0.0000 0 1.0755(Q3)ans = -0.9972 0.0718 0.0212 -0.0088 -0.0719 -0.9969 -0.0156 -0.0260 -0.0161 0.0236 -0.9583 -0.2844 -0.0159 -0.0194 -0.2846 0.9583結(jié)果分析：QR方法是冪法的一種推廣和變形，可以用來(lái)求任意矩陣的全部特征值。且對(duì)于可逆矩陣A，并上三角矩陣R的對(duì)角元都為正值，則

48、QR分解唯一。上機(jī)作業(yè)5目的：觀察lagrange插值的Runge現(xiàn)象。對(duì)于函數(shù)fx=11+25x2進(jìn)行l(wèi)agrange插值。取不同的等分?jǐn)?shù)n(n=5，n=10)，將區(qū)間-1，1n等分，取等距節(jié)點(diǎn)。把插值多項(xiàng)式的曲線畫(huà)在同一張圖上進(jìn)行比較。解：公式：Lnx=k=0nykln,kx程序：n=5 x-1-0.6-0.20.20.61f(x)0.0380.10.50.50.10.038n=10x-1-0.8-0.6-0.4-0.200.20.40.60.81f(x)0.0380.0590.10.1110.510.50.1110.10.0590.038M文件：function y=lagrange(x

49、0,y0,x)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(j)/(x0(k)-x0(j); end end s=p*y0(k)+s; end y(i)=s;end運(yùn)行結(jié)果如下：n=5：>>syms t>> x=-1:0.4:1;>> y=0.038,0.1,0.5,0.5,0.1,0.038;>> L1=lagrange(x,y,t)L1 =95/1536*(1/2*t+1/2)*(t+3/5)*(t+1/5)*(t-1/5)*(t-3/5)-125/192*(5/8*t+5/8)*(t+3/5)*(t+1/5)*(t-1/5)*(t