西安交通大学计算方法B上机报告(完整题目和程序).doc

西安交通大学计算方法B 上机编程报告学号：XXX姓名：XXX专业：工程热物理 班级：硕XXXX date2015/12/10注：本上机报告使用的程序语言均为Matlab语言，为本人独立完成！1. 计算以下和式：，要求：(1)若保留11个有效数字，给出计算结果，并评价计算的算法；(2)若要保留30个有效数字，则又将如何进行计算。a) 解题思想（1） 先根据精度要求估计所需累加的项数n，使用后验误差估计方法，条件为： (m为有效数字位数)。（2） 在该题中S的和式存在两个相近的数相减的问题，为了避免有效数字损失，在计算中改变了运算顺序，分别将正数和负数分别相加，然后再将其和相加。（3） 为了避免大数吃小数的问题，本题先计算出保留目标有效数字所需要的迭代次数，然后采用倒序相加的方法提高计算精度。b) 算法实现的结构 1. S1=s2=0； ;2. for n=0,1,2,iIf end 3. for n=i,i-1,i-2,0a1=4/(16n*(8*n+1);a2=2/(16n*(8*n+4);a3=1/(16n*(8*n+5);a4=1/(16n*(8*n+6);s1=a1+s1;s2=a4+a3+a2+s2;endS=s1-s2;c) 计算源程序clear; %清除工作空间变量clc; %清除命令窗口命令m=input(请输入有效数字的位数m=); %输入需要的有效数字位数S=0;s1=0;s2=0; %定义存储正数、负数和累加和的变量for n=0:1:1000t=(1/16n)*(4/(8*n+1)-(2/(8*n+4)-1/(8*n+5)-1/(8*n+6); if t=10(-m) %根据有效数字位数确定所需累加的n值 break; end k=n;end;for n=(k-1):-1:0 a1=4/(16n*(8*n+1); a2=2/(16n*(8*n+4); a3=1/(16n*(8*n+5); a4=1/(16n*(8*n+6); s1=a1+s1; %第一项倒序累加 s2=a4+a3+a2+s2; %后三项倒序累加endS=s1-s2; %正数累加值和负数累加值的和S=vpa(S,m) %控制S的精度d) 计算结果与评价当需保留11位有效数字时，需要将n值加到n=7， s=3.1415926536；当需保留30位有效数字时，需要将n值加到n=22, s=3.14159265358979311599796346854。由计算结果可以看出，采用从后往前进行计算的方式，避免了“大数吃小数”的问题，这种算法很好的保证了计算结果要求保留的准确数字有效位数的要求。另外，中有多个负数相加，按照绝对值递增的顺序求和，减少了舍入误差带来的影响。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)请用合适的曲线拟合所测数据点；(2)估算所需光缆长度的近似值，并作出铺设河底光缆的曲线图；a) 解题思想该题需要对所需光缆长度进行近似估计，通过对测量数据进行拟合并进行曲线积分即可得到河底光缆长度的近似值。由于数值点较多，如果使用多项式插值，会出现龙格现象导致误差较大，因此，用相对较少的插值数据点作插值，可以避免大的误差，但又希望将所得数据点都用上，且所用数据点越多越好，所以本题采用分段插值方式，即用分段多项式代替单个多项式作插值。分段多项式是由一些在相互连接的区间上的不同多项式连接而成的一条连续曲线，其中三次样条插值方法是一种具有较好“光滑性”的分段插值方法。在本题中，假设所铺设的光缆足够柔软，在铺设过程中光缆触地走势光滑，紧贴地面，并且忽略了水流对光缆的冲击，故本题使用三次样条插值进行计算精度是最高的。b) 算法实现的结构1 For 1.12 For 2.1 For 2.1.134 For 4.14.24.3567 获取M的矩阵元素个数，记为m8 For 8.18.28.3910 For 10.111 获取x的元素个数存入s1213 For 13.1 If then breakelse 14 15 计算光缆长度时，用如下公式： c) 计算源程序clear; %清除工作空间变量clc; %清除命令窗口命令x=0:1:20; %将分点位置以数组的形式存储于x中X=0:0.2:20; %将插值后要求的点存储在X中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)，见课本P110M=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);endM(1)=0; %补充自然边界条件M(n)=0;a(n)=0;b(1)=2;b(n)=2;c(1)=0;d(1)=0;d(n)=0;u(1)=b(1); %对三对角阵进行LU分解Y(1)=d(1);for k=2:n;l(k)=a(k)/u(k-1);u(k)=b(k)-l(k)*c(k-1);Y(k)=d(k)-l(k)*Y(k-1);endM(n)=Y(n)/u(n); %TSS算法求解样条参数M(i)，见课本P49for k=n-1:-1:1;M(k)=(Y(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)=x(i); k=i-1; break; else k=i; end end H=x(k+1)-x(k); %在各区间用三次样条插值函数计算待求X点处的s值 x1=x(k+1)-X(m); x2=X(m)-x(k);s(m)=(M(k)*(x13)/6+M(k+1)*(x23)/6+(y(k)-(M(k)*(H2)/6)*x1+(y(k+1)-(M(k+1)*(H2)/6)*x2)/H;end%计算所需光缆长度L=0;for i=2:NL=L+sqrt(X(i)-X(i-1)2+(s(i)-s(i-1)2); %利用第一类线积分求河底光缆的长度enddisp(所需光缆长度为L=);disp(L);figureplot(x,y,*,X,s,-) %绘制铺设河底光缆的曲线图xlabel(分点位置,fontsize,16); %标注坐标轴含义ylabel(测点深度/m,fontsize,16);title(铺设河底光缆的插值曲线图,fontsize,16);grid on;d) 计算结果与评价计算所得的运行结果：铺设海底光缆的插值曲线图如所示：上述插值结果表明，采用分段三次样条插值所得的拟合曲线能较准确地反映铺设光缆的走势，可以有效避免龙格现象。在本题的计算中采用自然三次样条函数作为边界条件，在解线性方程组时使用了TSS算法求解带状矩阵，其计算速度快，是一种求解线性方程组的有效手段，在估计河底光缆长度时使用第一类线积分，最终计算出所需光缆的长度为L=26.4844m。在计算编程过程中对三次样条插值的原理理解不透彻，导致编程逻辑思路不清晰，调试时间较长，体会到编程是一项仔细认真的工作，不能马虎！3. 假定某天的气温变化记录如下表所示，试用数据拟合的方法找出这一天的气温变化的规律;试计算这一天的平均气温,并试估计误差。时刻0123456789101112平均气温15141414141516182020232528时刻131415161718192021222324平均气温313431292725242220181716a) 解题思想该题要求拟合一天的气温变化规律，由于记录数据点的数目较多。需要采用很高次的多项式做数据近似，这会给计算带来困难，也会导致误差很大。用“插值样条函数”做数据近似，虽然有很好的数值性质，且计算量也不大，但存放参数Mi的量很大，且没有一个统一的数学公式来表示，也带来了一些不便。另一方面，在有的实际问题中，用插值方法并不合适，当数据点的数目很大时，要求通过所有数据点，可能会失去原数据所表示的规律，如果数据点是由测量而来的，必然带有误差，插值法要求准确通过这些不准确的数据点是不合适的。在这种情况下，不用插值标准而用选取使最小的标准，即采用最小二乘近似方法。本题采用“最小二乘法”找出这一天的气温变化的规律，考虑使用最小二乘的二次函数、三次函数、四次函数以及指数函数，计算相应的系数，估算误差，问题的难度是求解各种拟合函数的系数。由于利用法方程求解最小二乘系数时，方程的解不够稳定，所以本题利用正交化方法解最小二乘问题。b) 算法结构用正交化方法求数据的最小二乘近似问题，假定数据以用来生成了G，并将y作为其最后一列（第n+1列）存放，结果在a中，是误差是E2。(一) 使用高次多项式拟合1 将“时刻值”存入，数据点的个数存入N2 输入拟合多项式函数的最高项次数n-1，则拟合多项式函数为，。3 根据给定数据点确定G,并将y作为n+1列存入GFor For 4 For 4.1 形成矩阵4.1.14.1.24.1.3 For 4.1.3.14.1.44.2 变换到4.2.1For 4.2.24.2.3 For 4.2.3.15 解三角方程5.15.2 For 5.2.16 计算误差 (二) 使用指数函数拟合当为指数函数时，假定指数函数的形式为。只需将y值求对数，其主体部分不变，令，则可得多项式:。 c) 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;N=length(x); n=input(请输入所需拟合函数的最高项次数n=);plot(x,y,r*)grid;hold on;n=n+1;G=zeros(N,n+1); %定义零矩阵G(:,n+1)=y; %把y作为G的最后一列存放C=zeros(1,n); %建立C矩阵来存放E=0;F=0; %定义计算用的中间变量B=zeros(1,N); %建立B用来存放%根据已给数据生成矩阵Gfor j=1:n for i=1:N G(i,j)=x(1,i)(j-1); endend%形成矩阵Qk，见课本P123for k=1:n for i=k:N C(k)=G(i,k)2+C(k); end C(k)=-sign(G(k,k)*(C(k)0.5); %计算的值存入C中 w(k)=G(k,k)-C(k); %建立w来存放 for j=k+1:N w(j)=G(j,k); end B(k)=C(k)*w(k); %计算的值并存入B中 %变换矩阵Gk-1到Gk，参考课本P123 G(k,k)=C(k); for j=k+1:n+1 E=0; for i=k:N E=w(i)*G(i,j)+E; end t=E/B(k); for i=k:N G(i,j)=t*w(i)+G(i,j); end endend%求解三角方程：Rx=h1a(n)=G(n,n+1)/G(n,n);for i=n-1:(-1):1 for j=i+1:n F=G(i,j)*a(j)+F; end a(i)=(G(i,n+1)-F)/G(i,i); % a(i)存放各级系数 F=0;enddisp(各级系数分别为 ,num2str(a)%拟合后的回代求解过程p=zeros(1,N);for j=1:N for i=1:n p(j)=p(j)+a(i)*x(j)(i-1); endendplot(x,p,bx,x,p,-);xlabel(时刻/h,fontsize,16); %标注坐标轴含义ylabel(平均温度/,fontsize,16);title(拟合平均气温曲线图,fontsize,16);grid on;E2=0; %用E2来存放误差%计算拟合误差for i=n+1:N E2=G(i,n+1)2+E2;endE2=E20.5;disp(拟合所得误差为E2=,num2str(E2);Tm=0;for i=1:N Tm=Tm+p(i);endt=Tm/N;%计算平均温度disp(平均温度为t=,num2str(t),)d) 计算结果与分析（1） 采用二次函数拟合故可得函数形式为： 二次函数的最小二乘曲线如下图所示：（2） 采用三次函数拟合故可得函数形式为： 三次函数的最小二乘曲线如下图所示：（3） 采用四次函数拟合故可得函数形式为：四次函数的最小二乘曲线如下图所示：（4） 采用指数函数拟合故可得函数形式为： 指数函数的最小二乘曲线如下图所示：通过上述拟合结果可以发现，使用最小二乘法进行数据拟合时，多项式的次数越高，拟合曲线也更加接近数据点，计算拟合的效果越好，误差越小，结果越准确。同时，指数多项式拟合的次数虽然不高，但误差最小，是因为指数型二次函数是对数据点求对数进行拟合的，故结果最准确。在进行最小二乘法编程时发现，使用正交化方法求解最小二乘问题比使用法方程方法更稳定，结果更可靠。但编程是中间变量较多容易造成混乱，计算量也较大，但本题在编程过程中有一定的适用性，只需要对数据点进行处理就可以拟合出给定次数的多项式，为以后的数据拟合工作打下了基础。4. 设计算法，求出非线性方程的所有实根，并使误差不超过。a) 解题思想由题目所给的待求方程可知该方程为5次多项式，导数的求值比较容易，可采用牛顿法迭代求解。具体实现步骤为：首先，研究函数的形态，确定根的范围；通过剖分区间的方法确定根的位置，然后利用牛顿法的基本原理进行求解，找到满足精度要求的解。令，其迭代格式为要使牛顿迭代法收敛，则必须寻找根的合理区间a,b,使得在该区间内，即有根。在选定的区间内函数的一二阶导数均不变号。选定一个初值，使，则牛顿迭代求解过程完成。b) 算法结构1 For 1.11.2 If then1.2.1 输出“的值已足够小”1.2.2 ；stop else1.2.31.2.41.2.5 If then1.2.5.1 输出“已足够小”；stop else1.2.5.22 输出“次迭代后仍不收敛”，stopc) 计算源程序clear;clc;%根据函数曲线观察根的大致分布区间X=-8:0.1:8;for i=1:161Y(i)=6*X(i)5-45*X(i)2+20;endplot(X,Y,r-,LineWidth,2); %画出Y随X变化的曲线图%标注坐标轴含义及图像标题xlabel(X,fontsize,10);ylabel(Y,fontsize,10);title(多项式函数变化曲线图,fontsize,10);grid on;k=1;Y(162)=0;%使用二分法来确定不同根区间的迭代初值，并存于a中for i=1:161 if Y(i)*Y(i+1)0 a(k)=X(i); k=k+1; continue endendx=a;n=length(x);%牛顿迭代法计算方程根?E=1e-4;for i=1:n for j=1:50 f=x(i)-(6*x(i)5-45*x(i)2+20)/(30*x(i)4-90*x(i);b=abs(f-x(i);if bE %判断误差是否满足要求 breakendx(i)=f; %求得的根重新赋值给x endenddisp(计算结果为:,num2str(x);d) 计算结果与分析（1） 根的大致分布区间曲线图如下（2） 根的求解结果：（3） 分析总结由计算结果可知三个根分别为：-0.65455, 0.68118, 1.870，满足误差精度要求，在计算求解该题过程中，利用已知多项式函数的图像可以初步观察出根的大致分布区间，然后通过循环得出根区间的迭代初值，从而有效缩减了求根区间。为防止迭代进入死循环，设置了最高迭代次数为50次，而结果显示方程的三个根的迭代次数均在4次以内，可以看出牛顿迭代法的迭代求解多项式根的速度非常快，精度也高。5. 线性方程组求解。（1）编写程序实现大规模方程组的高斯消去法程序，并对所附的方程组进行求解。所附方程组的类型为对角占优的带状方程组。（2）针对本专业中所碰到的实际问题，提炼一个使用方程组进行求解的例子，并对求解过程进行分析、求解。(一)a) 解题思想高斯消去法的思想很简单，它的第一步是将n元方程组的n-1个方程实施“消元”，形成一个与原方程等价的新方程组，而被消元的n-1个方程实际上形成了一个n-1元方程组，逐步消元后得到一个上三角形状的方程组，它的从上到下的各个方程逐个减少一个未知量，最后一个方程是一元一次方程，然后，从最后一个方程解出一个未知量开始，逐次回代求得方程组的全部解。列主元消去法是当高斯消元到第k步时，从k列的以下（包括）的各元素中选出绝对值最大的，然后通过行交换将其交换到的位置上，这样就可以有效的避免大数除以小数而造成的上溢现象。根据教材提供的算法，编写高斯列主元消去法程序，编写列主元消去法的子函数与适应于超大规模超出系统内存的方程组的改编程序。b) 算法结构1.数据文件的文件名后缀为.dat,形式为：文件名.dat；2.数据文件中的数据均为二进制记录结构，因此必须使用二进制方式进行读取；3.数据文件的结构，分为以下四个部分：(1)文件标示部分,该部分用于存放数据文件的描述信息结构如下(用C语言格式进行描述)： typedef struct FileInfo long int id;/ 数据文件标示 long int ver;/ 数据文件版本号 long int id1;/ 备用标志 FILEINFO;其中：id:为该数据文件的标识，值为0xF1E1D1A0即为：十六进制的F1E1D1A0ver:为数据文件的版本号，值为16进制数据,版本号说明0x101系数矩阵为非压缩格式稀疏矩阵0x102系数矩阵为非压缩格式带状矩阵0x201系数矩阵为压缩格式稀疏矩阵0x202系数矩阵为压缩格式带状矩阵id1:为备用标志字段，暂时未用;(2)矩阵描述部分:此部分中包括矩阵的阶数和上下带宽,如果是稀疏矩阵，则上下带宽值为0。结构如下： typedef struct HeadInfo long int n; / 方程组的阶数 long int q; / 带状矩阵上带宽 long int p; / 带状矩阵下带宽 HEADINFO;(3) 系数矩阵数据部分：该部分存放方程组系数矩阵中的所有元素如存贮格式为非压缩格式,则按行方式顺序存贮系数矩阵中的每一个元素，元素总个数为n*n,每个元素的类型为float型；如果存贮格式是压缩方式,则同样是按行方式进行存贮,每行中只放上下带宽内的非零元素,即每行中存贮的元素都为p+q+1个；(4)右端系数部分：该部分存放方程组中的右端系数按顺序存贮右端系数的每个元素，个数为n个，每个系数的类型为float型4 数据文件的读取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与ver5 A的阶数6 For 1234566.1 找满足的下标6.2 For 6.2.16.36.4 For 6.4.16.4.2 For 6.4.2.16.4.3For c) 计算源程序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(读取矩阵的阶数);nq=head(4); %上带宽fprintf(读取矩阵的上带宽);qp=head(5); %下带宽fprintf(读取矩阵的下带宽);pdist=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(there is a mistake！) 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(there is a mistake！); 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(