整理西安交大计算方法b大作业.docx

1、整理西安交大计算方法b大作业西安交大计算方法b大作业 编辑整理：尊敬的读者朋友们：这里是精品文档编辑中心，本文档内容是由我和我的同事精心编辑整理后发布的，发布之前我们对文中内容进行仔细校对，但是难免会有疏漏的地方，但是任然希望（西安交大计算方法b大作业）的内容能够给您的工作和学习带来便利。同时也真诚的希望收到您的建议和反馈，这将是我们进步的源泉，前进的动力。本文可编辑可修改，如果觉得对您有帮助请收藏以便随时查阅，最后祝您生活愉快 业绩进步，以下为西安交大计算方法b大作业的全部内容。计算方法B上机实验报告学院: 机械工程学院 班级： 姓名： 学号： 2015年12月22日1。计算以下和式：，要求

2、:(1）若保留11个有效数字,给出计算结果，并评价计算的算法；（2）若要保留30个有效数字，则又将如何进行计算。实现思想：以上问题出现了近似数相减的问题，为了减小误差,可分别求得减数之和以及被减数之和，最后将两者相减.另外，减数与被减数求和均为同号计算,按照绝对值递增顺序相加可减小舍入误差。此题中对有效数字有要求,因而计算时首先需要根据有效数字位数计算得出迭代次数，以保证计算值的精度。源程序：m=input(输入有效数字个数m=)；s0=1；s1=0；s2=0;n=0；判断迭代次数while s0=0。510（m-1）s0=4/(16n*（8n+1）)-2/(16n*（8n+4）)1/(16n

3、（8n+5）)-1/（16n（8*n+6）)； n=n+1；end%分别求解各项并求和for k=n-1:1：0 a1=4/（16k*（8*k+1); a2=2/(16k*（8*k+4）); a3=1/(16k*（8k+5）)； a4=1/(16k*(8*k+6)）； s1=a1+s1; s2=a4+a3+a2+s2；endS=vpa(s1s2,m）实验结果：11位有效数字计算结果如图1所示；30为有效数字计算结果如图2所示。 图1.11位有效数字计算结果 图2.30为有效数字计算结果1.某通信公司在一次施工中,需要在水面宽度为20米的河沟底部沿直线走向铺设一条沟底光缆。在铺设光缆之前需要对沟

4、底的地形进行初步探测,从而估计所需光缆的长度,为工程预算提供依据。已探测到一组等分点位置的深度数据(单位：米）如下表所示：分点0123456深度9。018.967。967。978.029。0510.13分点78910111213深度11.1812。2613。2813.3212。6111。2910.22分点14151617181920深度9。157.907。958。869。8110.8010。93 （1）请用合适的曲线拟合所测数据点；（2)估算所需光缆长度的近似值，并作出铺设河底光缆的曲线图；算法思想：由于题中所给点数为20，若采用高次多项式插值将产生很大的误差，所以拉格朗日或牛顿并不适用。题中

5、光缆为柔性，可光滑铺设于水底，鉴于此特性，采用三次样条插值插值法较为合适.算法结构：三次样条算法结构见计算方法教程P110;光缆长度计算公式：源程序：clear；clc; x=0:20； y=9。01 8.96 7.96 7。97 8。02 9.05 10。13 11.18 12。26 13。28 13.32 12。61 11.29 10.22 9。15 7.90 7.95 8.86 9.81 10。80 10。93; d=y； plot（x,y，k.,markersize,15） hold on %计算差商 for k=1：2 for i=21:-1：（k+1) d(i）=(d(i）d（i1

6、)）/(x（i)-x(ik); end end%设定d的边界条件 for i=2：20 d（i)=6d（i+1）; end d（1)=0; d（21)=0; %带状矩阵求解（追赶法）a=0。5ones(1，21）；b=2ones(1,21)；c=0。5*ones(1,21）；a（1)=0；c(21）=0；u=ones(1，21）；u（1）=b(1）；r=c；yy（1)=d(1)； %追 for k=2:21 l（k)=a（k)/u(k1)； u（k）=b（k）l(k）r(k-1); yy（k）=d（k)l（k）yy（k-1); end %赶 m（21)=yy（21)/u（21）; for k=

7、20:-1:1 m(k)=(yy(k）-r（k）m(k+1）)/u(k); end %绘制曲线 k=1； nn=100; xx=linspace(0，20，nn); l=0； for j=1:nn for i=2：20 if xx(j)=x(i) k=i； break； else k=i+1； endendh=1；xbar=x（k）xx（j);xmao=xx（j）x(k-1）；s（j)=（m（k1)xbar3/6+m(k）*xmao3/6+(y(k-1）m（k-1）h2/6)xbar+(y（k）-m（k)h2/6）*xmao)/h;sp(j)=-m(k1）*（x（k)xx(j）2/（2*h）+

8、m(k)（xx（j)-x(k1)2/（2h)+（y(k)-y(k1）/h(m(k）m（k1)）h/6;l(j+1）=（1+sp（j）2）0。5(20/nn)+l(j）; %求解光缆长度end%绘图plot(xx，s，r-,linewidth,1.5)disp(光缆长度为，num2str（l（nn+1）),）曲线图如图2-1所示，计算光缆长度如图2-2所示。图2-1光缆插值曲线图图21光缆计算长度显示3。假定某天的气温变化记录如下表所示，试用数据拟合的方法找出这一天的气温变化的规律；试计算这一天的平均气温,并试估计误差。时刻0123456789101112平均气温1514141414151618

9、2020232528时刻131415161718192021222324平均气温313431292725242220181716实现思想：此题中所给数据点数目较多，采用拉格朗日插值法或者牛顿插值法需要很高次的多项式,计算困难，误差大；采用样条插值计算量虽然不大，但是存放参数Mi的量很大，且没有一个统一的数学公式来表示，也不是很方便。所以可考虑用最小二乘法进行拟合。计算过程中，分别使用二次函数、三次函数以及四次函数，计算其相应的系数，估算误差并作图比较各个函数之间的区别。算法结构：（参考课本P123）1.1形成矩阵Qk1.2变换Gk1到Gk2.求解三角方程3。计算误差源代码:clear；clc;

10、x=0：24;y=15141414141516182020232528313431292725242220181716;m=length(x）;n=input(请输入函数的次数)；plot（x,y，k。,x，y，-）grid;holdon;n=n+1；G=zeros（m,n+1)；G(:,n+1）=y；c=zeros(1,n);建立c来存放q=0;f=0;b=zeros(1,m)；建立b用来存放形成矩阵Gforj=1：nfori=1：mG（i，j）=x（1,i）(j-1）;endend%建立矩阵Qkfork=1：nfori=k：mc(k）=G(i，k）2+c(k)；endc(k)=sign（G

11、(k，k）)*（c（k）0.5）;w(k)=G（k，k）-c(k）；%建立w来存放forj=k+1:mw（j）=G（j，k）;endb(k)=c（k）*w(k);%变换矩阵Gk-1到GkG（k，k）=c(k)；forj=k+1:n+1q=0；fori=k：mq=w(i）G(i，j）+q;ends=q/b(k)；fori=k：mG(i,j)=s*w(i)+G(i，j）；endendend%求解三角方程Rx=h1a(n)=G(n,n+1)/G（n，n)；fori=n1:(-1）:1forj=i+1：nf=G（i，j）a（j)+f；enda(i)=（G（i，n+1)f)/G(i,i)；%a（i)存放

12、各级系数f=0；enda%回代过程p=zeros(1，m)；forj=1：mfori=1:np（j)=p(j)+a（i）x（j）(i1）；endendplot(x,p，r，x,p,);E2=0;%用E2来存放误差%误差求解fori=n+1：mE2=G（i，n+1）2+E2；endE2=E20。5;disp(误差为）;disp（E2）；t=0;fori=1：mt=t+p(i）;endt=t/m；%平均温度disp(平均温度为,num2str(t)，）实验结果：二次函数拟合，结果如下图所示图3-1 二次函数拟合结果三次函数拟合，结果如下图所示图3-2 三次函数拟合结果四次函数拟合，结果如下图所示图

13、33 四次函数拟合结果结果对比:将二次函数、三次函数和四次函数拟合结果绘制在同一个坐标内,如图34所示。其计算误差结果见表3-1所示。图3-4 拟合结果对比分析4。设计算法，求出非线性方程的所有实根，并使误差不超过。算法思想：本题可采用牛顿法迭代求解,令 ,得带格式为根据函数图像可以找出根的大致分布区间,带入不同的初值即可解出不同的根.源代码:function y=f2（x)y=6*x。545x.2+20；定义原函数function y=f3(x）y=30x490*x;定义原函数倒数i=5：0。1：5；y=f2(i）；plot(i,y）hold onplot(i，0,-）%画出原函数图像%Ne

14、wton法求根x1=input（输入初值);e=10（-4)；%误差设定Nmax=1000;迭代最大次数限定for n=1：Nmax f0=f2（x1）; if abs(f2（x1)e fprintf(输出的f(x)已经足够小); x=x1; break else F0=f3（x1)； x=x1f0/F0; if abs（x-x1)e break else x1=x; end endendfprintf（输出方程的根x=2f，x) 计算结果：函数图像如图4-1所示。计算结果分别见图42所示。图4-1 函数图像图42 计算结果根据带入不同的初值，可以求出不同的根，有图42可以看出，原函数的根大约

15、有三个,分别是0.654542、0。681174、1.870799。5。线性方程组求解.(1)编写程序实现大规模方程组的高斯消去法程序，并对所附的方程组进行求解。所附方程组的类型为对角占优的带状方程组。（2）针对本专业中所碰到的实际问题，提炼一个使用方程组进行求解的例子，并对求解过程进行分析、求解。算法思想:高斯消去法是利用现行方程组初等变换中的一种变换,将一个不为零的数乘到一个方程后加到另一个方程,使方程组变成同解的上三角方程组,然后再自下而上对上三角方程组求解.算法结构:1.读取二进制文件，存入计算矩阵2.对矩阵进行初等变换，然后求解(详见计算方法教程第2版高斯消去法以及列主元高斯消去法算

16、法）源代码：clear； clc； % 读取系数矩阵f，p=uigetfile（*.dat，选择数据文件）； %读取数据文件num=5; 输入系数矩阵文件头的个数name=strcat(p,f）;file=fopen(name，r)；head=fread(file，num，uint）； %读取二进制头文件id=dec2hex（head（1)）; %读取标识符fprintf(文件标识符为）；idver=dec2hex(head(2)； 读取版本号fprintf（文件版本号为）;vern=head(3）； %读取阶数fprintf(矩阵A的阶数）；nq=head（4); %上带宽fprintf（矩

17、阵A的上带宽）;q p=head(5）； 下带宽fprintf(矩阵A的下带宽）；p dist=4*num；fseek(file,dist,bof）; %把句柄值转向第六个元素开头处A,count=fread(file,inf,float）; 读取二进制文件，获取系数矩阵fclose（file）; %关闭二进制头文件% 对非压缩带状矩阵进行求解if ver=102, a=zeros(n,n)； for i=1：n， for j=1：n, a(i,j）=A(（i1)n+j）； 求系数矩阵a（i，j） end end b=zeros(n,1); for i=1：n, b(i)=A(nn+i）; e

18、nd for k=1：n-1， 列主元高斯消去法 m=k; for i=k+1：n, 寻找主元 if abs（a（m,k)abs(a（i,k）) m=i; end end if a(m,k)=0 %遇到条件终止 disp（错误！） return end for j=1：n， %交换元素位置得主元 t=a（k,j）； a（k,j)=a（m，j）; a(m,j）=t; t=b（k); b（k)=b(m）； b（m）=t; end for i=k+1:n, %计算l(i，k)并将其放到a(i,k)中 a(i，k）=a（i,k）/a（k，k）； for j=k+1:n a（i，j）=a(i,j)a(i

19、，k）*a(k，j); end b(i）=b(i）-a（i,k）b(k）； end end x=zeros（n,1)； 回代过程 x(n）=b（n）/a(n，n）; for k=n1：1：1， x(k）=（b(k）sum(a（k,k+1：n）*x（k+1：n)）/a(k,k）; endend% 对压缩带状矩阵进行求解if ver=202, 高斯消去法 m=p+q+1; a=zeros（n,m)； for i=1:1：n for j=1:1：m a（i，j）=A(（i1)m+j）； end end b=zeros（n，1）; for i=1：1：n b(i）=A(n*m+i)； %求b(i) e

20、nd for k=1：1:（n1) %开始消去 if a(k，(p+1）=0 disp(错误!）； break; end st1=n; if （k+p)n st1=k+p; end for i=（k+1)：1：st1 a（i，（k+pi+1)）=a(i，(k+pi+1)/a（k，(p+1）)； for j=（k+1)：1：（k+q） a（i,j+pi+1)=a（i，j+pi+1)-a(i,k+pi+1）a（k,j+pk+1）； end b（i）=b（i）a(i,k+pi+1）*b(k); end end x=zeros（n，1); 回代 x（n)=b(n）/a(n,p+1）； sum=0; f

21、or k=(n-1)：-1：1 sum=b(k）； st2=n; if （k+q）n st2=k+q; end for j=（k+1):1:st2 sum=suma(k，j+pk+1)x(j）; end x(k）=sum/a（k，p+1）; sum=0； endenddisp(方程组的的解为：） 输出解disp(x)求解结果对数据文件dat51求解,结果如下：文件标识符为id = F1E1D1A0文件版本号为ver = 102矩阵A的阶数n = 15矩阵A的上带宽q = 3矩阵A的下带宽p = 3方程组的的解为：1.0000 1。0000 1。0000 1.0000 1.0000 1。0000

22、 1。0000 1.0000 1。0000 1.0000 1。0000 1。0000 1.0000 1。0000 1。0000对数据文件dat52求解，结果如下：文件标识符为id = F1E1D1A0文件版本号为ver = 202矩阵A的阶数n = 20矩阵A的上带宽q = 5矩阵A的下带宽p = 5方程组的的解为：1.0000 1。0000 1.0000 1。0000 1.0000 1。0000 1.0000 1。0000 1.0000 1.0000 1.0000 1.0000 1.0000 1.0000 1.0000 1.0000 1。0000 1。0000 1.0000 1.0000对数据文件dat53求解,结果如下：文件标识符为id = F1E1D1A0文件版本号为ver = 102矩阵A的阶数n = 2160矩阵A的上带宽q = 5矩阵A的下带宽p = 5方程组的解截图如图51所示（由于矩阵阶数较大，计算结果未能截取完整)。图5-1数据文件dat53求解结果对数据文件dat53求解,结果如下：文件标识符为id = F1E1D1A0文件版本号为ver = 202矩阵A的阶数n = 43240矩阵A的上带宽q = 4矩阵A的下带宽p = 4方程组的解截图如图52所示（由于矩阵阶数较大，计算结果未能截取完整)。图52 数据文件dat54求解结果