计算方法上机作业(20210914024420)

1、手及交通丈李XrAMJIAOTOFKUNIVERSITY计算方法上机报告姓名:学号:班级:上课班级：说明:MATLAB本次上机实验使用的编程语言是Matlab语言，编译环境为7.11.0,运行平台为Windows7。14211一，、1.对以下和式计算：o凝8n18n48n58n6,要求：若只需保留11个有效数字，该如何进行计算；若要保留30个有效数字，则又将如何进行计算；(1)算法思想1、根据精度要求估计所加的项数，可以使用后验误差估计，通项为:1421114an-n-n-16n8n18n48n58n616n8n12、为了保证计算结果的准确性，写程序时，从后向前计算;3、使用Matlab时，可

2、以使用以下函数控制位数：digits(位数)或vpa(变量，精度为数)(2)算法结构1. s0;1421116n8n18n48n58n6if2. forn0,1,2,i10end;3. forni,i1,i2,0sst;(3)Matlab源程序clear;clc;m=input(,请输入有效数字的位数m=');s=0;forn=0:50t=(1/16An)*(4/(8*n+1)-2/(8*n+4)-1/(8*n+5)-1/(8*n+6);ift<=10A(-m)break;endend;fprintf('需要将n值加到n=%dn',n-1);fori=n-1:-1

3、:0t=(1/16Ai)*(4/(8*i+1)-2/(8*i+4)-1/(8*i+5)-1/(8*i+6);s=s+t;ends=vpa(s,m)(4)结果与分析当保留11位有效数字时，需要将n值加到s=3.1415926536当保留30位有效数字时，需要将n值加到%青除工作空间变量%青除命令窗口命令%俞入有效数字的位数喏U断通项与精度的关系%!要将n值加到的数值%求和运算%$制s的精度n=7,；n=22,o这种算法很好的保证了计算结20米的河沟底部沿直线s=3.14159265358979323846264338328通过上面的实验结果可以看出，通过从后往前计算，果要求保留的准确数字位数的要

4、求。2.某通信公司在一次施工中，需要在水面宽度为走向铺设一条沟底光缆。在铺设光缆之前需要对沟底的地形进行初步探测，从而估计所需光缆的长度，为工程预算提供依据。已探测到一组等分点位置的深度数据（单位：米）如下表所示：分点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)算法思想如果使用多

5、项式差值，则由于龙格现象，误差较大，因此，用相对较少的插值数据点作插值，可以避免大的误差，但是如果又希望将所得数据点都用上，且所用数据点越多越好，可以采用分段插值方式，即用分段多项式代替单个多项式作插值。分段多项式是由一些在相互连接的区间上的不同多项式连接而成的一条连续曲线，其中三次样条插值方法是一种具有较好“光滑性”的分段插值方法。在本题中，假设所铺设的光缆足够柔软，在铺设过程中光缆触地走势光滑，紧贴地面，并且忽略水流对光缆的冲击。海底光缆线的长度预测模型如下所示，光缆从A点铺至B点，在某点处的深度为hi。计算光缆长度时，用如下公式:20L0f(x)ds_'2f(x)dx019k1f

6、(x).1f'(x)2dxy)2(2)算法结构1.For0,1,2,2.3.4.5.6.7.8.9.1.1For2.1For4.14.24.3doVFor2.1.1x0hiMi1,2in,n1,k(MiMi1)/(xixik)Mi1,2,6Mi1M0;dnb。；n1,d1,n-1diMn；an;2Ci;1ca;2bbn获取M的矩阵元素个数，存入mFor8.18.28.3m/k2,3,mak/k1lkbk-lkCk1dk-lkk1mMm10. Forkm1,m2,110.1 (kckMki)/kMk11. 获取x的元素个数存入s12. 1k13. Fori1,2,s113.1 ifxXi

7、thenik；breakelsei1kh一八;XkXx;XXk1?X3x3h2h2Mk*Mk否(yk1Mk1x(ykMk石阂/h(3)Matlab源程序clear;clc;x=0:1:20;炉生从0到20含21个等分点的数组X=0:0.2: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;n=length(x);N=length(X);%求三次样条插值函数s(x)M=y;fork=2:3;fori=n:-1

8、:k;M(i)=(M(i)-M(i-1)/(x(i)-x(i-k+1);endendh(1)=x(2)-x(1);fori=2:n-1;h(i)=x(i+1)-x(i);c(i)=h(i)/(h(i)+h(i-1);a(i)=1-c(i);b(i)=2;d(i)=6*M(i+1);endM(1)=0;M(n)=0;嗡分点位置的深度数据%等分点的数目%十算二阶差商并存放在W%十算三对角阵系数a,b,c及右端向量d处择自然边界条件b(1)=2;b(n)=2;c(1)=0；a(n)=0;d(1)=0;d(n)=0;u(1)=b(1);y1(1)=d(1);fork=2:n;l(k)=a(k)/u(k

9、-1);u(k)=b(k)-l(k)*c(k-1);y1(k)=d(k)-l(k)*y1(k-1);endM(n)=y1(n)/u(n);fork=n-1:-1:1;M(k)=(y1(k)-c(k)*M(k+1)/u(k);ends=zeros(1,N);form=1:N;k=1;fori=2:n-1ifX(m)<=x(i);k=i-1;break;else%寸三对角阵进行LU分解%追赶法求解样条参数M(i)k=i;endendH=x(k+1)-x(k);%B各区间用三次样条插值函数计算X自处的值x1=x(k+1)-X(m);x2=X(m)-x(k);s(m)=(M(k)*(x1A3)/

10、6+M(k+1)*(x2A3)/6+(y(k)-(M(k)*(HA2)/6)*x1+(y(k+1)-(M(k+1)*(HA2)/6)*x2)/H;end%计算所需光缆长度L=0;fori=2:N%十算所需光缆长度L=L+sqrt(X(i)-X(i-1)A2+(s(i)-s(i-1)A2);enddisp('所需光缆长度为L=');disp(L);figureplot(x,y,'*',X,s,'-')线图%绘制铺设河底光缆的曲痴注坐标轴含义xlabel('位置，'fontsize',16);ylabel('深度/m&

11、#39;,'fontsize',16);title('铺设河底光缆的曲线图'，'fontsize',16);grid;(4)结果与分析铺设海底光缆的曲线图如下图所示：仿真结果表明，运用分段三次样条插值所得的拟合曲线能较准确地反映铺设光缆的走势图，计算出所需光缆的长度为L=26.4844m。3.假定某天的气温变化记录如下表所示，试用数据拟合的方法找出这一天的气温变化的规律；试计算这一天的平均气温，并试估计误差。时刻0123456789101112平均气温15141414141516182020232528时刻13141516171819202122

12、2324平均气温313431292725242220181716(1)算法思想在本题中，数据点的数目较多。当数据点的数目很多时，用“多项式插值”方法做数据近似要用较高次的多项式，这不仅给计算带来困难，更主要的缺点是误差很大。用“插值样条函数”做数据近似，虽然有很好的数值性质，且计算量也不大，但存放参数Mi的量很大，且没有一个统一的数学公式来表示，也带来了一些不便。另一方面，在有的实际问题中，用插值方法并不合适。当数据点的数目很大时，要求P(x)通过所有数据点，可能会失去原数据所表示的规律。如果数据点是由测量而来的，必然带有误差，插值法要求准确通过这些不准确的数据点是不合适的。在这种情况下，不用

13、插值标准而用其他近似标准更加合理。通常情况下，是选取i使E2最小，这就是最小二乘近似问题。在本题中，采用“最小二乘法”找出这一天的气温变化的规律，使用二次函数、三、2次函数、四次函数以及指数型函数Caeb(tc),计算相应的系数，估算误差，并作图比较各种函数之间的区别。(2)算法结构本算法用正交化方法求数据的最小二乘近似。假定数据以用来生成了g，并将y作为其最后一列(第n1列)存放。结果在a中，是误差E；。I、使用二次函数、三次函数、四次函数拟合时3.0.2 .将“时刻值”存入Xi,数据点的个数存入m4.0.2 .输入拟合多项式函数p(x)的最高项次数hn1，则拟合多项式函数为P(X)19(X

14、)2g2(X)ngn(x),根据给定数据点确定GForj0,1,2,n1Fori1,2,m4.1.2.1 Xijgi,j14.2.2.1 Vgi,n15.0.2 .Fork1,2,n5.1.2.1 形成矩阵Qk1. sgn(gkk)(gik)ik1. gkkkForjk1,kgjkk3.2变换Gk1到GkgkkForjk1,k2,m(ig"tikForik,k1,gjti.解三角方程Rah11gn.n1/gnnan.2Forin1,n2,1n21gi.n1gijxj/gii.计算误差E2m2gi,n1n1II、使用指数函数拟合时现将指数函数进行变形：2将Cy,tX代入Caeb()得:

15、对上式左右取对数得：2,m,n1,mgoaib(xc)2yae22InyInabc2bcxbx令zlny,olnabc2,12bc,则可得多项式:2zo1X2X(3)Matlab源程序clear;clc;x=0:24;%青除工作空间变量%青除命令窗口命令%各时刻值存入数组y=15,14,14,14,14,15,16,18,20,20,23,25,28,31,34,31,29,27,25,24,22,20,18,17,16;,m=size(x);T=sum(y(1:m)/m;fprintf('一天的平土!气温为T=%fn',T);%二次、三次、四次函数的最小二乘近似h=input

16、('请输入拟合多项式的最高项次数=');n=h+1;G=;forj=0:(n-1)g=x.Aj;G=vertcat(G,g);endG=G'fori=1:mG(i,n+1)=y(i);endG;fork=1:nifG(k,k)>0;sgn=1;elseifG(k,k)=0;sgn=0;elsesgn=-1;endsgm=-sgn*sqrt(sum(G(k:m,k).A2);w=zeros(1,n);w(k)=G(k,k)-sgm;forj=k+1:mw(j)=G(j,k);endbt=sgm*w(k);G(k,k)=sgm;forj=k+1:n+1t=sum(w(

17、k:m)*G(k:m,j)/bt;fori=k:m;%各数据点的个数存入m%求一天的平均气温%艮据给定数据点生成矩阵%g(x)按列排列%亚直连接G嗨专置得到矩阵G%各数据y作为GW最后一列嘛成矢I阵Q(k)射换Gk-1到Gk(n+1歹U)G(i,j)=G(i,j)+t*w(i);endendendA(n)=G(n,n+1)/G(n,n);fori=n-1:-1:1A(i)=(G(i,n+1)-sum(G(i,i+1:n).*A(i+1:n)/G(i,i);ende=sum(G(n+1:m,n+1).A2);嗨军三角方程求系数Afprintf(disp(A);fprintf(,%d次函数的系数是

18、：',h);%计算误差e%俞出系数a及误差e,使用歆函数拟合的误差是f:',h,e);t=0:0.05:24;A=fliplr(A);Y=poly2sym(A);subs(Y,'x',t);Y=double(ans);figure(1)plot(x,y,'k*',t,Y,'r-');xlabel('时刻);ylabel('平均气温');title(num2str(n-1),'次函数的最小二乘曲线');grid;%指数函数的最小二乘近似yy=log(y);n=3;G=;GG=;forj=0:(

19、n-1)g=x.Aj;G=vertcat(G,g);gg=t.Aj;GG=vertcat(GG,gg);endG=G'fori=1:mG(i,n+1)=yy(i);%各系数数组左右翻转%各系数数组转化为多项式侬制拟合多项式函数图形两注坐标轴含义%g(x)按列排列%亚直连接G%g(x)按列排列%我直连接G嗨专置得到矩阵G%各数据y作为GW最后一列(n+1歹U)endG;fork=1:n名形成矢I阵Q(k)ifG(k,k)>0;sgn=1;elseifG(k,k)=0;sgn=0;elsesgn=-1;endsgm=-sgn*sqrt(sum(G(k:m,k).A2);w=zeros

20、(1,n);w(k)=G(k,k)-sgm;forj=k+1:mw(j)=G(j,k);endbt=sgm*w(k);G(k,k)=sgm;forj=k+1:n+1t=sum(w(k:m)*G(k:m,j)/bt;fori=k:m;G(i,j)=G(i,j)+t*w(i);endendendA(n)=G(n,n+1)/G(n,n);fori=n-1:-1:1A(i)=(G(i,n+1)-sum(G(i,i+1:n).*A(i+1:n)/G(i,i);endb=-A(3);c=A(2)/(2*b);a=exp(A(1)+b*(cA2);G(n+1:m,n+1)=exp(sum(G(n+1:m,n

21、+1).A2);e=sum(G(n+1:m,n+1).A2);fprintf('n指数函数的系数是：a=%f,b=%f,c=%ffprintf('n使用指数函数拟合的误差是：f,e);t=0:0.05:24;YY=a.*exp(-b.*(t-c).A2);figure(2)plot(x,y,'k*',t,YY,'r-');xlabel('时刻);ylabel('平均气温');title('指数函数的最小二乘曲线');grid;侬换Gk-1到Gk%由三角方程求系数A%十算误差e,a,b,c);%俞出系数及误差

22、e%绘制拟合指数函数图形%标注坐标轴含义(4)结果与分析a、二次函数：一天的平均气温为：21.20002次函数的系数：8.30632.6064-0.0938使用2次函数拟合的误差是：280.339547.下载可编辑.二次函数的最小二乘曲线如下图所示:2次曲船为品子二暴曲孤20353O2SM佑to明lr£前b、三次函数:一天的平均气温为：21.20003次函数的系数：13.3880-0.22730.2075-0.0084使用3次函数拟合的误差是：131.061822三次函数的最小二乘曲线如下图所示：0ia3次瓯故的昂小二桀曲第c、四次函数：一天的平均气温为：21.20004次函数的系数

23、：16.7939-3.70500.8909-0.05320.0009使用4次函数拟合的误差是：59.04118四次函数的最小二乘曲线如下图所示：1生id05ID152025时荆4次打毁的我匚、乘曲绿5P2Zd、指数函数：一天的平均气温为：21.2000指数函数的系数是：a=26.160286,b=0.004442,c=14.081900使用指数函数拟合的误差是：57.034644指数函数的最小二乘曲线如下图所示：指却由荆刘易曾二乘曲姓L+*一*，、/7:/+*LAi*4y*/F1011110610IS2025时邦通过上述几种拟合可以发现，多项式的次数越高，计算拟合的效果越好，误差越小，说明结果

24、越准确；同时，指数多项式拟合的次数虽然不高，但误差最小，说明结果最准确。524.设计算法，求出非线性方程6x45x200的所有根，并使误差不超过10(1)算法思想首先，研究函数的形态，确定根的范围；通过剖分区间的方法确定根的位置，然后利用二分法的基本原理进行求解，找到满足精度要求的解。x,x(k°确定区间严1增勺二分法是产生一审区间，使新区间是旧区间I(k)的一个子区间，其长度是I(k)的一半，且有一个端点是|(k)的一个端点。由区间I的方法是计算区间I(k)的中点(k2)x-(x(k)x(k1)2(k)(k2)(k1)(k)(k2)(k1)(k2)(k1)右f(x)f(x)U,则取

25、Ix,x，畲则取Ix,x，生旻这也程即可。显然，每次迭代使区间长度减小一半，故二分法总是收敛的。(2)算法结构f(x(0)fo；f(x(1)flIff0f10thenstopIf|f0|2then输出x(0)作为根；stopIf|f1|2then输出x作为根；stop2(x(0)x(1)xIf|x(1)x|1|x(1)|then输出x作为根；stopf(x)fIf|f|2then输出x作为根；Iff1f0thenxx(0)；ff。else(1)xx;ff1goto5(3)Matlab源程序x=-100:100;线性方程组的表达式%确定根所在的区间%区间长度为1嘛定根的个数y=6*(x.A5)-

26、45*(x.A2)+20;g=;fori=-100:1:100k=i+1;if(y(x=i).*y(x=k)<eps)g=gi;endendsymsx;f=6*xA5-45*xA2+20;n=length(g);forj=1:nx0=g(j);%求根区间左端点x1=g(j)+1;9家根区间右端点while(x1-x0)>=10A(-4)ifsubs(f,x,x0)*subs(f,x,(x0+x1)/2)>epsx0=(x0+x1)/2;elsex1=(x0+x1)/2;endendroot=x0喻出方程的根end(4)结果与分析该非线性方程组有三个实根，分别为1.8708,0

27、.6812,-0.6545,且满足误差要求。5.编写程序实现大规模方程组的列主元高斯消去法程序，并对所附的方程组进行求解。针对本专业中所碰到的实际问题，提炼一个使用方程组进行求解的例子，并对求解过程进行分析、求解。(1)算法思想高斯消去法是利用现行方程组初等变换中的一种变换，即用一个不为零的数乘一个方程后加只另一个方程，使方程组变成同解的上三角方程组，然后再自下而上对上三角方程组求解。列主元消去法是当高斯消元到第k步时，从k列的akk以下(包括akk)的各元素中选出绝对值最大的，然后通过行交换将其交换到akk的位置上。交换系数矩阵中的两行(包括常数项)，只相当于两个方程的位置交换了，因此，列选

28、主元不影响求解的结果。程序的核心就是高斯列主元消去法。根据教材提供的算法，编写列主元消去法的子函数与适应于超大规模超出系统内存的方程组的改编程序。同时，在Gaussffi去过程中，适当交换方程的顺序对保证消去过程能顺利进行及计算解的精确度都是有必要的，交换方程的原则是使f1%k,k1,n)中，绝对值最大的一个换到(k,k)位置而成为第k步消去的主元，这就是列主元Gauss消去法。(2)算法结构1、数据文件的文件名为：文件名+.dat2、数据文件中的数据为二进制记录结构，分为以下四个部分：(1)文件头部分，其结构：typedefstruct(longintid;longintver;longin

29、tn;其中：id:为该数据文件的标识，值为0xF1E1D1A0即为：十六进制的F1E1D1A0ver:为数据文件的版本号，值为16进制数据，版本号说明0x101系数矩阵为非压缩格式稀疏矩阵0x102系数矩阵为非压缩格式带状对角阵0x201系数矩阵为压缩格式稀疏矩阵0x202系数矩阵为压缩格式带状对角阵n:表示方程的阶数(2)文件头2:此部分说明为条状矩阵的上下带宽，结构：typedefstruct(longintq;/为上带宽longintp;/为下带宽(3)系数矩阵a.如存贮格式非为压缩方式，则按行方式存贮系数矩阵中的每一个元素，个数为n*n,类型为float型；b.如果存贮格式是压缩方式，

30、则按行方式存贮，每行中只存放上下带宽内的非零元素，即，每行中存贮的最多元素为p+q+1个。(4)右端系数按顺序存贮右端系数的每个元素，个数为n个，类型为float型3、二进制文件的读取：f=fopen('fun003.dat','r');%打开文件，.dat文件放在m文件同一目录下,a=fread(f,3,'uint')%读取头文件，3-读取前3个，若读取压缩格式的，头文件为5个b=fread(f,inf,'float');%读取剩下白文件，float型id=dec2hex(a(1);ver=dec2hex(a(2);%这两句是进

31、行进制转换，读取id与verA的阶数nFork1,2,3,n12.1找满足kkmaxik的下标kForj1,2,n2.2.1kjkjk2.4Forik1,k2,n2.4.1ikikkk2.4.2Forjk1,k2,n2.4.2.1ijikkjij2.4.3i-ikkin/nnxnForkn1,n2,1n(kkjxj)/kkxk(3)Matlab源程序clear;clc;%读取系数矩阵f,p=uigetfile('*.dat','选择数据文件);num=5;name=strcat(p,f);file=fopen(name,'r');head=fread(f

32、ile,num,'uint');id=dec2hex(head(1);fprintf('文件标识符为)；idver=dec2hex(head(2);fprintf('文件版本号为)；vern=head(3);fprintf('矩阵A勺阶数');nq=head(4);fprintf('矩阵A勺上带宽');qp=head(5);fprintf('矩阵A勺下带宽');pdist=4*num;fseek(file,dist,'bof);头处A,count=fread(file,inf,'float'

33、);矩阵fclose(file);%对非压缩带状矩阵进行求解ifver='102',a=zeros(n,n);fori=1:n,forj=1:n,a(i,j)=A(i-1)*n+j);endendb=zeros(n,1);fori=1:n,b(i)=A(n*n+i);end%青除工作空间变量%青除命令窗口命令%卖取数据文件%俞入系数矩阵文件头的个数取二进制头文件%取标识符%卖取版本号%卖取阶数%上带宽%r带宽叫巴句柄值转向第六个元素开取二进制文件，获取系数%关闭二进制头文件%求系数矩阵a(i,j)%主元高斯消去法fork=1:n-1,m=k;fori=k+1:n,ifabs(a

34、(m,k)<abs(a(i,k)m=i;endendifa(m,k)=0disp('错误！)returnendforj=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;endfori=k+1:n,中a(i,k)=a(i,k)/a(k,k);forj=k+1:na(i,j)=a(i,j)-a(i,k)*a(k,j);endb(i)=b(i)-a(i,k)*b(k);endendx=zeros(n,1);x(n)=b(n)/a(n,n);fork=n-1:-1:1,x(k)=(b(k)-sum(a(k,k+1:n

35、)*x(k+1:n)/a(k,k);endend%对压缩带状矩阵进行求解ifver='202',m=p+q+1;a=zeros(n,m);fori=1:1:nforj=1:1:ma(i,j)=A(i-1)*m+j);endendb=zeros(n,1);fori=1:1:n%寻找主元%遇到条件终止眇换元素位置得主元%十算l(i,k)并将其放到a(i,k)刎代过程喻斯消去法麻a(i,j)麻b(i)%开始消去过程b(i)=A(n*m+i);endfork=1:1:(n-1)ifa(k,(p+1)=0disp('错误！)；break;endst1=n;if(k+p)<nst1=k+p;endfori=(k+1):1:st1a(i,(k+p-i+1)=a(i,(k+p-i+1)/a(k,(p+1);forj=(k+1):1:(k+q)a(i,j+p-i+1)=a(i,j+p-i+1)-a(i,k+p-i+1)*a(k,j+p-k+1);endb(i)=b(i)-a(i,k+p-i+1)*b(k);end%回代过程endx=zeros(n,1);x(n)=b(n)/a