计算方法(B)上机作业要点

?计算方法B?上机实习报告学院：电信学院姓名：许梦芸班级：硕5026 学号：3115032021代课老师：苏剑1.对以下和式计算：，要求：(1)假设只需保存11个有效数字，该如何进行计算；(2)假设要保存30个有效数字，那么又将如何进行计算；解答：A:算法实现的思想及依据首先对问题分析可知：当一般项116n(48n+1-28n+4-18n+5-18n+6)S01,这说明S要保存m个有效数字，那么小数点后位数应该取到m-1位，因此误差上限e=0.110-（m-1）。由题可得S(1/16n)*(4/(8*n+1)保存11个有效数字，即误差0.5*101-11保存30个有效数字，即误差0.5*101-30B:算法实现的结构分为两个局部：顺序结构和循环结构。由于不知道循环次数，因此采用if循环，利用通项是否大于等于误差上限来判断循环条件的真伪。Matalab中可以利用vpa函数控制最终结果的有效数位。C:源程序及相关的注释说明有效数字为11的情况digita=11; %有效数字位数for i=1:1:1000 S0=(1/(16(i-1)*(4/(8*i-7); S1=(1/16i)*(4/(8*i+1); if(S0-S1)0.5*10(1-digita) n=i; break; endendS=0;digits(digita); for i=0:1:n S=S+vpa(1/16i)*(4/(8*i+1)-2/(8*i+4)-1/(8*i+5)-1/(8*i+6);End1运行结果：算法同上所述%有效数字为30的情况digita=30;% 有效数字位数for i=1:1:1000 S0=(1/(16(i-1)*(4/(8*i-7); S1=(1/16i)*(4/(8*i+1); if(S0-S1)=x(j)&x1(i)=e %判断误差相邻两个近似值误差是否小于给定误差 temp=x2; x2=Fx(x1,x2); x1=temp; i=i+1; end xx=x2;end %采用割线法求解function f=Fx(x1,x2) Fx1=6*x15-45*x12+20; Fx2=6*x25-45*x22+20; f=x2-(x1-x2)*Fx2/Fx1/(1-Fx2/Fx1); %割线法迭代公式endC:结果及分析初值点零点迭代次数5x1=0.9,x2=14x1=2,x2=36由题意可知，一次导函数有两个零点，而原函数在正负无穷大分别取到正无穷和负无穷大，因此可知此题共有三个实根。分别选取初值进行逼近，可以看出迭代次数相对较小。同时迭代次数还与所选初值有关，初值离真实值越远，迭代次数越多。5.线性方程组求解。1编写程序实现大规模方程组的高斯消去法程序，并对所附的方程组进行求解。所附方程组的类型为对角占优的带状方程组。2针对本专业中所碰到的实际问题，提炼一个使用方程组进行求解的例子，并对求解过程进行分析、求解。附：数据文件格式说明1.数据文件的文件名后缀为.dat,形式为：文件名.dat；2.数据文件中的数据均为二进制记录结构，因此必须使用二进制方式进行读取；3.数据文件的结构，分为以下四个局部：(1)文件标示局部,该局部用于存放数据文件的描述信息结构如下(用C语言格式进行描述)： typedef struct FileInfo long int id;/ 数据文件标示 long int ver;/ 数据文件版本号 long int id1;/ 备用标志 FILEINFO; 其中：id:为该数据文件的标识，值为0xF1E1D1A0即为：十六进制的F1E1D1A0 ver:为数据文件的版本号，值为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型3.数据文件说明：(1)Dat51.dat 为非压缩带状方程组，阶数为15阶，该方程组供调试程序使用,该方程组的根都为1;(2)Dat52.dat 为压缩带状方程组，阶数为20阶，该方程组供调试程序使用,该方程组的根都为1;(2)Dat53.dat 为非压缩带状方程组，阶数在2000阶左右;(3)Dat54.dat 为压缩带状方程组，阶数在40000阶左右;解答：A:算法实现的思想及依据和算法实现的结构高斯消去法主要分为两大步骤，即消去过程和回带过程。算法借鉴课本GAUSSPP算法。由于题目中给出的系数矩阵是对角占优的矩阵，因此可以不选主元直接进行高斯消去；另外在非压缩格式带状矩阵中，存在着大量的非零元素，非零元素的运算毫无意义并且占用大量机器资源和时间，因此对课本中给出的GAUSSPP算法进行优化。对于n阶、上带宽q、下带宽p的带状矩阵，选取p与q较大者(赋值给c)，在第1行到第n-c行，第k列，只需要计算到，每一行也只需从到(i从k+1到k+q)；在第n-c+1到n行，第k列，计算到，每一行只需要从计算到(i从k+1到n)，这样可以减小运算量。对于压缩带状矩阵，其消去过程和非压缩带状矩阵根本一致，不同之处在于：压缩格式忽略了零元素，因此需要建立压缩格式带状矩阵与非压缩格式带状矩阵的对应关系。主元素对应关系： B(k,p+1)=A(k,k)B是压缩格式带状矩阵，A是非压缩格式带状矩阵，编写程序时需要根据此关系建立其元素间的对应关系。在对文件进行操作的时候，可以利用matlab的文件数据操作函数完成。要用到的文件函数如下：fopen文件翻开函数，文件格式选择r即二进制只读格式，同时获得文件指针fid(之后文件操作要用fid)；再用fread读取文件信息，函数中按照long类型读取文件头5个长整形数据。判断文件格式：将文件版本号信息数据(10进制)读出，利用dec2hex命令将10进制数据转换为16进制的字符串，然后将此字符串与202和102比拟，如果与202相等那么压缩格式带状矩阵，如果与102相等即说明其为非压缩格式带状矩阵。系数矩阵和右端向量读取时只能存放于一维数组中，因此需要根据格式要求将读取到的系数矩阵数据转换到相应的矩阵格式。将原始的dat51和dat52文件带入到编写好的程序所在文件夹(即matlab运行目录)下面，matlab运行程序会自动搜索到相应文件，观察运行结果与给定结果是否一致，假设不一致那么修改程序，否那么运行程序分别读取dat53和dat54，最终利用fprintf函数将运算结果保存到txt文件中。B: 源程序及相关的注释说明%处理文件dat51.datfname = dat51.dat;disp(文件名:);disp(fname);fid=fopen(fname,r); % r默认二进制模式只读数据fhead = fread(fid,5,long); %以long int形式读取文件前5个数据，获取系数矩阵信息geshi = dec2hex(fhead(2); %文件版本号转换为16进制对应字符串if strcmp(geshi,102) %比拟字符串，如果符合说明文件数据格式正确 disp(待求解矩阵为非压缩带状矩阵)elsedisp(待求解矩阵为压缩带状矩阵,数据与题目要求不符)endn = fhead(3); %读取矩阵的阶数q = fhead(4); %读取上带宽p = fhead(5); %读取下带宽d = fread(fid,n2,float); %非压缩格式，需要读取n2个浮点数，以一维格式存储到中间变量dm = fread(fid,n,float); %再读取右端向量的n个元素k = 1;for i = 1:nfor j = 1:nA(i,j) = d(k); k = k+1;endend %将读取到的数据放到阶数为n，上带宽为q，下带宽为p的系数矩阵中fclose(fid); %读取数据结束n = length(m); %求取右端向量的元素个数t = max(q,p); %比拟上下带宽的大小，将最大者赋值给变量cfor k=1:n-t %消去过程，由于消去进行到第n-c步之后，循环变量的上限、下限值不一致，因此需要分别讨论for i=k+1:k+pA(i,k)=A(i,k)/A(k,k);for j=k+1:k+q A(i,j)=A(i,j)-A(i,k)*A(k,j); endm(i)=m(i)-A(i,k)*m(k);endendfor k = n-t+1:n-1for i = (k+1):n 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 m(i)=m(i)-A(i,k)*m(k);endendX(n)=m(n)/A(n,n); %回带过程求解方程组的根for t=n-1:-1:1H=m(t);for i=t+1:nH=H-A(t,i)*X(i);endX(t)=H/A(t,t);enddisp(方程的前5个根为: ); %输出5个根，用于与正确解比照for j = 1:5fprintf(%5.5f ,X (j) %输出小数点后保存5位数的浮点数end%处理文件dat52.datfilename = dat52.dat;disp(文件名:);disp(filename);fid=fopen(filename,r); % r默认二进制模式只读数据inform = fread(fid,5,long); %以long int形式读取文件前5个数据，获取系数矩阵信息geshi = dec2hex(inform(2); %文件版本号转换为16进制对应字符串if strcmp(geshi,202) %比拟字符串，如果符合说明文件数据格式正确disp(待求解矩阵为压缩带状矩阵)else disp(待求解矩阵为非压缩带状矩阵,数据与题目要求不符)endn = inform(3); %读取矩阵的阶数q = inform(4); %读取上带宽p = inform(5); %读取下带宽d = fread(fid,n*(p+q+1),float); %压缩格式一共要读取n*(p+q+1)个数m = fread(fid,n,float); %再读取右端向量的n个元素t = 1;for i = 1:n %将读取到的数据放到n行、p+q+1列的系数矩阵中for j = 1:(q+p+1)A(i,j) = d(t);t = t+1;endendfclose(fid); %读取数据结束n = length(m); %求取右端向量元素个数for t=1:n-qfor i=1:pA(t+i,p+1-i)=A(t+i,p+1-i)/A(t,p+1);for j=1:qA(t+i,p+1-i+j)=A(t+i,p+1-i+j)-A(t+i,p+1-i)*A(t,p+1+j);endm(t+i)=m(t+i)-m(t)*A(t+i,p+1-i);endendfor k = n-q+1:n-1for i = 1:n-kA(k+i,p+1-i)=A(k+i,p+1-i)/A(k,p+1); for j = 1:n-k A(k+i,p+1-i+j)=A(k+i,p+1-i+j)-A(k+i,p+1-i)*A(k,p+1+j); end m(k+i)=m(k+i)-A(k+i,p+1-i)*m(k);endendX(n)=m(n)/A(n,p+1);for t=n-1:-1:n-q+1S=m(t);for i = t+1:n S=S-A(t,p+1+i-t)*X(i);end X(t)=S/A(t,p+1);endfor t=n-q:-1:1S=m(t);for j = t+1:t+qS=S-A(t,p+1+j-t)*X(j);endX(t)=S/A(t,p+1);enddisp(方程的前5个根:); %输出5个根，用于与正确解比照for j = 1:5fprintf(%5.5f ,X (j) %输出小数点后保存5位数的浮点数end%处理文件dat53.datfname = dat53.dat;disp(文件名:);disp(fname);fid=fopen(fname,r); % r默认二进制模式只读数据inform = fread(fid,5,long); %以long int形式读取文件前5个数据，获取系数矩阵信息geshi = dec2hex(inform (2); %文件版本号转换为16进制对应字符串if strcmp(geshi,102) %比拟字符串，如果符合说明文件数据格式正确 disp(待求解矩阵为非压缩带状矩阵)else disp(待求解矩阵为压缩带状矩阵,数据与题目要求不符)endn = inform(3); %读取矩阵的阶数q = inform(4); %读取上带宽p = inform(5); %读取下带宽d = fread(fid,n2,float); %非压缩格式，需要读取n2个浮点数，以一维格式存储到中间变量dm=ones(1,n);m = fread(fid,n,float); %再读取右端向量的n个元素A=ones(n,n); %预定义矩阵Ak = 1;for i = 1:nfor j = 1:nA(i,j) = d(k); k = k+1;endend %将读取到的数据放到阶数为n，上带宽为q，下带宽为p的系数矩阵中fclose(fid); %读取数据结束n = length(m); %求取右端向量元素个数c = max(q,p); %比拟上下带宽的大小，将最大者赋值给变量cfor k=1:n-c %消去过程，由于消去进行到第n-c步之后，循环变量的上限、下限值不一致，因此需要分别讨论for i=k+1:k+pA(i,k)=A(i,k)/A(k,k);for j=k+1:k+qA(i,j)=A(i,j)-A(i,k)*A(k,j);endm(i)=m(i)-A(i,k)*m(k);endendfor k = n-c+1:n-1for i = (k+1):n 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 m(i)=m(i)-A(i,k)*m(k);endendX(n)=m(n)/A(n,n); %回带过程求解方程组的根for t=n-1:-1:1H=m(t);for j=t+1:nH=H-A(t,j)*X(j);endX(t)=H/A(t,t);enddisp(方程的前5个根:);for i = 1:5fprintf(%5.5f ,X (i) %在matlab的命令窗口显示前5个根endfile = strcat(f:,filename,矩阵方程组的解,.txt); %生成要创立文件的文件路径和文件名fid = fopen(file,a+); %以读写方式翻开指定文件，将文件指针指向文件末尾fprintf(fid,%s的所有解,filename);fprintf(fid, ,X); %向指定文件写入数据，保存5位小数fclose(fid); %结束文件操作%处理文件dat54.datfname = dat54.dat;disp(处理的文件名:);disp(fname);fid=fopen(fname,r); % r默认二进制模式只读数据fhead = fread(fid,5,long); %以long int形式读取文件前5个数据，获取系数矩阵信息geshi = dec2hex(fhead(2); %文件版本号转换为16进制对应字符串if strcmp(geshi,202)disp(待求解矩阵为压缩带状矩阵)else disp(待求解矩阵为非压缩带状矩阵,数据与题目要求不符)endn = fhead(3); %读取矩阵的阶数q = fhead(4); %读取上带宽p = fhead(5); %读取下带宽u=n*(q+p+1);d = fread(fid, u,float); %压缩格式一共要读取n*(p+q+1)个数m = fread(fid,n,float); %再读取右端向量的n个元素m=ones(1,n);A=ones(n,q+p+1); %预定义矩阵Ak = 1;for i = 1:nfor j = 1:(u)A(i,j) = d(k);k = k+1;endend %将读取的数据存储到系数矩阵fclose(fid); %读取数据结束n = length(m); %求取右端向量元素个数for k=1:n-qfor i=1:pA(k+i,p+1-i)=A(k+i,p+1-i)/A(k,p+1); for j=1:qA(k+i,p+1-i+j)=A(k+i,p+1-i+j)-A(k+i,p+1-i)*A(k,p+1+j);endm(k+i)=m(k+i)-A(k+i,p+1-i)*m(k);end endfor k = n-q+1:n-1for i = 1:n-kA(k+i,p+1-i)=A(k+i,p+1-i)/A(k,p+1);for j = 1:n-kA(k+i,p+1-i+j)=A(k+i,p+1-i+j)-A(k+i,p+1-i)*A(k,p+1+j);endm(k+i)=b(k+i)-m(k+i,p+1-i)*m(k);endendX(n)=m(n)/A(n,p+1);for t=n-1:-1:n-q+1S=m(t);for j = t+1:n S=S-A(t,p+1+j-t)*X(j);end X(t)=S/A(t,p+1);endfor t=n-q:-1:1S=m(t);for i = t+1:t+qS=S-A(t,p+1+i-t)*X(i);endX(t)=S/A(t,p+1);enddisp(方程的前5个根为:);for i = 1:5fprintf(%5.5f ,X (i) %在matlab的命令窗口显示前5个根endfile = strcat(f:,filename,矩阵方程组的解,.txt); %生成要创立文件的文件路径和文件名fid = fopen(file,a+); %以读写方式翻开指定文件，将文件指针指向文件末尾fprintf(fid,%s的所有解,filename);fprintf(fid, ,X); %向指定文件写入数据，保存5位小数fclose(fid); %结束文件操作C:运行结果分析及总结在编写此程序中，出过3次错误，分别是：fread函数格式调用错误；压缩格式与非压缩格式带状矩阵之间的元素对应关系弄错；dat54非压缩格式中，没有对矩阵A预定义，循环条件下将读取的系数矩阵数据d向矩阵A转换时直接对A进行不断扩充，导致这一步运行很慢，依次解决了这3个问题。运行文件和时，获得解均为1，与给出的正确解一致，说明程序正确。再运行和，发现dat53阶数低，但是用时比dat54要多，说明压缩格式能够大大减小文件操作时间，缩短运行时间。古今名言敏而好学，不耻下问孔子业精于勤，荒于嬉；行成于思，毁于随韩愈兴于?诗?，立于礼，成于乐孔子己所不欲，勿施于人孔子读书破万卷，下笔如有神杜甫读书有三到，谓心到，眼到，口到朱熹立身以立学为先，立学以读书为本欧阳修读万卷书，行万里路刘彝黑发不知勤学早，白首方悔读书迟颜真卿书卷多情似故人，晨昏忧乐每相亲于谦书犹药也，善读之可以医愚刘向莫等闲，白了少年头，空悲切岳飞发奋识遍天下字，立志读尽人间书苏轼鸟欲高飞先振翅，人求上进先读书李苦禅立志宜思真品格，读书须尽苦功夫阮元非淡泊无以明志，非宁静无以致远诸葛亮熟读唐诗三百首，不会作诗也会吟孙洙?唐诗三百首序?书到用时方恨少，事非经过不知难陆游问渠那得清如许，为有源头活水来朱熹旧书不厌百回读，熟读精思子自知苏轼书痴者文必工，艺痴者技必良蒲松龄 声明访问者可将本资料提供的内容用于个人学习、研究或欣赏，以及其他非商业性或非盈利性用途，但同时应遵守著作权法及其他相关法律的规定，不得侵犯本文档及相关权利人的合法权利。谢谢合作！