第一次&第二次上机作业上机作业：1.在Matlab上执行：>> 5.1-5-0.1和>> 1.5-1-0.5给出执行结果，并简要分析一下产生现象的原因。

解：执行结果如下：在Matlab中，小数值很难用二进制进行描述。

由于计算精度的影响，相近两数相减会出现误差。

2.（课本181页第一题）解：（1）n=0时，积分得I0=ln6-ln5，编写如下图代码从以上代码显示的结果可以看出，I 20的近似值为0.012712966517465（2）I n =∫x n 5+x 10dx,可得∫x n 610dx ≤∫x n 5+x 10dx ≤∫x n 510dx,得 16(n+1)≤I n ≤15(n+1),则有1126≤I 20≤1105, 取I 20=1105,以此逆序估算I 0。

代码段及结果如下图所示结果是从I 19逆序输出至I 0，所以得到I 0的近似值为0.088392216030227。

（3）从I 20估计的过程更为可靠。

首先根据积分得表达式是可知，被积函数随着n 的增大，其所围面积应当是逐步减小的，即积分值应是随着n 的递增二单调减小的，（1）中输出的值不满足这一条件，（2）满足。

设S n 表示I n 的近似值，S n -I n =(−5)n (S 0−I 0) 根据递推公式可以导出此式),可以看出，随着n 的增大，误差也在增大，所以顺序估计时，算法不稳定性逐渐增大，逆序估计情况则刚好相反，误差不断减小，算法逐渐趋于稳定。

2.（课本181页第二题）（1）上机代码如图所示求得近似根为0.09058（2）上机代码如图所示得近似根为0.09064；（3）牛顿法上机代码如下计算所得近似解为0.09091第三次上机作业上机作业181页第四题线性方程组为[1.1348 3.83260.5301 1.78751.1651 3.40172.5330 1.54353.41294.93171.2371 4.99988.7643 1.314210.67210.0147][x1x2x3x4]=[9.53426.394118.423116.9237]（1）顺序消元法A=[1.1348,3.8326,1.1651,3.4017;0.5301,1.7875,2.5330,1.5435;3.4129,4. 9317,8.7643,1.3142;1.2371,4.9998,10.6721,0.0147];b=[9.5342;6.3941;18.4231;16.9237];上机代码（函数部分）如下function [b] = gaus( A,b )%用b返回方程组的解B=[A,b];n=length(b);RA=rank(A);RB=rank(B);dif=RB-RA;if dif>0disp('此方程组无解');returnendif RA==RBif RA==nformat long;disp('此方程组有唯一解');for p=1:n-1for k=p+1:nm=B(k,p)/B(p,p);B(k,p:n+1)=B(k,p:n+1)-m*B(p,p:n+1);endend %顺序消元形成上三角矩阵b=B(1:n,n+1);A=B(1:n,1:n);b(n)=b(n)/A(n,n);for q=n-1:-1:1b(q)=(b(q)-sum(A(q,q+1:n)*b(q+1:n)))/A(q,q);end %回代求解elsedisp('此方程组有无数组解');endend上机运行结果为>>A=[1.1348,3.8326,1.1651,3.4017;0.5301,1.7875,2.5330,1.5435;3.4129,4. 9317,8.7643,1.3142;1.2371,4.9998,10.6721,0.0147];b=[9.5342;6.3941;18.4231;16.9237];>> X=gaus A,b)此方程组有唯一解X =1.0000000000000001.0000000000000001.0000000000000001.000000000000000>>（2）列主元消元法A=[1.1348,3.8326,1.1651,3.4017;0.5301,1.7875,2.5330,1.5435;3.4129,4. 9317,8.7643,1.3142;1.2371,4.9998,10.6721,0.0147];b=[9.5342;6.3941;18.4231;16.9237];函数部分代码如下function [b] = lieZhu(A,b )%用b返回方程组的解B=[A,b];RA=rank(A);RB=rank(B);n=length(b);dif=RB-RA;format long;if dif>0disp('该方程组无解');returnendif dif==0if RA==ndisp('该方程组有唯一解');c=zeros(1,n);for i=1:n-1max=abs(A(i,i));m=i;for j=i+1:nif max<abs(A(j,i))max=abs(A(j,i));m=j;endend %求出每一次消元时绝对值最大的一行的行号if m~=ifor k=i:nc(k)=A(i,k);A(i,k)=A(m,k);A(m,k)=c(k);endd1=b(i);b(i)=b(m);b(m)=d1;%函数值跟随方程一起换位置endfor k=i+1:nfor j=i+1:nA(k,j)=A(k,j)-A(i,j)*A(k,i)/A(i,i); endb(k)=b(k)-b(i)*A(k,i)/A(i,i);A(k,i)=0;endend %完成消元操作，形成上三角矩阵b(n)=b(n)/A(n,n);for i=n-1:-1:1sum=0;for j=i+1:nsum=sum+A(i,j)*b(j);endb(i)=(b(i)-sum)/A(i,i) ;%回代求解其他未知数endendelsedisp('此方程组有无数组解');endend上机运行，结果为>>A=[1.1348,3.8326,1.1651,3.4017;0.5301,1.7875,2.53 30,1.5435;3.4129,4.9317,8.7643,1.3142;1.2371,4.99 98,10.6721,0.0147];b=[9.5342;6.3941;18.4231;16.9237];X=lieZhu(A,b)该方程组有唯一解X =1.0000000000000001.0000000000000020.9999999999999990.999999999999999>>根据两种方法运算结果，顺序消元法得到结果比列主元消元法得到的结果精度更高。

（注：matlab使用的是2015b版本，不知道是保留小数位数太少，还是程序原因，顺序消元输出结果总是等于准确解，请老师指正）第四次上机作业7.分析用下列迭代法解线性方程组[ 4−1−1 40−1−1 00 0−1 00−1−1 04−1−1 40 −1−1 00−10 00−1−1 04−1−1 4][x1x2x3x4x5x6]=[5−25−26]的收敛性，并求出使‖X(k+1)−X(k)‖2≤0.0001的近似解及相应的迭代次数。

（1）雅可比迭代法解：上机编写的函数如下function [] = Jacobi(X,b)%雅可比迭代法D=diag(diag(X));%得到对角线元素组成的矩阵B=inv(D)*(D-X);f=inv(D)*b;b(:,:)=0;x1=B*b+f;num=1;while(norm(x1-b)>0.0001)%判断当前的解是否达到精度要求b=x1;x1=B*b+f;num=num+1;end;fprintf('求得的解为：\n');disp(x1);fprintf('迭代次数：%d次\n',num);end上机运行，结果如下>> A=[4,-1,0,-1,0,0;-1,4,-1,0,-1,0;0,-1,4,-1,0,-1;-1,0,-1,4,-1,0;0,-1,0,-1,4,-1;0,0,-1,0,-1,4];>> b=[0;5;-2;5;-2;6];>> Jacobi(A,b)求得的解为：0.9999817650743811.999950181256740.9999750906283681.999950181256740.9999750906283681.99996353014876迭代次数：28次满足要求的解如输出结果所示，总共迭代了28次（2）高斯-赛德尔迭代法上机程序如下所示function [] =Gauss_Seidel( A,b )%高斯赛德尔迭代法D=diag(diag(A));L=D-tril(A);U=D-triu(A);B=inv(D-L)*U;f=inv(D-L)*b;b(:,:)=0;x0=B*b+f;num=1;while(norm(x0-b)>0.0001)num=num+1;b=x0;x0=B*b+f;end;fprintf('结果为\n');disp(x0);fprintf('迭代次数为：%d次\n',num);end>> A=[4,-1,0,-1,0,0;-1,4,-1,0,-1,0;0,-1,4,-1,0,-1;-1,0,-1,4,-1,0;0,-1,0,-1,4,-1;0,0,-1,0,-1,4];>> b=[0;5;-2;5;-2;6];>> Gauss_Seidel(A,b)结果为0.9999601438106581.999956761521390.9999635082998331.999966621628740.9999725271797151.99998400886989迭代次数为：15次满足精度要求的解如上述程序打印输出所示，迭代了15次（3）S OR迭代法（w依次取1.334,1.95,0.95）上机代码如下function [] = SOR(A,b,w )%SOR迭代法¨D=diag(diag(A));L=D-tril(A);U=D-triu(A);B=inv(D-w*L)*((1-w)*D+w*U);f=w*inv(D-w*L)*b;b(:,:)=0;x0=B*b+f;num=1;while(norm(x0-b)>0.0001)num=num+1;b=x0;x0=B*b+f;end;fprintf('结果为\n');disp(x0);fprintf('迭代次数为%d\n',num);end上机运行>> A=[4,-1,0,-1,0,0;-1,4,-1,0,-1,0;0,-1,4,-1,0,-1;-1,0,-1,4,-1,0;0,-1,0,-1,4,-1;0,0,-1,0,-1,4]; >> b=[0;5;-2;5;-2;6];>> SOR(A,b,1.334)结果为1.000018784810091.999986983228581.000018150130682.00000041318053 0.999991858543476 2.0000068413569 迭代次数为13>> SOR(A,b,1.95)结果为0.999984952088107 2.00000960832604 0.999959115182729 2.0000168426006 1.00000443526697 1.99997885113446 迭代次数为241>> SOR(A,b,0.95)结果为0.9999615183093511.999957068252310.9999630548384531.999965805720330.9999711417275891.9999827300678 迭代次数为17由以上输出得到w取值不同的情况下，得到的满足精度要求的结果，迭代次数分别如输出所示第五次上机作业8.从函数表出发，用下列方法计算f(0.15),f(0.31)及f(0.47)的近似值（1）分段线性插值（2）分段二次插值（3）全区间上拉格朗日插值解：（1）线性插值编写函数如下function [R] = div_line( x0,y0,x )%线性插值p=length(x0);n=length(y0);m=length(x);if(p~=n)%x的个数与y的个数不等error('数据输入有误，请重新输入');return;elsefprintf('线性插值\n');for t=1:mz=x(t);if(z<x0(1)||z>x0(p))fprintf('x[%d]不在所给自变量范围内，无法进行插值',t);continue;elsefor i=1:p-1if(z<x0(i+1))break;end;end;R(t)=y0(i)*(x(t)-x0(i+1))/(x0(i)-x0(i+1))+y0(i+1)*(x(t)-x0(i))/(x0(i+1)-x0(i));end;end;end;end上机运行如下>> x0=[0.0 0.1 0.195 0.3 0.401 0.5];>> y0=[0.39894 0.39695 0.39142 0.38138 0.36812 0.35206]; >> x=[0.15 0.31 0.47];>> div_line(x0,y0,x)线性插值ans =0.39404 0.38007 0.35693即结果为f(0.15)≈0.39404,f(0.31) ≈0.38007,f(0.47) ≈0.35693 （2）分段二次插值编写的函数如下function [R] = div2line(x0,y0,x)%分段二次插值p=length(x0);m=length(y0);n=length(x);if(p~=m)error('输入错误，请重新输入数据');end;for t=1:nif(x(t)<x0(1)||x(t)>x0(p))fprintf('x[%d]不在所给区间上',t);continue;end;tag=2;m=abs(x(t)-x0(1))+abs(x(t)-x0(2))+abs(x(t)-x0(3));for i=3:p-1sum=abs(x(t)-x0(i-1))+abs(x(t)-x0(i))+abs(x(t)-x0(i+1));if(sum<m)m=sum;tag=i;endend;fprintf('tag=%d\n',tag);R(t)=y0(tag-1)*(x(t)-x0(tag))*(x(t)-x0(tag+1))/((x0(tag-1)-x0(tag))*(x0(tag-1)-x0(tag+1)))+y0(tag)*(x(t)-x0(tag-1))*(x(t)-x0(tag+1))/((x0(tag)-x0(tag-1))*(x0(tag)-x0(tag+1)))+y0(tag+1)*(x(t)-x0(tag-1))*(x(t)-x0(tag))/((x0(tag+1)-x0(tag-1))*(x0(tag+1)-x0(tag)));End上机运行，执行结果为>> x0=[0.0 0.1 0.195 0.3 0.401 0.5];>> y0=[0.39894 0.39695 0.39142 0.38138 0.36812 0.35206]; >> x=[0.15 0.31 0.47];>>div2line(x0,y0,x)ans =0.39448 0.38022 0.35725即分段二次插值方法下，f(0.15)≈0.39448,f(0.31) ≈0.38022,f(0.47) ≈0.35725（3）上机编写的程序如下function [R] = lagrange(x0,y0,x)%全区间上拉格朗日插值p=length(y0);n=length(x0);m=length(x);%计算函数表和x的长度if p ~= n error('数据输入有误，请重新输入');%若函数表的x与y长度不一致则输入有误else fprintf('拉格朗日插值\n');for t=1:m%利用循环计算每个x的插值s=0.0;z=x(t);for k=1:np=1;for i=1:nif i~=kp=p*(z-x0(i))/(x0(k)-x0(i));endends=s+y0(k)*p;end%根据拉格朗日插值公式求解yR(t)=s;%输出插值结果endend上机运行结果为>> x0=[0.0 0.1 0.195 0.3 0.401 0.5];>> y0=[0.39894 0.39695 0.39142 0.38138 0.36812 0.35206]; >> x=[0.15 0.31 0.47];>> lagrange(x0,y0,x)拉格朗日插值ans =0.39447 0.38022 0.35722即分段二次插值方法下，f(0.15)≈0.39447,f(0.31) ≈0.38022,f(0.47) ≈0.357229.解：上机程序如下，为方便起见，将所有操作分在四个函数中进行入口函数function [] =spline X,Y,xx,y1_0,y1_18 )%输出自变量所对应的函数值M=getM X,Y,y1_0,y1_18);%先得到Ms=xx;k=length xx);for a=1:ks xx a))=getExp X,Y,M,xx a));%计算自变量所在小区间对应曲线的表达式，并根据表达式计算对应的函数值fprintf 's %d)=%f\n',xx a),s xx a))); %输出打印函数值endend获取Mfunction [M] = getM X,Y,y1_0,y1_1)%得到Mn=length X);a=0*X;b=a;c=a;h=a;f=a;b=b+2;h 2:n)=X 2:n)-X 1:n-1); % h 1)不用a 2:n-1)=h 2:n-1)./ h 2:n-1)+h 3:n));c 2:n-1)=1-a 2:n-1);a n)=1;c 1)=1;Yf 2:n)=Y 2:n)-Y 1:n-1);f 2:n-1)=6* Yf 3:n)./h 3:n)-Yf 2:n-1)./h 2:n-1))./ h 2:n-1)+h 3:n));f 1)=6* Yf 2)/h 2)-y1_0)/h 2);f n)=6* y1_1-Yf n)/h n))/h n);M=CalM n,a,b,c,f);%计算Mend计算Mfunction [f] = CalM n,a,b,c,f)% 追赶法求解Meps=0.1e-15; %防止参数过小，是的计算误差过大if abs b 1))<epsdisp '除数为0，停止计算');returnelsec 1)=c 1)/b 1);end%追赶法：根据递推算式计算βfor i=2:n-1b i)=b i)-a i)*c i-1);if abs b i))<epsdisp '除数为0，停止计算');returnelsec i)=c i)/b i);endendb n)=b n)-a n)*c n-1);%追赶法:根据递推算式计算f 1)=f 1)/b 1);for i=2:nf i)= f i)-a i)*f i-1))/b i);end%以下求解Ux=y, x的值存入ffor i=n-1:-1:1f i)=f i)-c i)*f i+1);endreturnend得到自变量所在区间的表达式，并求自变量对应的函数值function [y] = getExp X,Y,M,x)%%根据X、Y、M计算表达式，并根据表达式计算对应的函数值n=length X);h 2:n)=X 2:n)-X 1:n-1);%%判断x落在哪个小区间n1=1;n2=n;while n2~=n1+1n5=fix n1+n2)/2);if x>X n5)n1=n5;elsen2=n5;endend%%%%计算yy=M n1)* X n2)-x)^3/ 6*h n2))+ M n2)* x-X n1))^3/ 6*h n2));y=y+ Y n1)-M n1)*h n2)*h n2)/6)* X n2)-x)/h n2);y=y+ Y n2)-M n2)*h n2)*h n2)/6)* x-X n1))/h n2);%%%结束end上机试运行，过程如下>> X=[0.52 3.1 8.0 17.95 28.65 39.62 50.65 78 104.6 156.6 208.6 260.7 312.5 364.4 416.3 468 494 507 520];>> Y=[5.28794 9.4 13.84 20.2 24.9 28.44 31.1 35 36.5 36.6 34.6 31.0 26.34 20.9 14.8 7.8 3.7 1.5 0.2];>> xx=[2 4 6 12 16 30 60 110 180 280 400 515];>> y1_0=1.86548;>> y1_18=-0.046115;spline(X,Y,xx,y1_0,y1_18)s(2)=7.825123s(4)=10.481311s(6)=12.363477s(12)=16.575574s(16)=19.091594s(30)=25.386402s(60)=32.804283s(110)=36.647878s(180)=35.917148s(280)=29.368462s(400)=16.799784s(515)=0.542713根据上述程序运行结果，可得所有自变量对应的函数值，如上输出结果所示第六次上机作业10.已知一组实验数据试用最小二乘法求它的多项式拟合曲线，并求出最低点位置。