计算方法上机作业集合

第一次&第二次上机作业上机作业：1.在Matlab上执行：5.1-5-0.1和1.5-1-0.5给出执行结果，并简要分析一下产生现象的原因。解：执行结果如下：在Matlab中，小数值很难用二进制进行描述。由于计算精度的影响，相近两数相减会出现误差。2.（课本181页第一题）解：（1）n=0时，积分得𝐼0=ln6-ln5，编写如下图代码从以上代码显示的结果可以看出，𝐼20的近似值为0.012712966517465（2）𝐼𝑛=∫𝑥𝑛5+𝑥10dx,可得∫𝑥𝑛610dx≤∫𝑥𝑛5+𝑥10dx≤∫𝑥𝑛510dx,得16(𝑛+1)≤𝐼𝑛≤15(𝑛+1),则有1126≤𝐼20≤1105,取𝐼20=1105,以此逆序估算𝐼0。代码段及结果如下图所示结果是从𝐼19逆序输出至𝐼0，所以得到𝐼0的近似值为0.088392216030227。（3）从𝐼20估计的过程更为可靠。首先根据积分得表达式是可知，被积函数随着n的增大，其所围面积应当是逐步减小的，即积分值应是随着n的递增二单调减小的，（1）中输出的值不满足这一条件，（2）满足。设𝑆𝑛表示𝐼𝑛的近似值，𝑆𝑛-𝐼𝑛=(−5)𝑛(𝑆0−𝐼0)(根据递推公式可以导出此式),可以看出，随着n的增大，误差也在增大，所以顺序估计时，算法不稳定性逐渐增大，逆序估计情况则刚好相反，误差不断减小，算法逐渐趋于稳定。2.（课本181页第二题）（1）上机代码如图所示求得近似根为0.09058（2）上机代码如图所示得近似根为0.09064；（3）牛顿法上机代码如下计算所得近似解为0.09091第三次上机作业上机作业181页第四题线性方程组为[1.13483.83260.53011.78751.16513.40172.53301.54353.41294.93171.23714.99988.76431.314210.67210.0147][𝑥1𝑥2𝑥3𝑥4]=[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;ifdif0disp('此方程组无解');returnendifRA==RBifRA==nformatlong;disp('此方程组有唯一解');forp=1:n-1fork=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);forq=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;formatlong;ifdif0disp('该方程组无解');returnendifdif==0ifRA==ndisp('该方程组有唯一解');c=zeros(1,n);fori=1:n-1max=abs(A(i,i));m=i;forj=i+1:nifmaxabs(A(j,i))max=abs(A(j,i));m=j;endend%求出每一次消元时绝对值最大的一行的行号ifm~=ifork=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;%函数值跟随方程一起换位置endfork=i+1:nforj=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);fori=n-1:-1:1sum=0;forj=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.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=lieZhu(A,b)该方程组有唯一解X=1.0000000000000001.0000000000000020.9999999999999990.999999999999999根据两种方法运算结果，顺序消元法得到结果比列主元消元法得到的结果精度更高。（注：matlab使用的是2015b版本，不知道是保留小数位数太少，还是程序原因，顺序消元输出结果总是等于准确解，请老师指正）第四次上机作业7.分析用下列迭代法解线性方程组[4−1−140−1−1000−100−1−104−1−140−1−100−1000−1−104−1−14][𝑥1𝑥2𝑥3𝑥4𝑥5𝑥6]=[05−25−26]的收敛性，并求出使‖𝑋(𝑘+1)−𝑋(𝑘)‖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);endA=[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）SOR迭代法（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.000000413180530.9999918585434762.0000068413569迭代次数为13SOR(A,b,1.95)结果为0.9999849520881072.000009608326040.9999591151827292.00001684260061.000004435266971.99997885113446迭代次数为241SOR(A,b,0.95)结果为0.9999615183093511.999957068252310.9999630548384531.999965805720330.9999711417275891.9999827300678迭代次数为17由以上输出得到w取值不同的情况下，得到的满足精度要求的结果，迭代次数分别如输出所示第五次上机作业8.从函数表x0.00.10.20.30.4010.5f(x)0.398940.396950.391420.381380.368120.35206出发，用下列方法计算f(0.15),f(0.31)及f(0.47)的近似值（1）分段线性插值（2）分段二次插值（3）全区间上拉格朗日插值解：（1）线性插值编写函数如下functi