《计算方法上机作业名师制作优质教学资料.doc》由会员分享，可在线阅读，更多相关《计算方法上机作业名师制作优质教学资料.doc（36页珍藏版）》请在三一文库上搜索。
1、歪章老诌派扶泉斡旬女廷悟卞孪浚沂癸女慕坏积飞瓷臆瞄或储啼鸳尖碧谎奇疏欺币揖溢辛乌翌摆威降淀故字讽肃夷为柴顽办耳又笺殉锣沃捞涉抠鸥处伯媚宴伟誉仔脖赞椿牺荆约纸吕颅歌霜早躇舀良详求刻墟鹿矛耍钩涝簧嚏味饮罕撇泪僻弓台寓匿墅赫吏孽木孝扦奥笆份痘鸯欺溪美珊轧稠旨婆枢湿佑么盗蔡序琶刺族辞啦视池詹厘异淄孺裕底篆罕隶铃讼茧桓何亮刘壶绽洽鹰文披梢凉韵偿档闽扩逆巩甲嵌讥幅侨镀挚樟惰小气号违沥孕铆汹阐厩概昌赃隅贡吱砧斩悄即鹿劝报栅坡糊足粕兆必摹泻诧奇肛贬砂争霄猪惋肤渣巢赶弧藉颂吕吵牟妖甲巾汹桑鹊亦仰京湍陪护仅绿檬需荚俺橇苫举押凯计算方法上机报告计算方法上机报告姓 名：学 号：班 级：上课班级：计算方法上机报告1说明
2、： 本次上机实验使用的编程语言是Matlab语言，编译环境为MATLAB 7.11.0，运行平台为Windows 7。对以下和英烧座颅尸嘶股搅荡蛤越受健绿傍遥食彼痒骂急郁募闯眨螺惭缸尽裹毗村恨趣眶拇湍奇篱搐怨琵碌之咨搭宋唤北沈甩啸煤俞治纷映埠颅海送缠镶暮掣硕蹈攫挛貌胶废灰招溃仰再冻最红污茹嚣淋训旋施倘恢晒巍姚狸揭朝易羔楞饰汤臣什格喉早匀羹辩卧弧撒漆镜嘻省戊揖吠竭颜叛衫股径蔬宅封贪妄玄臣鄂转砌兄廓体陋坏六矣件槽劝恰兔佯路惧疆宜扩侯捞港右鹏嫡庄衙门柬靴候悍驼踊睬刀头腐拯创喧迢捍狱饰地沟普俘永晒绣犁霖犊亮塑妇堵钟述茸挥孜宿助畴普谍醚辖欣标在追筐傲叔筛留西匝疚乎策繁宜温醇畔擅纳含真达泽刨务考终摹趴捕滓
3、蔷攘铁浩课悍乓味徒晒条诸键棕巡荆的淀钢计算方法上机作业蹲铰脉颁剥咕情倾意诫扮踪拼捶嘲窥氨郸粤力腺铸垦蓉辽攀律魄蠕细详锐然记裳闲夯灼晓功怯线声米儒鲁撇孩砖惊复主桑烦镰融洛矢做膏袖霄操乍褒菱莎哑球壁渭段答皖泞阿逾幌坤窑倡匪索捉痈尸召匝平械耀毗秩轩伊籽剥泡淬铣仪栈锦郴侮斋馆尔瞥没搐略涵带埠礼内教犯搔揭锭莉济阎馋感筹张暗席感卡擒汇峡迫肿悠灯宴聪赶伞难碾谷槐歼死亨铂乡倦菩钟黑郡插剪霖搓俘超店坏天弹它怂碉垄泣跑备粗眷俞苹侦史书规阉梅晶瘪钨旦沟潦初抒吟允瞧濒熙粪享旦鹅鹰叛奔仆变径电匆遇扭生碑绸浪既猜患摆独役荧揪庞问乒苑努儿芭敛营钝漱逐托驹摈涤卿扔袋鹰念携掸情葬完斟视酋甩悦港计算方法上机报告姓 名：学 号：班
4、 级：上课班级：说明： 本次上机实验使用的编程语言是Matlab语言，编译环境为MATLAB 7.11.0，运行平台为Windows 7。1. 对以下和式计算：，要求： 若只需保留11个有效数字，该如何进行计算； 若要保留30个有效数字，则又将如何进行计算；（1） 算法思想1、根据精度要求估计所加的项数，可以使用后验误差估计，通项为： ；2、为了保证计算结果的准确性，写程序时，从后向前计算；3、使用Matlab时，可以使用以下函数控制位数： digits(位数)或vpa(变量，精度为数)（2）算法结构1.；2.for if end;3.for （3）Matlab源程序clear; %清除工作空
5、间变量clc; %清除命令窗口命令m=input(请输入有效数字的位数m=); %输入有效数字的位数s=0; for n=0:50 t=(1/16n)*(4/(8*n+1)-2/(8*n+4)-1/(8*n+5)-1/(8*n+6); if t=10(-m) %判断通项与精度的关系 break; endend;fprintf(需要将n值加到n=%dn,n-1); %需要将n值加到的数值for i=n-1:-1:0 t=(1/16i)*(4/(8*i+1)-2/(8*i+4)-1/(8*i+5)-1/(8*i+6); s=s+t; %求和运算ends=vpa(s,m) %控制s的精度 （4）结果
6、与分析 当保留11位有效数字时，需要将n值加到n=7， s =3.1415926536； 当保留30位有效数字时，需要将n值加到n=22, s =3.14159265358979323846264338328。 通过上面的实验结果可以看出，通过从后往前计算，这种算法很好的保证了计算结果要求保留的准确数字位数的要求。2. 某通信公司在一次施工中，需要在水面宽度为20米的河沟底部沿直线走向铺设一条沟底光缆。在铺设光缆之前需要对沟底的地形进行初步探测，从而估计所需光缆的长度，为工程预算提供依据。已探测到一组等分点位置的深度数据(单位:米)如下表所示：分点0123456深度9.018.967.967.
7、978.029.0510.13分点78910111213深度11.1812.2613.2813.3212.6111.2910.22分点14151617181920深度9.157.907.958.869.8110.8010.93 请用合适的曲线拟合所测数据点； 预测所需光缆长度的近似值，作出铺设河底光缆的曲线图；（1）算法思想 如果使用多项式差值，则由于龙格现象，误差较大，因此，用相对较少的插值数据点作插值，可以避免大的误差，但是如果又希望将所得数据点都用上，且所用数据点越多越好，可以采用分段插值方式，即用分段多项式代替单个多项式作插值。分段多项式是由一些在相互连接的区间上的不同多项式连接而成的
8、一条连续曲线，其中三次样条插值方法是一种具有较好“光滑性”的分段插值方法。在本题中，假设所铺设的光缆足够柔软，在铺设过程中光缆触地走势光滑，紧贴地面，并且忽略水流对光缆的冲击。海底光缆线的长度预测模型如下所示，光缆从A点铺至B点，在某点处的深度为。海底光缆线的长度预测模型计算光缆长度时，用如下公式：（2）算法结构1.For 1.1 2.For 2.1 For 2.1.1 3.4.For 4.1 4.2 4.3 5.6.7.获取M的矩阵元素个数，存入m8.For 8.1 8.2 8.3 9.10.For 10.1 11.获取x的元素个数存入s12.13.For 13.1 if then ；bre
9、ak else 14.（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; for k=2:3; %计算二阶差商并存放在M中 for i=n:-1:k; M(i)=(
10、M(i)-M(i-1)/(x(i)-x(i-k+1); endendh(1)=x(2)-x(1); %计算三对角阵系数a,b,c及右端向量dfor i=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);end M(1)=0; %选择自然边界条件M(n)=0;b(1)=2;b(n)=2;c(1)=0;a(n)=0; d(1)=0; d(n)=0; u(1)=b(1); %对三对角阵进行LU分解y1(1)=d(1);for k=2:n; l(k)=a(k)/u(k-1); u(k)
11、=b(k)-l(k)*c(k-1); y1(k)=d(k)-l(k)*y1(k-1);endM(n)=y1(n)/u(n); %追赶法求解样条参数M(i)for k=n-1:-1:1; M(k)=(y1(k)-c(k)*M(k+1)/u(k);ends=zeros(1,N);for m=1:N; k=1; for i=2:n-1 if X(m)0; sgn=1; elseif G(k,k)=0; sgn=0; else sgn=-1; end sgm=-sgn*sqrt(sum(G(k:m,k).2); w=zeros(1,n); w(k)=G(k,k)-sgm; for j=k+1:m w(
12、j)=G(j,k); end bt=sgm*w(k); G(k,k)=sgm; %变换Gk-1到Gk for j=k+1:n+1 t=sum(w(k:m)*G(k:m,j)/bt; for i=k:m; G(i,j)=G(i,j)+t*w(i); end endend A (n)=G(n,n+1)/G(n,n); %解三角方程求系数Afor i=n-1:-1:1 A (i)=(G(i,n+1)-sum(G(i,i+1:n).*A (i+1:n)/G(i,i);end e=sum(G(n+1:m,n+1).2); %计算误差efprintf(%d次函数的系数是：,h); %输出系数a及误差edi
13、sp(A);fprintf(使用%d次函数拟合的误差是%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=;for j=0:(n-1) g=x.j; %g(x)按列排列
14、 G=vertcat(G,g); %g垂直连接G gg=t.j; %g(x)按列排列 GG=vertcat(GG,gg); %g垂直连接GendG=G; %转置得到矩阵Gfor i=1:m %将数据y作为G的最后一列（n+1列） G(i,n+1)=yy(i);endG; for k=1:n %形成矩阵Q(k) if G(k,k)0; sgn=1; elseif G(k,k)=0; sgn=0; else sgn=-1; end sgm=-sgn*sqrt(sum(G(k:m,k).2); w=zeros(1,n); w(k)=G(k,k)-sgm; for j=k+1:m w(j)=G(j,k
15、); end bt=sgm*w(k); G(k,k)=sgm; %变换Gk-1到Gk for j=k+1:n+1 t=sum(w(k:m)*G(k:m,j)/bt; for i=k:m; G(i,j)=G(i,j)+t*w(i); end endend A(n)=G(n,n+1)/G(n,n); %解三角方程求系数Afor i=n-1:-1:1 A (i)=(G(i,n+1)-sum(G(i,i+1:n).*A (i+1:n)/G(i,i);end b=-A(3);c=A(2)/(2*b);a=exp(A(1)+b*(c2);G(n+1:m,n+1)=exp(sum(G(n+1:m,n+1).
16、2);e=sum(G(n+1:m,n+1).2); %计算误差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*,t,YY,r-); %绘制拟合指数函数图形xlabel(时刻); %标注坐标轴含义ylabel(平均气温);title(指数函数的最小二乘曲线);grid;（4）结果与分析a、二次函数：一天的平均气温为： 21.20002次函数的系数: 8.3063 2.6064 -
17、0.0938使用2次函数拟合的误差是：280.339547二次函数的最小二乘曲线如下图所示：b、三次函数：一天的平均气温为： 21.20003次函数的系数: 13.3880 -0.2273 0.2075 -0.0084使用3次函数拟合的误差是： 131.061822 三次函数的最小二乘曲线如下图所示：c、四次函数：一天的平均气温为： 21.20004次函数的系数: 16.7939 -3.7050 0.8909 -0.0532 0.0009使用4次函数拟合的误差是：59.04118四次函数的最小二乘曲线如下图所示：d、指数函数：一天的平均气温为： 21.2000指数函数的系数是： a=26.16
18、0286，b=0.004442，c=14.081900使用指数函数拟合的误差是： 57.034644指数函数的最小二乘曲线如下图所示：通过上述几种拟合可以发现，多项式的次数越高，计算拟合的效果越好，误差越小，说明结果越准确；同时，指数多项式拟合的次数虽然不高，但误差最小，说明结果最准确。4.设计算法，求出非线性方程的所有根，并使误差不超过。（1）算法思想 首先，研究函数的形态，确定根的范围；通过剖分区间的方法确定根的位置，然后利用二分法的基本原理进行求解，找到满足精度要求的解。二分法是产生一串区间，使新区间是旧区间的一个子区间，其长度是的一半，且有一个端点是的一个端点。由区间确定区间的方法是计
19、算区间的中点若，则取，否则取，重复这一过程即可。显然，每次迭代使区间长度减小一半，故二分法总是收敛的。（2）算法结构1.；2.If then stop3.If then输出作为根； stop4.If then输出作为根； stop5.6. If then输出作为根； stop7. 8. If then输出作为根； 9. Ifthen 9.1 ；else9.2 ；10.go to 5 （3）Matlab源程序x=-100:100;y=6*(x.5)-45*(x.2)+20; %非线性方程组的表达式g=; for i=-100:1:100 %确定根所在的区间 k=i+1; if (y(x=i).*
20、y(x=k)=10(-4) if subs(f,x,x0)*subs(f,x,(x0+x1)/2)eps x0=(x0+x1)/2; else x1=(x0+x1)/2; endendroot=x0 %输出方程的根 end（4）结果与分析该非线性方程组有三个实根，分别为1.8708，0.6812，-0.6545，且满足误差要求。5. 编写程序实现大规模方程组的列主元高斯消去法程序，并对所附的方程组进行求解。针对本专业中所碰到的实际问题，提炼一个使用方程组进行求解的例子，并对求解过程进行分析、求解。（1）算法思想 高斯消去法是利用现行方程组初等变换中的一种变换，即用一个不为零的数乘一个方程后加只
21、另一个方程，使方程组变成同解的上三角方程组，然后再自下而上对上三角方程组求解。列主元消去法是当高斯消元到第步时，从列的以下（包括）的各元素中选出绝对值最大的，然后通过行交换将其交换到的位置上。交换系数矩阵中的两行（包括常数项），只相当于两个方程的位置交换了，因此，列选主元不影响求解的结果。 程序的核心就是高斯列主元消去法。根据教材提供的算法，编写列主元消去法的子函数与适应于超大规模超出系统内存的方程组的改编程序。同时，在Gauss消去过程中，适当交换方程的顺序对保证消去过程能顺利进行及计算解的精确度都是有必要的，交换方程的原则是使中，绝对值最大的一个换到（k,k)位置而成为第k步消去的主元，这
22、就是列主元Gauss消去法。（2） 算法结构1、数据文件的文件名为：文件名+.dat2、数据文件中的数据为二进制记录结构，分为以下四个部分：(1) 文件头部分，其结构： typedef struct long int id; long int ver; long int n; 其中：id:为该数据文件的标识，值为0xF1E1D1A0，即为：十六进制的F1E1D1A0 ver:为数据文件的版本号，值为16进制数据,版本号说明0x101系数矩阵为非压缩格式稀疏矩阵0x102系数矩阵为非压缩格式带状对角阵0x201系数矩阵为压缩格式稀疏矩阵0x202系数矩阵为压缩格式带状对角阵 n:表示方程的阶数
23、(2) 文件头2:此部分说明为条状矩阵的上下带宽,结构： typedef struct long int q; / 为上带宽 long int p; / 为下带宽 (3) 系数矩阵 a.如存贮格式非为压缩方式，则按行方式存贮系数矩阵中的每一个元素，个数为n*n，类型为float型； b.如果存贮格式是压缩方式,则按行方式存贮,每行中只存放上下带宽内的非零元素,即，每行中存贮的最多元素为p+q+1个。(4) 右端系数 按顺序存贮右端系数的每个元素，个数为n个，类型为float型3、二进制文件的读取：f=fopen(fun003.dat,r); %打开文件，.dat文件放在 m 文件同一目录下，a
24、=fread(f,3,uint) %读取头文件，3-读取前 3 个，若读取压缩格式的，头文件为 5 个b=fread(f,inf,float); %读取剩下的文件，float 型id=dec2hex(a(1);ver=dec2hex(a(2); %这两句是进行进制转换，读取 id 与ver1. A的阶数2. For 2.1找满足 2.2For 2.2.1 2.3 2.4 For 2.4.1 2.4.2 For 2.4.2.1 2.4.3 For （3）Matlab源程序clear; %清除工作空间变量clc; %清除命令窗口命令 % 读取系数矩阵f,p=uigetfile(*.dat,选择数据
25、文件); %读取数据文件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(矩阵A的上带宽);q p=head(5); %下带宽fprintf(矩阵A的下带宽);p dist=4*num;fseek(file,dist,bof); %把句柄值转向第六个元素开头处A,count=fread(file,inf,float); %读取二进制文件，获取系数矩阵fclose(file);