计算方法上机作业重点讲义资料

计算方法上机报告计算方法上机报告姓 名：学 号：班 级：上课班级：计算方法上机报告说明： 本次上机实验使用的编程语言是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; %清除工作空间变量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）结果与分析 当保留11位有效数字时，需要将n值加到n=7， s =3.； 当保留30位有效数字时，需要将n值加到n=22, s =3.。 通过上面的实验结果可以看出，通过从后往前计算，这种算法很好的保证了计算结果要求保留的准确数字位数的要求。2. 某通信公司在一次施工中，需要在水面宽度为20米的河沟底部沿直线走向铺设一条沟底光缆。在铺设光缆之前需要对沟底的地形进行初步探测，从而估计所需光缆的长度，为工程预算提供依据。已探测到一组等分点位置的深度数据(单位:米)如下表所示：分点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）算法思想 如果使用多项式差值，则由于龙格现象，误差较大，因此，用相对较少的插值数据点作插值，可以避免大的误差，但是如果又希望将所得数据点都用上，且所用数据点越多越好，可以采用分段插值方式，即用分段多项式代替单个多项式作插值。分段多项式是由一些在相互连接的区间上的不同多项式连接而成的一条连续曲线，其中三次样条插值方法是一种具有较好“光滑性”的分段插值方法。在本题中，假设所铺设的光缆足够柔软，在铺设过程中光缆触地走势光滑，紧贴地面，并且忽略水流对光缆的冲击。海底光缆线的长度预测模型如下所示，光缆从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 ；break 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)=(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)=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(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及误差edisp(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)按列排列 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); 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).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 -0.0938使用2次函数拟合的误差是：280.二次函数的最小二乘曲线如下图所示：b、三次函数：一天的平均气温为： 21.20003次函数的系数: 13.3880 -0.2273 0.2075 -0.0084使用3次函数拟合的误差是： 131. 三次函数的最小二乘曲线如下图所示：c、四次函数：一天的平均气温为： 21.20004次函数的系数: 16.7939 -3.7050 0.8909 -0.0532 0.0009使用4次函数拟合的误差是：59.04118四次函数的最小二乘曲线如下图所示：d、指数函数：一天的平均气温为： 21.2000指数函数的系数是： a=26.，b=0.，c=14.使用指数函数拟合的误差是： 57.指数函数的最小二乘曲线如下图所示：通过上述几种拟合可以发现，多项式的次数越高，计算拟合的效果越好，误差越小，说明结果越准确；同时，指数多项式拟合的次数虽然不高，但误差最小，说明结果最准确。4.设计算法，求出非线性方程的所有根，并使误差不超过。（1）算法思想 首先，研究函数的形态，确定根的范围；通过剖分区间的方法确定根的位置，然后利用二分法的基本原理进行求解，找到满足精度要求的解。二分法是产生一串区间，使新区间是旧区间的一个子区间，其长度是的一半，且有一个端点是的一个端点。由区间确定区间的方法是计算区间的中点若，则取，否则取，重复这一过程即可。显然，每次迭代使区间长度减小一半，故二分法总是收敛的。（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).*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）算法思想 高斯消去法是利用现行方程组初等变换中的一种变换，即用一个不为零的数乘一个方程后加只另一个方程，使方程组变成同解的上三角方程组，然后再自下而上对上三角方程组求解。列主元消去法是当高斯消元到第步时，从列的以下（包括）的各元素中选出绝对值最大的，然后通过行交换将其交换到的位置上。交换系数矩阵中的两行（包括常数项），只相当于两个方程的位置交换了，因此，列选主元不影响求解的结果。 程序的核心就是高斯列主元消去法。根据教材提供的算法，编写列主元消去法的子函数与适应于超大规模超出系统内存的方程组的改编程序。同时，在Gauss消去过程中，适当交换方程的顺序对保证消去过程能顺利进行及计算解的精确度都是有必要的，交换方程的原则是使中，绝对值最大的一个换到（k,k)位置而成为第k步消去的主元，这就是列主元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:表示方程的阶数 (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=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,选择数据文件); %读取数据文件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); %关闭二进制头文件% 对非压缩带状矩阵进行求解if ver=102, a=zeros(n,n); for i=1:n, for j=1:n, a(i,j)=A(i-1)*n+j); %求系数矩阵a(i,j) end end b=zeros(n,1); for i=1:n, b(i)=A(n*n+i); end 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,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=n-1:-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(i-1)*m+j); %求a(i,j) end end b=zeros(n,1); for i=1:1:n b(i)=A(n*m+i); %求b(i) end for k=1:1:(n-1) %开始消去过程 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+p-i+1)=a(i,(k+p-i+1)/a(k,(p+1); for j=(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); end b(i)=b(i)-a(i,k+p-i+1)*b(k); end end x=zeros(n,1); %回代过程 x(n)=b(n)/a(n,p+1); sum=0; for 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=sum-a(k,j+p-k+1)*x(j); end x(k)=sum/a(k,p+1); sum=0; endenddisp(方程组的的解为：) %输出方程组的解disp(x)（3） 结果与分析方程组一：文件标识符为 id =F1E1D1A0文件版本号为 ver =102矩阵A的阶数 n = 15矩阵A的上带宽 q =6矩阵A的下带宽 p =6方程组的的解为： 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方程