计算方法(B)上机作业

1、计算方法B上机实习报告学院：电信学院姓名：许梦芸班级：硕5026 学号：3115032012代课老师：苏剑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*10（1-11）保留30个有效数字，即误差0.5*10（1-3

2、0）B:算法实现的结构分为两个部分：顺序结构和循环结构。由于不知道循环次数，因此采用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

3、/(8*i+1)-2/(8*i+4)-1/(8*i+5)-1/(8*i+6);End（1）运行结果：算法同上所述%有效数字为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; F

4、x2=6*x25-45*x22+20; f=x2-(x1-x2)*Fx2/Fx1/(1-Fx2/Fx1); %割线法迭代公式endC:结果及分析初值点零点迭代次数x1=-1,x2=-0.9-0.6545423905x1=0.9,x2=10.6811741074x1=2,x2=31.8707990736由题意可知，一次导函数有两个零点，而原函数在正负无穷大分别取到正无穷和负无穷大，因此可知本题共有三个实根。分别选取初值进行逼近，可以看出迭代次数相对较小。同时迭代次数还与所选初值有关，初值离真实值越远，迭代次数越多。5.线性方程组求解。（1）编写程序实现大规模方程组的高斯消去法程序，并对所附的方程

5、组进行求解。所附方程组的类型为对角占优的带状方程组。（2）针对本专业中所碰到的实际问题，提炼一个使用方程组进行求解的例子，并对求解过程进行分析、求解。附：数据文件格式说明1.数据文件的文件名后缀为.dat,形式为：文件名.dat；2.数据文件中的数据均为二进制记录结构，因此必须使用二进制方式进行读取；3.数据文件的结构，分为以下四个部分：(1)文件标示部分,该部分用于存放数据文件的描述信息结构如下(用C语言格式进行描述)： typedef struct FileInfo long int id;/ 数据文件标示 long int ver;/ 数据文件版本号 long int id1;/ 备用标

6、志 FILEINFO; 其中：id:为该数据文件的标识，值为0xF1E1D1A0即为：十六进制的F1E1D1A0 ver:为数据文件的版本号，值为16进制数据,版本号说明0x101系数矩阵为非压缩格式稀疏矩阵0x102系数矩阵为非压缩格式带状矩阵0x201系数矩阵为压缩格式稀疏矩阵0x202系数矩阵为压缩格式带状矩阵 id1:为备用标志字段，暂时未用; (2)矩阵描述部分:此部分中包括矩阵的阶数和上下带宽,如果是稀疏矩阵，则上下带宽值为0。结构如下： typedef struct HeadInfo long int n; / 方程组的阶数 long int q; / 带状矩阵上带宽 long

7、int p; / 带状矩阵下带宽 HEADINFO;(3) 系数矩阵数据部分：该部分存放方程组系数矩阵中的所有元素如存贮格式为非压缩格式,则按行方式顺序存贮系数矩阵中的每一个元素，元素总个数为n*n,每个元素的类型为float型；如果存贮格式是压缩方式,则同样是按行方式进行存贮,每行中只放上下带宽内的非零元素,即每行中存贮的元素都为p+q+1个； (4)右端系数部分：该部分存放方程组中的右端系数 按顺序存贮右端系数的每个元素，个数为n个，每个系数的类型为float型3.数据文件说明：(1)Dat51.dat 为非压缩带状方程组，阶数为15阶，该方程组供调试程序使用,该方程组的根都为1;(2)D

8、at52.dat 为压缩带状方程组，阶数为20阶，该方程组供调试程序使用,该方程组的根都为1;(2)Dat53.dat 为非压缩带状方程组，阶数在2000阶左右;(3)Dat54.dat 为压缩带状方程组，阶数在40000阶左右;解答：A:算法实现的思想及依据和算法实现的结构高斯消去法主要分为两大步骤，即消去过程和回带过程。算法借鉴课本GAUSSPP算法。由于题目中给出的系数矩阵是对角占优的矩阵，因此可以不选主元直接进行高斯消去；另外在非压缩格式带状矩阵中，存在着大量的非零元素，非零元素的运算毫无意义并且占用大量机器资源和时间，因此对课本中给出的GAUSSPP算法进行优化。对于n阶、上带宽q、

9、下带宽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的文件数据操作函数完成。要用

10、到的文件函数如下：fopen文件打开函数，文件格式选择r（即二进制只读格式），同时获得文件指针fid(之后文件操作要用fid)；再用fread读取文件信息，函数中按照long类型读取文件头5个长整形数据。判断文件格式：将文件版本号信息数据(10进制)读出，利用dec2hex命令将10进制数据转换为16进制的字符串，然后将此字符串与202和102比较，如果与202相等则压缩格式带状矩阵，如果与102相等即说明其为非压缩格式带状矩阵。系数矩阵和右端向量读取时只能存放于一维数组中，因此需要根据格式要求将读取到的系数矩阵数据转换到相应的矩阵格式。将原始的dat51和dat52文件带入到编写好的程序所在

11、文件夹(即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)

12、; %文件版本号转换为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:n

13、A(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

14、);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

15、) %输出小数点后保留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(待求解矩阵为非压

16、缩带状矩阵,数据与题目要求不符)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); %求取

17、右端向量元素个数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)

18、; 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

19、位数的浮点数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(待求解矩阵为压缩带状矩阵,数据与题目要求不符)end

20、n = 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); %读取数据结束

21、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)=

22、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+);

23、 %以读写方式打开指定文件，将文件指针指向文件末尾fprintf(fid,%s的所有解,filename);fprintf(fid,%5.5f ,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

24、进制对应字符串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 =

25、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