1.给定n 阶方程组Ax b =，其中6186186186A ⎛⎫ ⎪ ⎪ ⎪= ⎪ ⎪ ⎪⎝⎭ ，7151514b ⎛⎫ ⎪ ⎪ ⎪= ⎪⎪⎪⎝⎭则方程组有解(1,1,,1)Tx = 。

对10n =和84n =，分别用Gauss 消去法和列主元消去法解方程组，并比较计算结果。

Gauss 消去法：Matlab 编程（建立GS.m 文件）：function x=GS(n) A=[];b=[]; for i=1:n-1 A(i,i)=6; A(i,i+1)=1; A(i+1,i)=8; b(i)=15; endA(n,n)=6;b(1)=7;b(n)=14;b=b'; for k=1:n-1 for i=k+1:nm(i,k)=A(i,k)/A(k,k);A(i,k:n)=A(i,k:n)-m(i,k)*A(k,k:n); b(i)=b(i)-m(i,k)*b(k); end endb(n)=b(n)/A(n,n); for i=n-1:-1:1b(i)=(b(i)-sum(A(i,i+1:n).*b(i+1:n)'))/A(i,i); end clear x; x=b;disp( 'AX=b 的解x 是') end计算结果：在matlab 命令框里输出GS （10）得： >> GS(10)AX=b 的解x 是 ans =1.0000 1.00001.00001.00001.00001.00001.00001.00001.0000在matlab命令框里输出GS（84）得：>> GS(84)AX=b的解x是ans =1.0e+008 *0.0000………0.0000-0.00000.0000-0.00000.0000-0.00000.0000-0.00000.0000-0.00000.0000-0.00000.0000-0.00010.0002-0.00030.0007-0.00130.0026-0.00520.0105-0.02090.0419-0.08360.1665-0.3303-1.25822.3487-4.02635.3684列主元消去法：Matlab编程（建立GLZX.m文件）：function x=GLZX(n)A=[];b=[];eps=10^-2;for i=1:n-1A(i,i)=6;A(i,i+1)=1;A(i+1,i)=8;b(i)=15;endA(n,n)=6;b(1)=7;b(n)=14;b=b';for k=1:n-1[mainElement,index]=max(abs(A(k:n,k)));index=index+k-1;%indexif abs(mainElement)<epsdisp('列元素太小！！');break;elseif index>ktemp=A(k,:);temp1=b(k);A(k,:)=A(index,:);b(k)=b(index);A(index,:)=temp;b(index)=temp1;endfor i=k+1:nm(i,k)=A(i,k)/A(k,k);%A(k,k) ;A(i,k:n)=A(i,k:n)-m(i,k)*A(k,k:n);b(i)=b(i)-m(i,k)*b(k);endendb(n)=b(n)/A(n,n);for i=n-1:-1:1b(i)=(b(i)-sum(A(i,i+1:n).*b(i+1:n)'))/A(i,i); endclear x;x=b;disp('AX=b的解x是')end在matlab命令框里输出GLZX（10）得：>> GLZX(10)AX=b的解x是ans =1111111111在matlab命令框里输出GLZX（84）得：>> GLZX(84)AX=b的解x是ans =1.00001.00001.00001.00001.00001.0000………1.00001.00001.00001.00001.00001.00001.00001.0000分析：n较小时，两种方法均能得到正确解，当n较大后，Gauss消去法计算结果严重偏离准确值成为错解，列主元消去法依然能得到正确解。

这是因为Gauss消去法中有小主元做除数，在计算过程中的舍入误差会对解产生极大影响，从而导致错误。

列主元消去法则避免了这种情况。

2. 设有方程组Ax b =，其中2020A R ⨯∈，30.50.250.530.50.250.50.50.250.530.50.250.53A --⎛⎫⎪-- ⎪ ⎪-- ⎪= ⎪ ⎪-- ⎪-- ⎪ ⎪--⎝⎭(a) 选取不同的初始向量(0)x和不同的右端向量b ，给定迭代误差要求，用Jacobi 和Gauss-Seidel 迭代法计算，观测得出的迭代向量序列是否收敛。

若收敛，记录迭代次数，分析计算结果并得出你的结论。

(b) 选定初始向量初始向量(0)x和不同的右端向量b ，如取(0)0,,(1,1,1)T xb Ae e === 。

将A 的主对角线元素成倍增长若干次，非主对角元素不变，每次用Jacobi 法计算，要求迭代误差满足(1)()610k k xx +-∞-<，比较收敛速度，分析现象并得出你的结论。

(a)Matlab 编程： Jacobi 迭代法：function x=Jacobi(b,x0) A=zeros(20); for i=1:18A(i,i)=3;A(i+1,i)=-0.5;A(i,i+1)=-0.5;A(i+2,i)=-0.25;A(i,i+2)=-0.25; endA(19,19)=3;A(20,20)=3;A(19,20)=-0.5;A(20,19)=-0.5; D=diag(diag(A));L=tril(A-D);U=triu(A-D); B=D\(L+U);f=D\b; x=x0;i=0;ep=10^-5; while i<1000 y=x;x=B*x+f;i=i+1; if norm(x-y,inf)<ep break; end endif i<1000disp('迭代收敛且迭代次数为') i elsedisp('迭代不收敛')endendGauss-Seidel迭代法：function x=G_S(b,x0)A=zeros(20);for i=1:18A(i,i)=3;A(i+1,i)=-0.5;A(i,i+1)=-0.5;A(i+2,i)=-0.25;A(i,i+2)=-0.25;endA(19,19)=3;A(20,20)=3;A(19,20)=-0.5;A(20,19)=-0.5;D=diag(diag(A));L=tril(A-D);U=triu(A-D);B=(D-L)\U;f=(D-L)\b;x=x0;i=0;ep=10^-5;while i<1000y=x;x=B*x+f;i=i+1;if norm(x-y,inf)<epbreak;endendif i<1000disp('迭代收敛且迭代次数为')ielsedisp('迭代不收敛')endend计算结果：x和不同的右端向量b：设误差要求为1×10^-5任选了两组不同的初始向量(0)b=[1,1,1,1,1,1,1,1,2,3,4,1,4,2,3,1,3,1,4,6]' x0=[2,3,1,4,4,5,5,3,1,3,1,1,1,1,1,1,1,1,1,1]' b=[1,4,5,12,42,0,0,0,0,0,0,0,0,0,0,0,0,0,0,0]' x0=[0,1,2,1,1,1,1,1,1,1,1,1,1,1,1,1,1,1,1,1]' 在matlab内运行两个程序得到：第一组：>> Jacobi(b,x0)迭代收敛且迭代次数为i =19>> G_S(b,x0)迭代收敛且迭代次数为i =9第二组：>> Jacobi(b,x0)迭代收敛且迭代次数为i =18>> G_S(b,x0)迭代收敛且迭代次数为i =9分析：下Gauss-Seidel迭代法要比Jacobi迭代法收敛速度快。

因为Gauss-Seidel 迭代法对已经前面计算出来的信息进行了充分利用。

(b)Matlab编程：function x=Jacobi2(n)A=zeros(20);x=diag(A);e=diag(eye(20));for i=1:18A(i,i)=3*n;A(i+1,i)=-0.5;A(i,i+1)=-0.5;A(i+2,i)=-0.25;A(i,i+2)=-0.25;endA(19,19)=3*n;A(20,20)=3*n;A(19,20)=-0.5;A(20,19)=-0.5;b=A*e;D=diag(diag(A));L=tril(A-D);U=triu(A-D);B=D\(L+U);f=D\b;i=0;ep=10^-6;while i<1000y=x;x=B*x+f;i=i+1;if norm(x-y,inf)<epbreak;endendif i<1000disp('迭代收敛且迭代次数为')ielsedisp('迭代不收敛')endend计算结果：n代表主元增长倍数，选取不同的n得：>> Jacobi2(1)迭代收敛且迭代次数为i =20>> Jacobi2(2)迭代收敛且迭代次数为i =11>> Jacobi2(3)迭代收敛且迭代次数为i =9>> Jacobi2(4)迭代收敛且迭代次数为i =8>> Jacobi2(10)迭代收敛且迭代次数为i =6分析：矩阵主对角元增长倍数的变大，其收敛速度变快。

3.用迭代法求方程32⨯。

0.510-310+-=的全部根，要求误差限为8x x首先经过简单计算并结合图形得知该方程的三个单根区间为[-3,-1],[-1,0],[0,1],然后利用二分法求解。

Matlab编程：function x=f(a,b)fa=a^3+3*a^2-1;fb=b^3+3*b^2-1;if fa*fb>0disp('区间不是单根区间')return;endep=0.5*10^-8;while abs(a-b)>epc=(a+b)/2;f1=a^3+3*a^2-1; f2=c^3+3*c^2-1; j=f1*f2; if j==0 x=c; break; elseif j<0 b=c; else a=c; end x=a; end enddisp('该区间内根为') end计算结果： >> f(-3,-1) 该区间内根为 ans =-2.879385244101286>> f(-1,0)该区间内根为 ans =-0.652703646570444>> f(0,1)该区间内根为 ans =0.5320888832211494. 给定数据如下表：编制程序求三次样条插值函数在插值中点的样条函数插值，并作点集{,}i i x y 和样条插值函数的图形，满足的边界条件为（1）(0)0.8,(10)0.2.s s ''== （2）(0)(10)0.s s ''''==(1)Matlab编程：clear allx=[0 1 2 3 4 5 6 7 8 9 10];f1=0.8;f2=0.2;y=[0 0.79 1.53 2.19 2.71 3.03 3.27 2.89 3.06 3.19 3.29];n=length(x);h=[];t=[];g=[];u=[];s=[];A=2*eye(n-2);yn=[];xn=[];for k=2:nhk=x(n)-x(n-1);h=[h hk];endfor k=2:length(h)tk=h(k)/(h(k)+h(k-1));uk=h(k-1)/(h(k)+h(k-1));gk=3*(uk*(y(k+1)-y(k))/h(k)+tk*(y(k)-y(k-1))/h(k-1));g=[g gk];t=[t,tk];u=[u uk];endfor k=2:n-2A(k-1,k)=u(k-1);A(k,k-1)=t(k);endg(1)=g(1)-t(1)*f1;g(n-2)=g(n-2)-u(n-2)*f2;m=A\g';m=m';m=[f1 m f2];for k=1:n-1a1=2*y(k)/h(k)^3;a2=-2*y(k+1)/h(k)^3;a3=m(k)/h(k)^2; a4=m(k+1)/h(k)^2;b1=y(k)/h(k)^2; b2=y(k+1)/h(k)^2;c1=[1 -x(k)];c2=[1 -x(k+1)];c3=conv(c1,c1);c4=conv(c2,c2);sk1=b1*c4+b2*c3;sk2=(a1+a3)*conv(c1,c4)+(a2+a4)*conv(c2,c3);sk=[0 sk1]+sk2;s=[s;sk];Pn=poly2str(sk,'x');Y=polyval(s(k,:),(x(k)+x(k+1))/2)endfor k=1:n-1xi=[x(k):0.02:x(k+1)-0.02];yi=polyval(s(k,:),xi);yn=[yn yi];xn=[xn xi];endxn=[xn x(n)];yn=[yn polyval(s(n-1,:),x(n))];plot(xn,yn, 'x');title('3次样条函数插值结果')计算结果：在Matlab命令框内输入上述程序得到插值中点的样条函数值Y和样条插值函数的图形：Y =0.398564254606255Y =1.168428726968726Y =1.871470837518841Y =2.478187922956004Y =2.873277470657021Y =3.213702194414573Y =3.084413751685133Y =2.919892798844714Y =3.149765052935493Y =3.222296989411632(2)Matlab编程：clear allx=[0 1 2 3 4 5 6 7 8 9 10];f1=0;f2=0;y=[0 0.79 1.53 2.19 2.71 3.03 3.27 2.89 3.06 3.19 3.29];n=length(x);h=[];t=[];g=[];u=[];s=[];A=2*eye(n);yn=[];xn=[];for k=2:nhk=x(n)-x(n-1);h=[h hk];endfor k=2:length(h)tk=h(k)/(h(k)+h(k-1));uk=h(k-1)/(h(k)+h(k-1));gk=3*(uk*(y(k+1)-y(k))/h(k)+tk*(y(k)-y(k-1))/h(k-1));g=[g gk];t=[t,tk];u=[u uk];endfor k=1:n-2A(k+1,k+2)=u(k);A(k+1,k)=t(k);endA(1,2)=1;A(n,n-1)=1;g0=3*(y(2)-y(1))/h(1)-h(1)*f1/2;gn=3*(y(n)-y(n-1))/h(n-2)+h(n-2)*f2/2;g=[g0,g,gn];m=A\g';m=m';for k=1:n-1a1=2*y(k)/h(k)^3;a2=-2*y(k+1)/h(k)^3;a3=m(k)/h(k)^2; a4=m(k+1)/h(k)^2;b1=y(k)/h(k)^2; b2=y(k+1)/h(k)^2;c1=[1 -x(k)];c2=[1 -x(k+1)];c3=conv(c1,c1);c4=conv(c2,c2);sk1=b1*c4+b2*c3;sk2=(a1+a3)*conv(c1,c4)+(a2+a4)*conv(c2,c3);sk=[0 sk1]+sk2;s=[s;sk];Pn=poly2str(sk,'x');Y=polyval(s(k,:),(x(k)+x(k+1))/2)endfor k=1:n-1xi=[x(k):0.02:x(k+1)-0.02];yi=polyval(s(k,:),xi);yn=[yn yi];xn=[xn xi];endxn=[xn x(n)];yn=[yn polyval(s(n-1,:),x(n))];plot(xn,yn, 'x');title('3次样条函数插值结果')计算结果：在命令框内输入上述程序得到：Y =0.398428148708663Y =1.168465553874013Y =1.871459635795274Y =2.478195902944736Y =2.873256752425371Y =3.213777087353492Y =3.084134898160016Y =2.920933320005332Y =3.145881821815571Y =3.236789392730741x y的图形再次输入plot(x,y, 'x');可以得到点集{,}i i5.对下列数据作三次多项式拟合，取权系数1w ，给出拟合多项式的系数、平方误差，ix y和拟合多项式的图形。