计算方法上机作业.docx

《计算方法上机作业.docx》由会员分享，可在线阅读，更多相关《计算方法上机作业.docx（70页珍藏版）》请在课桌文档上搜索。

1、计算方法上机报告姓名：学号：班级：上课班级：说明：本次上机实验使用的编程语言是Matlab语言，编译环境为MATLAB1.11.0,运行平台为WindOWS7。1.对以下和式计算：s=y(-i?1Zj16n8n18n+48n+58n+6)TTTr八，要求：若只需保留11个有效数字，该如何进行计算；若要保留30个有效数字，则又将如何进行计算；(1)算法思想1、根据精度要求估计所加的项数，可以使用后验误差估计，通项为:(211148n+48n+5-8n+6Jcc16n8n+lf2、为了保证计算结果的准确性，写程序时，从后向前计算；3、使用MatIab时，可以使用以下函数控制位数：digits（位数

2、）或vpa（变量，精度为数）（2）算法结构t=_L（_l21i_1s-0;16n8n+18n+48n+58n+6/12. for11-0,1,2,Jif10end;3. for11=j-lj-2,0s-5t;O)MatIab源程序clear;%清除工作空间变量clc;%清除命令窗口命令m=input(谛输入有效数字的位数m三,);%输入有E数字的位数s=0;forn=0:50t=(l16n)*(4(8*n+l)-2(8*n+4)-l(8*n+5)-l(8*n+6);ift4/)4*M1不+如石+(JI-MI石)x+(A-Mlrr)x/力=yMatlab源程序clear;clc;x=0:l:20

3、;%产生从0到20含21个等分点的数组X=O: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=Iength(X);%求三次样条插值函数S(X)M=y;fork=2:3;%计算二阶差商并存放在M中fori=n:-l:k;M(i)=(M(i)-M(i-l)(i)-x(i-k+l);endendh(l)=(2)-(l);%计算三对角阵系与a

4、,b,c及右端向量dfori=2:n-l;h(i)=x(i+l)-x(i);c(i)=h(i)(h(i)+h(i-l);a()=l-c(i);b(i)=2;d(i)=6*M(i+l);endM(I)=O;%选择自然边界条件M(n)=0;b(l)=2;b(n)=2;c(D=0;a(n)=0;d(D=0;d(n)=0;u(l)=b(l);%对三对角阵进行LU分解yl(l)=d(l);fork=2:n;l(k)=a(k)u(k-l);u(k)=b(k)-l(k)*c(k-l);yl(k)=d(k)-l(k)*yl(k-l);endM(n)=yl(n)u(n);%追赶法求解样条参数M(i)fork=n

5、-l:-l:l;M(k)=(yl(k)-c(k)*M(k+l)u(k);ends=zeros(l,N);form=l:N;k=l;fori=2:n-lifX(m)0;sgn=l;elseifG(kzk)=;sgn=0;elsesgn=-l;endsgm=-sgn*sqrt(sum(G(k:m/k).A2);w=zeros(l,n);w(k)=G(k,k)-sgm;forj=k+l:mw(j)=G(j,k);endbt=sgm*w(k);G(k,k)=sgm;%变换Gk-I到Gkforj=k+l:n+lt=sum(w(k:m)*G(k:m,j)/bt;fori=k:m;G(i,j)=G(i,j)

6、+t*w(i);endendendA(n)=G(n,n+l)G(n,n);%解三角方程求系数Afori=n-l:-l:lA(i)=G(i,n+l)-sum(G(iJ+lrn).*A(i+l:n)/G(iJ);ende=sum(G(n+lmzn+l).2);%计算误差efprintf(,%d次函数的系数是：,h);%输出系数a及idisp（八）;fprintf(使用d次函数拟合的误差是f：,h,e);t=0:0.05:24;A=fliplr（八）;%将系数数组左右翻转Y=poly2sym（八）;%将系数数组转化为多项式SUbS(Yx,t);Y=double(ans);figure(l)plot(

7、,yk*X,r-);%绘制拟合多项式函数图形XlabeIr时期;%标注坐标轴含义Wabel。平均气温上title(wm2str(nl)J次函数的最小二乘曲线grid;%指数函数的最小二乘近似y=og(y);n=3;G=;GG=U;forj=O:(n-l)g=.j;%g(x)按列排列G=vertcat(Gzg);%g垂直连接Ggg=t.j;%g(x)按列排列GG=vertcat(GG,gg);%g垂直连接GendG=G;%转置得到矩阵Gfori=l:m%将数据y作为G的最后一列(n+1列)G(i,n+l)=yy；endG;fork=l:n%形成矩阵Q(k)ifG(k,k)O;sgn=l;else

8、ifG(kzk)=O;sgn=O;elsesgn=-l;endsgm=-sgn*sqrt(sum(G(k:m/k).A2);w=zeros(l,n);w(k)=G(k,k卜Sgm;forj=k+l:mw(j)=G(j,k);endbt=sgm*w(k);G(k,k)=sgm;%变换Gk-I到Gkforj=k+l:n+lt=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+l)G(nfn);%解三角方程求系数Afori=n-l:-l:lA(i)=(G(i,n+l)-sum(G(iJ+lrn).*A(i

9、+l:n)/G(iJ);endb=-A；c=A(2)(2*b);a=ep(A(l)+b*(c2);G(n+l:m/n+l)=exp(sum(G(n+l:mzn+l).A2);e=sum(G(n+l:mzn+l).A2);%计算误差efprintf(n指数函数的系数是：a=%f,b=%f,c=%f,a,b,c);%输出系数及误差efprintf(,n使用指数函数拟合的误差是:%f,e);t=0:0.05:24;YY=a.*exp(-b.*(t-c).2);figure(2)PlOt(X,y-k*力YYk,);%绘制拟合指数函数图形Xlabel(时亥J);%标注坐标轴含义ylabe(平均气温，);

10、titled指数函数的最小二乘曲线,);grid;(4)结果与分析a、二次函数：一天的平均气温为：21.20002次函数的系数：8.30632.6064-0.0938使用2次函数拟合的误差是：280.339547二次函数的最小二乘曲线如下图所示:302520b、三次函数：一天的平均气温为：21.20003次函数的系数：13.3880-0.22730.2075-0.0084使用3次函数拟合的误差是：131.061822三次函数的最小二乘曲线如下图所示:3；欠函数的最小二乘曲线15时刻2025302520明r口H-C、四次函数：一天的平均气温为：21.20004次函数的系数：16.7939-3.7

11、0500.8909-0.05320.0009使用4次函数拟合的误差是：59.04118四次函数的最小二乘曲线如下图所示:2；，d、指数函数：一天的平均气温为：21.2000指数函数的系数是：a=26.160286,b=0.004442,c=14.081900使用指数函数拟合的误差是：57.034644指数函数的最小二乘曲线如下图所示:通过上述几种拟合可以发现，多项式的次数越高，计算拟合的效果越好，误差越小，说明结果越准确；同时，指数多项式拟合的次数虽然不高，但误差最小，说明结果最准确。4.设计算法，求出非线性方程6/-45/+20=0的所有根,并使误差不超过10”。（1）算法思想首先，研究函数

12、的形态，确定根的范围；通过剖分区间的方法确定根的位置，然后利用二分法的基本原理进行求解，找到满足精度要求的解。二分法是产生一串区间，使新区间P是旧区间的一个子区间，其长度是/的一半，且有一个端点是/的一个端点。由区间确定区间/一的方法是计算区间的中点XE=I(Xw)若f(/(x(k+2)0thenstop3. IfL:I2then输出X-作为根；stop4. IfI/1then输出XT作为根；stop5. 2(O)+)=*6.Ifxm-xx(叫then输出X作为根；StoP7. f(x)=f8. IfI/IAhen输出X作为根；9. Iffjthen9.1 X=X%=clse9.2X=XC)J

13、nJIio.g0to5(3)MatIab源程序x=-100:100;y=6*(x.A5)-45*(x，2)+20;%非线性方程组的表达式g=；fori=-100:l:100%确定根所在的区间k=i+l;if(y(x=i).*y(=k)=10(-4)ifsubs(fzxzx)*subs(f,(x+xl)2)eps0=(0+l)2;elsexl=(x+xl)2;endendroot=xO%输出方程的根end(4)结果与分析该非线性方程组有三个实根，分别为1.8708,0.6812,-0.6545,且满足误差要求。5.编写程序实现大规模方程组的列主元高斯消去法程序，并对所附的方程组进行求解。针对本专

14、业中所碰到的实际问题，提炼一个使用方程组进行求解的例子，并对求解过程进行分析、求解。（1）算法思想高斯消去法是利用现行方程组初等变换中的一种变换，即用一个不为零的数乘一个方程后加只另一个方程，使方程组变成同解的上三角方程组，然后再自下而上对上三角方程组求解。列主元消去法是当高斯消元到第k步时，从A列的”以下（包括11O的各元素中选出绝对值最大的，然后通过行交换将其交换到On的位置上。交换系数矩阵中的两行（包括常数项），只相当于两个方程的位置交换了，因此，列选主元不影响求解的结果。程序的核心就是高斯列主元消去法。根据教材提供的算法，编写列主元消去法的子函数与适应于超大规模超出系统内存的方程组的改

15、编程序。同时，在GaUSS消去过程中，适当交换方程的顺序对保证消去过程能顺利进行及计算解的精确度都是有必要的，交换方程的原则是使寸=A/nL）中，绝对值最大的一个换到（k,k）位置而成为第k步消去的主元，这就是列主元Gauss消去法。（2）算法结构1、数据文件的文件名为：文件名+dat2、数据文件中的数据为二进制记录结构，分为以下四个部分：（1）文件头部分，其结构：typedefstruct（longintid;longintver;longintn;其中：id:为该数据文件的标识，值为OXFlElDlA0,即为：十六进制的FlEn)IAO说明系数矩阵为非压缩格式稀疏矩阵 系数矩阵为非压缩格式

16、带状对角阵 系数矩阵为压缩格式稀疏矩阵 系数矩阵为压缩格式带状对角阵ver:为数据文件的版本号，值为16进制数据，版本号0x1010x1020x2010x202n:表示方程的阶数(2)文件头2:此部分说明为条状矩阵的上下带宽,结构：typedefstruct(longintq;/为上带宽longintp;/为下带宽(3)系数矩阵a.如存贮格式非为压缩方式，则按行方式存贮系数矩阵中的每一个元素，个数为n*n,类型为float型；b.如果存贮格式是压缩方式,则按行方式存贮,每行中只存放上下带宽内的非零元素,即，每行中存贮的最多元素为p+q+1个。(4)右端系数按顺序存贮右端系数的每个元素，个数为n

17、个，类型为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(l);ver=dec2hex(a(2);%这两句是进行进制转换，读取id与Ver1. A的阶数2.For-l,2,3,n-l2.1找满足I%a=maxa脑的下标，k2.2For/=1,2,112.2.1。ffjHt)2.3=SRLRa业kA2.4For=fc+l+2

18、,n2.4.12.4.2For=k+l+2,112.4.2.1atakakaj2.4.3卜厂k=BtnPMnn11fok=n-l,n-2,1=*+Matlab源程序clear;%清除工作空间变量clc;%清除命令窗I命令%读取系数矩阵f,p=uigetfileC*.dat;选择数据文件。读取数据文件num=5;%输入系数矩阵文件头的个数name=strcat(pj);file=fopen(namer,);head=fread(file,numuint,);%读取二进制头文件id=dec2hex(head);读取标识符fprintf。文件标识符为少idver=dec2hex(head(2);%读

19、取版本，fprint用文件版本号为少vern=head(3);%读取阶数fprint矩阵A的阶数,);nq=head(4);%上带宽fprint年矩阵A的上带宽,);qp=head(5);%卜带宽fprin用矩阵A的下带宽,);Pdist=4*num;fseek(filezdistbof,);%把句柄值转向第六个兀素开头处Azcount=fread(filejnffloat,);%读取二进制文件,获取系数矩阵fclose(file);%关闭二进制头文件%对非压缩带状矩阵进行求解ifver=102a=zeros(nzn);fori=ln,forj=l:n,a(ij)=A(i-l)*n+j);%求

20、系数矩阵a(ij)endendb=zeros(nzl);fori=ln,b(i)=A(n*n+i);endfork=ln-lz%列主元高斯消去法m=k;fori=k+l:nz%寻找主元ifabs(a(mzk)abs(a(izk)m=i;endendifa(m,k)=0%遇到条件终止disp(错误！,)returnendforj=l:n,%交换元素位置得t=a(k,j);a(kj)=a(m,j);a(mj)=t;t=b(k);b(k)=b(m);b(m)=t;endfori=k+lm,%计算l(i,k)并将其放到a(izk),l,a(i)=a(ijk)a(kzk);forj=k+l:na(ij)

21、=a(ij)-a(izk)*a(k,j);endb(i)=b(i&a(i,k)*b(k);endendx=zeros(n,l);%回代过程x(n)=b(n)a(nzn);fork=n-l:-l:l.x(k)=(b(k)-sum(a(k,k+ln)*x(k+lrn)a(k,k);endend%对压缩带状矩阵进行求解讦ver=202,%高斯消去法m=p+q+l;a=zeros(nzm);fori=l:l:nforj=l:l:ma(ij)=A(i-l)*m+j);%求a(i,j)endendb=zeros(n,l);b(i)=A(n*m+i);%求b(i)endfork=l:l:(n-l)%开始消去

22、过程ifa(kz(p+l)=Odisp(错误！,)；break;endstl=n;if(k+p)nstl=k+p;endfori=(k+l):l:stla(iz(k+p-i+l)=a(i,(k+p-i+l)a(kz(p+l);forj=(k+l):l:(k+q)a(ij+p-i+l)=a(ij+p-i+l)-a(bk+p-i+l)*a(kzj+p-k+l);endb(i)=b(i)-a(i,k+p-i+l)*b(k);endend=zeros(n7l);%回代过程x(n)=b(n)a(nzp+l);sum=0;fork=(n-l):-l:lsum=b(k);st2=n;if(k+q)nst2=

23、k+q;endforj=(k+l):l:st2sum=sum-a(kj+p-k+l)*(j);endx(k)=suma(k,p+l);sum=0;endenddisp(,方程组的的解为：，)输出方程组的解disp()(3)结果与分析方程组一：文件标识符为id=Fieidiao文件版本号为Ver=102矩阵A的阶数n=15矩阵A的上带宽q=6矩阵A的下带宽p=6方程组的的解为：1.00001.00001.OOOO1.00001.00001.00001.00001.00001.00001.00001.00001.00001.00001.0000方程组二：文件标识符为id=Fieidiao文件版本号

24、为ver=102矩阵A的阶数n=2050矩阵A的上带宽q=6矩阵A的下带宽p=5方程组的的解为：1.96001.96001.96001.96001.96001.96001.96001.9600方程组三：文件标识符为id=Fieidiao文件版本号为ver=202矩阵A的阶数n=43512矩阵A的上带宽q=4矩阵A的下带宽p=3方程组的的解为：3.14133.14133.14133.14133.14133.14133.14133.14133.14133.1413实际应用:以求得三阶系统的第二阶固有频率，现通过特征方程来求解其主振型，其中A=7227;113011;7227,b=0O0,根据振动理

25、论归一化理论，取x(2)=l,计算出x(l)、X(2)o即可得到第二阶的主振型。利用GaUSS法求解可得x=T0U收获与体会首先，非常感谢老师百忙之中给我们讲授计算方法这门课程，使我对数值计算有了一个全新的认识，在课堂的学习中，我对数值计算方法有了一个基本的了解，但是这门课程要经过上机练习才能很好的掌握。根据老师给定的题目，我运用MatIab语言对题目进行了编程，并对计算结果进行了详细的分析。由于这是我第一次深入接触Matlab语言，在编程的过程中也遇到了不少困难，于是就去图书馆找资料，或者从网上查询每一条语句的用法，一句一句的编写程序，最终编写完了所有的程序，我的自信一下子提高了，享受到了劳

26、动成果的滋味。这次编程实践，是我学完理论课程之后对自己该方面的能力的一次很好的检验，从开始的算法思路，到运行调试，以及另人兴奋的可用程序，都是一个很好的学习和锻炼的过程。使我巩固了原有的理论知识，培养了我灵活运用和组合集成所学过知识及技能来分析解决实际问题的能力。使我体会到自身知识和能力能在实际中的应用和发挥。不但可以激发创新意识，还可以开发创造能力、培养沟通能力。报告中存在的不当之处，还请老师批评指正。计算方法上机报告姓名：学号：班级：上课班级：说明:本次上机实验使用的编程语言是Matlab语言，编译环境为MATLAB7.11.0,运行平台为WindOWS7。1.对以下和式计算：-（8nl8

27、n+4-8n+5-8n+6）-、，要求：若只需保留11个有效数字，该如何进行计算；若要保留30个有效数字，则又将如何进行计算;（1）算法思想1、根据精度要求估计所加的项数，可以使用后验误差估计，通项为：=1/421114l67rl8n+l8n+4-8n+5-811+616n8n+lxcf.2、为了保证计算结果的准确性，写程序时，从后向前计算；3、使用MatIab时，可以使用以下函数控制位数：digits（位数）或vpa（变量，精度为数）（2）算法结构1. 5=0；=TM8n+l-8n+4-8n+5-8n+6）l2. for=0/2Jift10end;3. for11-j-lj-2,s-s+t;

28、(3)Matlab源程序dear;%清除工作空间变量clc;%清除命令窗口命令m=inputC请输入有效数字的位数m=);%输入有效数字的位数s=0;forn=0:50t=(l16n)*(4(8*n+l)-2(8*n+4)-l(8*n+5)-l(8*n+6);ift=10(-m)%判断通项与精度的关系break;endend;fprint需要将n值加到n=%dN,nL);%需要将n值加到的数值fori=n-l:-l:0t=l16i)*(4(8*i+l)-2(8*i+4)-l(8*i+5)-l(8*i+6);s=s+t;%求和运算endS=VPa(S,m)%控制s的精度(4)结果与分析当保留11

29、位有效数字时，需要将n值加到n=7,s=3.1415926536；当保留30位有效数字时，需要将n值加到n=22,s=3.14159265358979323846264338328o通过上面的实验结果可以看出，通过从后往前计算，这种算法很好的保证了计算结果要求保留的准确数字位数的要求。2.某通信公司在一次施工中，需要在水面宽度为20米的河沟底部沿直线走向铺设一条沟底光缆。在铺设光缆之前需要对沟底的地形进行初步探测，从而估计所需光缆的长度，为工程预算提供依据。已探测到一组等分点位置的深度数据（单位：米）如下表所示：分点O123456深度9.018.967.967.978.029.0510.13分

30、点78910111213深度11.1812.2613.2813.3212.6111.2910.22分点14151617181920深度9.157.907.958.869.8110.8010.93请用合适的曲线拟合所测数据点；预测所需光缆长度的近似值，作出铺设河底光缆的曲线图；（1）算法思想如果使用多项式差值，则由于龙格现象，误差较大，因此，用相对较少的插值数据点作插值，可以避免大的误差，但是如果又希望将所得数据点都用上，且所用数据点越多越好，可以采用分段插值方式，即用分段多项式代替单个多项式作插值。分段多项式是由一些在相互连接的区间上的不同多项式连接而成的一条连续曲线，其中三次样条插值方法是一

31、种具有较好“光滑性”的分段插值方法。在本题中，假设所铺设的光缆足够柔软，在铺设过程中光缆触地走势光滑,紧贴地面，并且忽略水流对光缆的冲击。海底光缆线的长度预测模型如下所示，光缆从A点铺至B点，在某点处的深度为10海底光缆线的长度预测模型计算光缆长度时，用如下公式:L 2。/(x)ds0=f2(x)Vl(x)2dxJo=r%w+)2d =0 kV(x)2+(y)2(2)算法结构1.FOri0,12,Ll尸2.Fork-1,22.1Fori=ntn-lftk2,Li(MLMI)/(一&_女)=3.81-刈=力14.FOri=12/4.1储+1一万户力*14.2力Hd/=C|；1一。产。,；2=。4.36M114=dr5获取M的矩阵元素个数，存入In8. ForA-2,3,，m&lI=，8.2*,*尸乙8.3九/丫k-1=乙9.丫mm=Mm10.ForA加一1,】一2,110.1(.c卜Mk41)/UlnMkv.*H,3一七人34、*+1k11.获取