西安交通大学计算方法大作业.pdf

《西安交通大学计算方法大作业.pdf》由会员分享，可在线阅读，更多相关《西安交通大学计算方法大作业.pdf（27页珍藏版）》请在淘文阁 - 分享文档赚钱的网站上搜索。

1、计算方法上机报告计算方法上机报告姓名：学号：班级：目录目录题目一-4-1.1 题目内容-4-1.2 算法思想-4-1.3Matlab 源程序-5-1.4 计算结果及总结-5-题目二-7-2.1 题目内容-7-2.2 算法思想-7-2.3 Matlab 源程序-8-2.4 计算结果及总结-9-题目三-11-3.1 题目内容-11-3.2 算法思想-11-3.3Matlab 源程序-13-3.4 计算结果及总结-14-题目四-15-4.1 题目内容-15-4.2 算法思想-15-4.3Matlab 源程序-15-4.4 计算结果及总结-16-题目五-18-5.1 题目内容-18-5.2 算法思想-

2、18-5.3 Matlab 源程序-18-5.3.1 非压缩带状对角方程组-18-5.3.2 压缩带状对角方程组-20-5.4 实验结果及分析-22-5.4.1Matlab 运行结果-22-5.4.2 总结分析-24-5.5 本专业算例-24-学习感悟-27-题目一题目一1.11.1 题目内容题目内容计算以下和式：S 1 4211，要求：n8n 48n58n6n0168n1(1)若保留 11 个有效数字，给出计算结果，并评价计算的算法；(2)若要保留 30 个有效数字，则又将如何进行计算。1.21.2 算法思想算法思想在程序编写中需要把握以下几点：随着 n 值的增加，和式的项递减速度很快，因此

3、我们可以认为，在确定为某一精度的前提下，n 达到一定的值，加下一项将不会对最终的加和产生影响，首先我们应找到 n 值。根据精度要求估计所加的项数，可以使用后验误差估计，通项为：an1421114nn168n 18n 48n 58n 6168n 1为减小舍入误差，在计算s 时所采用的方法是逆序相加，其依据是：两个数量级相差较大的数字相加减时，较小数的有效数字会被丧失，从而导致最后的运算结果失真。为避免“大数吃小数”现象的发生，采用逆序相加。对于实现 30 位有效数字，则调用从工具箱中 digits(位数)或 vpa（变量，精度位数）即可实现。-4-1.3Matlab1.3Matlab 源程序源程

4、序clear;clc;m=input(输入需要求的有效数字位数 m=);s=0for n=0:200%寻找满足条件的最小 ns=(1/16n)*(4/(8*n+1)-2/(8*n+4)-1/(8*n+5)-1/(8*n+6);if s0.0001续循环x1=x1-(6*x15-45*x12+20)/(30*x14-90*x1);end-15-%初值为-1 时%误差大于 0.0001 时继%Newton 迭代法x2=0.5;%初值为 0.5 时while abs(6*x25-45*x22+20)0.0001x2=x2-(6*x25-45*x22+20)/(30*x24-90*x2);endx3=

5、1.5;%初值为 1.5 时while abs(6*x35-45*x32+20)0.0001x3=x3-(6*x35-45*x32+20)/(30*x34-90*x3);endfprintf(x1=%.6f,x2=%.6f,x3=%.6f,x1,x2,x3)绘图程序x=-2:0.1:2y=6*x.5-45*x.2+20plot(x,y)xlabel(X 轴);ylabel(Y 轴);grid on;axis on;4.44.4 计算结果及总结计算结果及总结x1=-0.654542x2=0.681174x3=1.870799所得曲线-16-%输出三个根%坐标轴表示对象标签%显示网格线%显示坐标轴

6、所给的方程形式较为简单，导数较为好求所以用牛顿法非常方便，先作图知道根的大致范围后再取初值可以有效的减小迭代次数。-17-题目五题目五5.15.1 题目内容题目内容线性方程组求解。（1）编写程序实现大规模方程组的高斯消去法程序，并对所附的方程组进行求解。所附方程组的类型为对角占优的带状方程组。（2）针对本专业中所碰到的实际问题，提炼一个使用方程组进行求解的例子，并对求解过程进行分析、求解。5.25.2 算法思想算法思想消去法的中心是“降维”，即求解 n 元方程组的问题转化为先解 n-1 元方程组，一旦这个 n-1 元方程组的解取得，则剩余的一个未知量自然可以求得。这样逐步减少未知量个数的方法，

7、是求解多元方程组的一个重要思想。Gauss 消去过程中，适当交换方程的顺序对保证消去过程能够顺利进行及计算解的精确度都是必要的。消去过程中产生的数称为第 k 步消去的主元。交换方程的原则是使中，绝对值最大的一个换到(k,k)位置而成为第 k 步消去的主元，带有这种交换的 Gauss 消去法为列主元 Gauss 消去法。5.3 Matlab5.3 Matlab 源程序源程序5.3.15.3.1 非压缩带状对角方程组非压缩带状对角方程组（1）dat51.dat、dat52.dat 主程序:%非压缩格式 dat51.dat;dat52.dat-18-clear;clc;format short%读出

8、数据fid=fopen(E:2015 上机题目实验五dat52.dat,r);D=fread(fid,3,long,0);n=D(3)%矩阵维数offset=20;fseek(fid,offset,bof);A=zeros(n);B=zeros(1,n);C=fread(fid,n*(n+1),float,0)for i=1:nfor j=1:nA(i,j)=C(i-1)*n+j);endendA;for i=1:nB(i)=C(n*n+i);endB;x=GAUSSPP(A,B)（2）列主元高斯消去法 Matlab 函数程序：%列主元高斯消去法function x=GAUSSPP(A,B)n

9、=length(B);%-消去过程-%找出列主元for k=1:n-1MAX=abs(A(k,k);Rmax=k;for i=k+1:nif abs(A(i,k)MAXMAX=abs(A(i,k);Rmax=i;endend-19-if A(Rmax,k)=0break;end%交换当前行和主元所在行for j=1:ntemp=A(k,j);A(k,j)=A(Rmax,j);A(Rmax,j)=temp;endtemp=B(k);B(k)=B(Rmax);B(Rmax)=temp;for i=k+1:nA(i,k)=A(i,k)/A(k,k);for j=k+1:nA(i,j)=A(i,j)-

10、A(i,k)*A(k,j);endB(i)=B(i)-A(i,k)*B(k);endendA%-回代过程-x=zeros(1,n);x(n)=B(n)/A(n,n);for k=n-1:-1:1S=B(k);for j=k+1:nS=S-A(k,j)*x(j);endx(k)=S/A(k,k);end5.3.25.3.2 压缩带状对角方程组压缩带状对角方程组（1）dat53.dat 主程序:%压缩格式 dat143.datclear;format short%读出数据fid=fopen(E:2015 上机题目实验五dat53.dat,r);-20-D=fread(fid,5,long,0);n

11、=D(3);%矩阵维数p=D(5);%上带宽q=D(4);%下带宽A=fread(fid,p+q+1,n,float);B=fread(fid,n,1,float);C=A;D=B;%采用求解压缩型的列主元 Gauss 消去法x=GAUSSPP1(C,D,n,p,q);（2）针对压缩格式的带状矩阵，尤其是大型矩阵的列主元 Gauss 消去法子程序GAUSSPP1 如下：%针对压缩格式的带状方程组的列主元 Gauss 消去法function x=GAUSSPP1(A,B,n,p,q)%-消去过程-for k=1:n-1pmax=k+p;if pmax npmax=n;endfor i=k+1:p

12、maxA(i,k-i+p+1)=A(i,k-i+p+1)/A(k,p+1);qmax=k+q;if qmax nqmax=n;endfor j=k+1:qmaxA(i,j-i+p+1)=A(i,j-i+p+1)-A(i,k-i+p+1)*A(k,j-k+p+1);endB(i)=B(i)-B(k)*A(i,k-i+p+1);endend%-回代过程-x(n)=B(n)/A(n,p+1);for k=n-1:-1:1sum=0;-21-qmax=k+q;if qmax nqmax=n;endfor j=k+1:qmaxsum=sum+A(k,j-k+p+1)*x(j);endx(k)=(B(k)

13、-sum)/A(k,p+1);endend5.45.4 实验结果及分析实验结果及分析5.4.1Matlab5.4.1Matlab 运行结果运行结果方程组一：文件标识符为id=F1E1D1A0文件版本号为ver=102矩阵 A 的阶数n=15矩阵 A 的上带宽q=6矩阵 A 的下带宽p=6方程组的的解为：1.00001.00001.00001.00001.00001.00001.00001.00001.00001.00001.00001.00001.00001.00001.0000-22-方程组二：文件标识符为id=F1E1D1A0文件版本号为ver=102矩阵 A 的阶数n=2050矩阵 A

14、的上带宽q=6矩阵 A 的下带宽p=5方程组的的解为：1.96001.96001.96001.96001.96001.96001.96001.96001.96001.9600方程组三：文件标识符为文件版本号为矩阵 A 的阶数矩阵 A 的上带宽矩阵 A 的下带宽方程组的的解为：3.14133.14133.14133.14133.14133.14133.1413id=F1E1D1A0ver=202n=43512q=4p=3-23-3.14133.14133.14133.14135.4.25.4.2 总结分析总结分析对于任意一个方程组，所有根都相等，因此 Matlab 求解结果如下表所示：方程dat

15、141.datdat142.datdat143.dat阶数10205043512根1.00001.96003.1413由上面结果可以看出，所编写程序可以解决正常系数矩阵方程组，同时能解决带状系数矩阵方程组，并且在解决该类系数矩阵时具有很高的效率，说明了列主元消去法在这一类方程组的有效性。如果要继续提高计算的精度，可以采用全主元消去法。但是此时的工作量要大的多。5.55.5 本专业算例本专业算例本专业算例：一底边绝热的正方形导热物体，内无热源，其余三边温度如下，求解该正方形内节点 1、2、3、4 的温度。-24-解：导热物体无内热源，物性为常数的控制方程可以写成2T2T 0 x2y2考虑到空间步

16、长相同时，离散方程可以写成：aETEaWTWaNTNaSTSaPTP 0式中 aE=1,aW=1,aN=1,aS=1,aP=4从而有：TETWTNTS 4TP 0所有的节点方程可整理为:用列主元高斯消去法对上述矩阵进行求解，其程序如下：clear;clc;A=-4 1 1 0-70;1-4 0 1-50;1 0-3 1-15;0 1 1-3-10;m,n=size(A);for i=1:(m-1)max=i;for k=i+1:mif abs(A(k,i)abs(A(max,i)max=k;endendtemp=A(i,i:m);A(i,i:m)=A(max,i:m);A(max,i:m)=t

17、emp;for j=(i+1):mA(j,:)=A(j,:)-A(i,:)*A(j,i)/A(i,i);endenddisp(回代)x(m)=A(m,n)/A(m,m);-25-for k=(m-1):-1:1x(k)=(A(k,n)-A(k,k+1:m)*x(k+1:m)/A(k,k);endx运算结果为 T1=28.7368，T2=24.2632，T3=20.6842，T4=18.3158。运算结果和正确结果吻合。-26-学习感悟学习感悟学习是一件忙碌的事情，但主旋律是充实和快乐。在学期末，首先要感谢的是老师这学期对我们的教学。您板书的教学方式，能让我们更快的理解知识，面对面交流，会让我们对知识更深一步的掌握。编程的过程很乏味，刚入学的我也许还没有特别深入的掌握一门编程语言，所以我去借阅书籍，咨询师长，通过这样的过程，自身的编程水平也得到了一定的提升，同时也进一步掌握了我们所学的知识，这个过程让人快乐。值此元旦来临之际，祝老师节日快乐。报告中如有不足之处，还希望老师批评指正。-27-