西安交大计算方法B大作业Word格式.docx

1、（2）估算所需光缆长度的近似值，并作出铺设河底光缆的曲线图;1.2实现题目的思想及算法依据首先在题目（1）中要实现的是数据的拟合，显然用到的是我们在第三章中数 据近似的知识内容。多项式插值时，这里有 21个数据点，则是一个20次的多 项式，但是多项式插值随着数据点的增多， 会导致误差也会随之增大，插值结果 会出现龙格现象，所以不适用于该题目中点数较多的情况。 为了避免结果出现大 的误差，同时又希望尽可能多地使用所提供的数据点， 提高数据点的有效使用率， 这里选择分段插值方法进行数据拟合。分段插值又可分为分段线性插值、分段二 次插值和三次样条插值。由于题目中所求光缆的现实意义，而前两者在节点处的

2、 光滑性较差，因此在这里选择使用三次样条插值。根据课本SPLINEM算法和TSS算法，采用第三种真正的自然边界条件，在 选定边界条件和选定插值点等距分布后，可以先将数据点的二阶差商求出并赋值 给右端向量d，再根据TSS解法求解三对角线线性方程组从而解得 M值。求出 M后，对区间进行加密，计算200个点以便于绘图以及光缆长度计算。对于问题（2），使用以下的公式（x ）ds1.3算法结构1. For i 0,1,2, , n1.1yiM iFork1,22.1in, n1, ,k2.1.1 （Mi Mi 1）/（ Xi Xixh11,2,n-4.1X 1Xhi 14.2hi 1/（ hihi 1）

3、Ci;1 C4.36MididoM0;dnMn;0 C0b0;nan;bn1, d12,3,m!获取17.1aJk 1lk7.2bk-lk C7.3dk -m/mMrm 1, m 2,19.1（kCkMk 1）/k M kFor i12,s 111.1ifXi then i kbreakai；2 bk） MiM的矩阵元素个数，存入 m获取x的元素个数存入s2.3.4.5.6.7.8.9.10.11.12. xkelse i 1 kxk1 一 h; x x x; x x=x（j）&x1（i）%误差允许值%初始区间%判断方程是否有解disp（该方程在a,b上无解elseif fO=Ox=a;els

4、eif f1=0x=b; %判断方程解是否在区间两边界上else%二分法求解方程得解a0=a;b0=b;while abs（b0-a0）/2）=errorhalf=（a0+b0）/2;fa=hsz（aO）;fb=hsz（bO）;fhalf=hsz（half）; %计算中点处的函数值，用以判断解的位置if fhalf=0x=half;break;elseif fa*fhalfb0=half; %定义新区间，为原区间的一半a0=half; %定义新区间x=（b0+a0）/2; %方程组的解方程组的解为：）disp（x）2.5结果与说明n的取值是很关由于是利用梯形公式来求解，需要将区间划分为 n个区

5、间，而 键的，需要取得适当的n值才能满足误差精度。 th2方程的解为； 2.4048这里将区间分成1OOOO份，可以得到结果图2.1 , x=2.4048图2. 1 n=10000 方程的解由于变量的嵌套很多及自身水平的不足，所以我将很多变量通过调用子程序的方 式实现，分别为 theta.m : g.m; hsz.m 。3.1题目内容：编写程序实现大规模方程组的高斯消去法程序，并对所附的方程组进行求解。 所附方程组的类型为对角占优的带状方程组。数据说明：dat20171.dat 为非压缩带状方程组，阶数为 10阶，该方程组供调试程序使用，该方程组的根都为1;（2）dat20172.dat 为压

6、缩带状方程组，阶数为 20阶，该方程组供调试程序使用，该方程组的根都为1;dat20173.dat 为非压缩带状方程组，阶数在 3000阶左右；dat20174.dat 为压缩带状方程组，阶数在 50000阶左右；3.2实现题目的思想及算法依据高斯消去法主要分为消去和回代过程， 消去过程形成上三角矩阵，回代过程 就是自下而上实现方程组求解。题目中系数矩阵是严格对角占优矩阵， 所以无需 进行列主元，极大的提高了程序的效率。对于n阶、上带宽q、下带宽p的非压缩格式带状矩阵，通过比较 k+p、k+q 及n的大小，对于行消去，选择k+p和n中的最小值作为循环终点。在用第 k 行进行消去时，只需要计算第

7、 k列的ak+1，k到3k+p，k，每一行只需计算ak，k+1到 ak，k+q ；从第n-c+1到n行，第k列计算ak+1, k到an，k，每一行只需计算ak，k+1到 3k，n，于是减小运算量。压缩格式忽略了零元素，存储方式与非压缩格式有所不同，所以需要建立两种格 式中元素之间的对应关系，即压缩 Ak，p+1 =非压缩人，k，根据此关系进行程序编 写。3.3算法结构1.A的阶数 nFor i k 1,k 2,n2.2 For j k 1,k 2,n2.3nnXn3.For k n 1,n 2,1（k kjXj）/ kk xkj k 13.4 Matlab 源程序%首先计算非压缩带状方程组fn

8、ame=dat20171.dat：%读取文件，可改为 dat20173.dat目标文件为:disp（fname）;fid=fope n（ fname,r%二进制模式读取数据fhead=fread（fid,5,l on g%以long int形式读取文件前5个数据，获取系数矩阵信息bbh=dec2hex（fhead（2）;%将版本号转换为16进制%判断所读取的文件是否为目标文件if strcmp（bbh,1O2%比较版本号，相同则说明为非压缩带状矩阵目标文件系数矩阵为非压缩带状矩阵，读取正确 读取错误n=fhead （3）;%读取矩阵阶数nq=fhead（4）;%读取上带宽qp=fhead（5）

9、;%读取下带宽p系数矩阵阶数为:disp（ n）上带宽为:disp（q）下带宽为:d=fread（fid, nA2,float b=fread（fid ,n ,%生成系数矩阵A=zeros （n,n）;k=1;A（i,j）=d（k）;（i,j）k=k+1;fclose（fid）;%非压缩格式，需要读取 nA2个浮点数,%读取右端系数%因为系数矩阵按行方式顺序存储在 d%读取结束n=len gth（b）;按行方式顺序存储中，所以依次形成 A%高斯消去法，书本 33页到34页算法GAUSSPP%消去过程 n-1for i=k+1: min（ k+p ,n）%仅对不为零的区域作高斯消去， 边界时k+

10、p有可能大A（i,k）=A（i,k）/A（k,k）; %用第 k 行消去for j=k+1: min（ k+q ,n）%只需处理右边k+q个数据，边界时k+q可能大于nA（i,j）=A（i,j）-A（i,k）*A（k,j）; % 更新各列b（i）=b（i）-A（i,k）*b（k）; %更新右端系数 b 列%回代过程x（ n）=b（ n）/A（ n,n）;for k=n-1:sum=0;sum=sum+A（k,j）*x（ j）;x（k）=（b（k）-sum）/A（k,k）;end m=10;方程的前d个根:n,m）; %输出前m个根fprintf（1,%0.5f ,x（j）; %小数点后保留5位

11、小数%将方程出的解输出为txt文件if strcmp（fname,dat20173.datfile=strcat（f:,fname, 矩阵方程组的解,.txt %生成要创建文件的文件路径和文件名fid=fope n（file,a+ %以读写方式打开指定文件，将文件指针指向文件末尾fprintf（fid,%s 的所有解 n,fname）;,x）; %向指定文件写入数据，保留 5位小数end fprin tf（%第二步处理压缩带状方程组，即数据文件（ 2 ）和（4），（ 2 ）供调试使用dat20172.dat%读取文件，可改为 dat20174.dat2O2%比较版本号，相同则说明文件数据是压缩

12、带状矩阵目标文件系数矩阵为压缩带状矩阵，读取正确系数矩阵阶数为 disp（ n）d=fread（fid, n*（p+q+1）,%压缩格式，以n*（p+q+1）矩阵存储，同样按行存储b=fread（fid ,n ,A=zeros （n, p+q+1）; k=1;p+q+1 %因为系数矩阵按行方式顺序存储在 d中，所以依次形成 A % 读取结束min（p,n-k） %用第 k 行消去A（k+i,p+1-i）=A（k+i,p+1-i）/A（k,p+1）; %压缩格式 p+1 列即为原来第 k 列min（q,n-k） %只需处理右边 q个数据，但边界时 q可能大于（n-k）A（k+i,p+1-i+j）

13、=A（k+i,p+1-i+j）-A（k+i,p+1-i）*A（k,p+1+j）;b（k+i）=b（k+i）-b（k）*A（k+i,p+1-i）;x（ n）=b（ n）/A（ n,p+1）;for i=n-1:sum=O;mi n（q, n-i）sum=sum+x（i+j）*A（i,p+1+j）;x（i）=（b（i）-sum）/A（i,p+1）;m=20;方程的前%d个根:dat20174.dat3.5结果与说明目标dat目标女件系埶年阵淘非压缩帯状拒阵读取IF确系数世解阶數为 10上带宽为： 3下帯克为. 3方程的前汕吓更1.00000 1, COOOO 1. OOCOO L000D0 LOO

14、OCO目标文件为：datUl, v. dat目标交件系数柜阵为圧缩帚我矩陳，读則正确垂皴柜解撕數为 201.000001. D0000L dOOOO L 00000下带竞为：方程的前创个艰1, 00000 l.COOOO L. 00000 1. C0050 1. OOOCO1. 000001. 00000 1.00000 1. 04000首先使用dat20171和dat20172 验证程序的正确性,如图3.1图3. 1验证dat20171 和dat20172 的正确性由上图可知dat20171.dat中为上带宽为3，下带宽为3，阶数为10的非压缩矩阵，解都为dat20172.dat中为上带宽为3，下带宽为3，阶数为20的压缩矩阵，解都为1与题目中所给的条件相符图3. 2 dat20173 和dat20174的方程组解目棕文件为:dat201T3.dat目标文件系数矩阵为非压縮带状矩0车读现止确慕埶拒晦翫數为: 29E0上帯宽为： 4下带竟为：方程的前1Q吓很.1.01900 1.01300 1. 01900 L 01900 1. 01900 1. 01S&0 1. 01900 1. 01900 1. 019&0 1. 01900目标