1、西安交通大学计算方法 B 上机编程报告学号：XXX姓名：XXX专业：工程热物理 班级：硕 XXXX date2015/12/10注：本上机报告使用的程序语言均为 Matlab 语言，为本人独立完成！11. 计算以下和式： ，要求：0142168856nSn (1)若保留 11 个有效数字，给出计算结果，并评价计算的算法；(2)若要保留 30 个有效数字，则又将如何进行计算。a) 解题思想（1） 先根据精度要求估计所需累加的项数 n，使用后验误差估计方法，条件为： (m 为有效数-m014211068856nS 字位数)。（2） 在该题中 S 的和式存在两个相近的数相减的问题，为了避免有效数字损
2、失，在计算中改变了运算顺序，分别将正数和负数分别相加，然后再将其和相加。（3） 为了避免大数吃小数的问题，本题先计算出保留目标有效数字所需要的迭代次数，然后采用倒序相加的方法提高计算精度。b) 算法实现的结构 1. S1=s2=0；;142168856nann ()2. for n=0,1,2,iIf 0mnend 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; 2s2=a4+a3+a2+s2;end S=s1-s2;c) 计算源程
3、序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;endk=n;end;for n=(k-1):-1:0a1=4/(16n*(8*n+1);a2=2/(16n*(8*n+4);a3=1/(16n*(8*n+5);a4=1/(16n*(8*n+6
4、);s1=a1+s1; %第一项倒序累加s2=a4+a3+a2+s2; %后三项倒序累加endS=s1-s2; %正数累加值和负数累加值的和S=vpa(S,m) %控制 S 的精度d) 计算结果与评价当需保留 11 位有效数字时，需要将 n 值加到3n=7， s=3.1415926536； 当需保留 30 位有效数字时，需要将 n 值加到n=22, s=3.14159265358979311599796346854。 由计算结果可以看出，采用从后往前进行计算的方式，避免了“大数吃小数”的问题，这种算法很好的保证了计算结果要求保留的准确数字有效位数的要求。另外， 中有多个负数相加，按照绝对值递增
5、的顺序求和，减少了舍入误差带ns来的影响。2. 某通信公司在一次施工中，需要在水面宽度为 20 米的河沟底部沿直线走向铺设一条沟底光缆。在铺设光缆之前需要对沟底的地形进行初步探测，从而估计所需光缆的长度，为工程预算提供依据。已探测到一组等分点位置的深度数据(单位:米)如下表所示：分点 0 1 2 3 4 5 6深度 9.01 8.96 7.96 7.97 8.02 9.05 10.13分点 7 8 9 10 11 12 13深度 11.18 12.26 13.28 13.32 12.61 11.29 10.22分点 14 15 16 17 18 19 20深度 9.15 7.90 7.95 8
6、.86 9.81 10.80 10.93(1)请用合适的曲线拟合所测数据点；(2)估算所需光缆长度的近似值，并作出铺设河底光缆的曲线图；a) 解题思想该题需要对所需光缆长度进行近似估计，通过对测量数据进行拟合并进行曲线积分即可得到河底光缆长度的近似值。由于数值点较多，如果使用多项式插值，会出现龙格现象导致误差较大，因此，用相对较少的插值数据点作插值，可以避免大的误差，但又希望将所得数据点都用上，且所用数据点越多越好，所以本题采用分段插值方式，即用分段多项式代替单个多项式作插值。分段多项式是由一些在相互连接的区间上的不同多项式连接而成的一条连续曲线，其中三次样条插值方法是一种具有较好“光滑性”的
7、分段插值方法。 在本题中，假设所铺设的光缆足够柔软，在铺设过程中光缆触地走势光滑，紧贴地面，并且忽略了水流对光缆的冲击，故本题使用三次样条插值进行计算精度是最高的。4b) 算法实现的结构1 For 0,12,in1.1 iiyM2 For ,k2.1 For 1,ink2.1.1 ()/()iiikixM3 101xh4 For ,2in4.1 11iiix4.2 /();2iiiiiihcab4.3 16iiM5 000;n nndc6 11,bY7 获取 M 的矩阵元素个数，记为 m8 For 2,3k8.1 1/kkal8.2 bc8.3 1kkdlY9 /mM10 For ,2,k10
8、.1 1()/kkkYcM11 获取 x 的元素个数存入 s 512 1k13 For ,21in13.1 If then breakix;ikelse k14 1;kh;x1;kx332211()()/666kkkkkxhhMyMyxs15 计算光缆长度时，用如下公式： 20Lds1920222201()()()kfxfxdxyc) 计算源程序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,1
9、1.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);6endendh(1)=x(2)-x(1);%计算三对角阵系数 a,b,c 及右端向量 dfor i=2:n-1
10、;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=
11、n-1:-1:1;M(k)=(Y(k)-c(k)*M(k+1)/u(k);ends=zeros(1,N);7for m=1:N;k=1;for i=2:n-1if X(m)=x(i);k=i-1;break;elsek=i;endendH=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:
12、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;8d) 计算结果与评价计算所得的运行结果：铺设海底光缆的插值曲线图如所示：0 2 4 6 8 10 12 14 16 18 20789101112131
13、4/m上述插值结果表明，采用分段三次样条插值所得的拟合曲线能较准确地反映铺设光缆的走势，可以有效避免龙格现象。在本题的计算中采用自然三次样条函数作为边界条件，在解线性方程组时使用了 TSS 算法求解带状矩阵，其计算速度快，是一种求解线性方程组的有效手段，在估计河底光缆长度时使用第一类线积分，最终计算出所需光缆的长度为 L=26.4844m。在计算编程过程中对三次样条插值的原理理解不透彻，导致编程逻辑思路不清晰，调试时间较长，体会到编程是一项仔细认真的工作，不能马虎！3. 假定某天的气温变化记录如下表所示，试用数据拟合的方法找出这一天的9气温变化的规律;试计算这一天的平均气温,并试估计误差。时刻
14、 0 1 2 3 4 5 6 7 8 9 10 11 12平均气温 15 14 14 14 14 15 16 18 20 20 23 25 28时刻 13 14 15 16 17 18 19 20 21 22 23 24平均气温 31 34 31 29 27 25 24 22 20 18 17 16a) 解题思想该题要求拟合一天的气温变化规律，由于记录数据点的数目较多。需要采用很高次的多项式做数据近似，这会给计算带来困难，也会导致误差很大。用“插值样条函数”做数据近似，虽然有很好的数值性质，且计算量也不大，但存放参数 Mi 的量很大，且没有一个统一的数学公式来表示，也带来了一些不便。另一方面，
15、在有的实际问题中，用插值方法并不合适，当数据点的数目很大时，要求 通过所有数据点，可能会失去原数据所表示的规律，如果数据点是由()px测量而来的，必然带有误差，插值法要求准确通过这些不准确的数据点是不合适的。在这种情况下，不用插值标准而用选取 使 最小的标准，即采用最小i2E二乘近似方法。本题采用“最小二乘法”找出这一天的气温变化的规律，考虑使用最小二乘的二次函数、三次函数、四次函数以及指数函数 ，计算相应的系213axye数，估算误差，问题的难度是求解各种拟合函数的系数。由于利用法方程求解最小二乘系数时，方程的解不够稳定，所以本题利用正交化方法解最小二乘问题。b) 算法结构用正交化方法求数据
16、的最小二乘近似问题，假定数据以用来生成了 G，并将 y 作为其最后一列（第 n+1 列）存放，结果在 a 中，是误差是 E2。(一) 使用高次多项式拟合1 将“时刻值”存入 ，数据点的个数存入 Nix2 输入拟合多项式函数 的最高项次数 n-1，则拟合多项式函数为()p， 。12()()()npxagxagx103 根据给定数据点确定 G,并将 y 作为 n+1 列存入 GFor 1,23,jnFor im1,jiijxg,1iny4 For 2k4.1 形成矩阵 kQ4.1.1 21/sgn()Nkik4.1.2 kk4.1.3 For 1,2,j4.1.3.1 jkjg4.1.4 4.2
17、变换 到 1kGk4.2.1 gFor ,2,1jkn4.2.2 ()/Nijikgt4.2.3 For ,1,4.2.3.1 ijiijgt5 解三角方程 1Rxh5.1 ,1/nna5.2 For ,2,i5.2.1 ,1/niijiijgxg116 计算误差 2E,1minig2E(二) 使用指数函数拟合当为指数函数时，假定指数函数的形式为 。只需将 y 值求对数，213axye其主体部分不变，令 ，则可得多项式: 。 lnyz2ozc) MATLAB 源程序clear; %清除工作空间变量clc; %清除命令窗口命令x=0:24; %存入时刻值y=15 14 14 14 14 15 1
18、6 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:nfor i=1:N12G(i,j)=x(1,i)(j-1
19、);endend%形成矩阵 Qk，见课本 P123for k=1:nfor i=k:NC(k)=G(i,k)2+C(k);endC(k)=-sign(G(k,k)*(C(k)0.5); %计算 的值存入 C 中w(k)=G(k,k)-C(k); %建立 w 来存放 for j=k+1:Nw(j)=G(j,k);endB(k)=C(k)*w(k); %计算 的值并存入 B 中%变换矩阵 Gk-1 到 Gk，参考课本 P123G(k,k)=C(k);for j=k+1:n+1E=0;for i=k:NE=w(i)*G(i,j)+E;endt=E/B(k);for i=k:NG(i,j)=t*w(i
20、)+G(i,j);endendend%求解三角方程：Rx=h1a(n)=G(n,n+1)/G(n,n);13for i=n-1:(-1):1for j=i+1:nF=G(i,j)*a(j)+F;enda(i)=(G(i,n+1)-F)/G(i,i); % a(i)存放各级系数F=0;enddisp(各级系数分别为 ,num2str(a)%拟合后的回代求解过程p=zeros(1,N);for j=1:Nfor i=1:np(j)=p(j)+a(i)*x(j)(i-1);endendplot(x,p,bx,x,p,-);xlabel(时刻/h,fontsize,16); %标注坐标轴含义ylabe
21、l(平均温度 /,fontsize,16);title(拟合平均气温曲线图,fontsize,16);grid on;E2=0; %用 E2 来存放误差%计算拟合误差for i=n+1:NE2=G(i,n+1)2+E2;endE2=E20.5;disp(拟合所得误差为 E2=,num2str(E2);Tm=0;for i=1:N14Tm=Tm+p(i);endt=Tm/N;%计算平均温度disp(平均温度为 t=,num2str(t),)d) 计算结果与分析（1） 采用二次函数拟合故可得函数形式为： 2()8.3062.40.937pxxx二次函数的最小二乘曲线如下图所示：0 5 10 15
22、20 255101520253035/h/（2） 采用三次函数拟合15故可得函数形式为： 23()13.80.2739.043.0869pxxxx三次函数的最小二乘曲线如下图所示：（3） 采用四次函数拟合故可得函数形式为： 234()16.793.046.890.54180.98pxxxxx0 5 10 15 20 25101520253035/h/16四次函数的最小二乘曲线如下图所示：0 5 10 15 20 25101520253035/h/（4） 采用指数函数拟合故可得函数形式为： 22.3850.19.0416x()16.793()xpxpxe指数函数的最小二乘曲线如下图所示：170
23、5 10 15 20 252.22.42.62.833.23.43.63.8/h/通过上述拟合结果可以发现，使用最小二乘法进行数据拟合时，多项式的次数越高，拟合曲线也更加接近数据点，计算拟合的效果越好，误差越小，结果越准确。同时，指数多项式拟合的次数虽然不高，但误差最小，是因为指数型二次函数是对数据点求对数进行拟合的，故结果最准确。在进行最小二乘法编程时发现，使用正交化方法求解最小二乘问题比使用法方程方法更稳定，结果更可靠。但编程是中间变量较多容易造成混乱，计算量也较大，但本题在编程过程中有一定的适用性，只需要对数据点进行处理就可以拟合出给定次数的多项式，为以后的数据拟合工作打下了基础。 4.
24、 设计算法，求出非线性方程 的所有实根，并使误差不超52640x过 。410a) 解题思想由题目所给的待求方程可知该方程为 5 次多项式，导数的求值比较容易，可采用牛顿法迭代求解。具体实现步骤为：首先，研究函数的形态，确定根的范围；通过剖分区间的方法确定根的位置，然后利用牛顿法的基本原理进行求解，找到满足精度要求的解。18令 ，其迭代格式为52()640fxx()1kkfxx要使牛顿迭代法收敛，则必须寻找根的合理区间a,b,使得在该区间内，即有根。在选定的区间内函数 的一二阶导数 均()0fab()fx(),fx不变号。选定一个初值 ，使 ，则牛顿迭代求解过程完成。0x(1)()kkb) 算法
25、结构1 For max,2nN1.1 (0)ff1.2 If then02|1.2.1输出“ 的值已足够小”()fx1.2.2 ；stop(0)else1.2.3 (0)0fxf1.2.4 ()/x1.2.5 If then(0)1|x1.2.5.1 输出“ 已足够小 ”；stop |nxelse1.2.5.2 (0)x2 输出“ 次迭代后仍不收敛” ，stopmaNc) 计算源程序clear;clc;%根据函数曲线观察根的大致分布区间X=-8:0.1:8;for i=1:161Y(i)=6*X(i)5-45*X(i)2+20;19endplot(X,Y,r-,LineWidth,2); %画
26、出 Y 随 X 变化的曲线图%标注坐标轴含义及图像标题xlabel(X,fontsize,10);ylabel(Y,fontsize,10);title(多项式函数变化曲线图,fontsize,10);grid on;k=1;Y(162)=0;%使用二分法来确定不同根区间的迭代初值，并存于 a 中for i=1:161if Y(i)*Y(i+1)0a(k)=X(i);k=k+1;continueendendx=a;n=length(x);%牛顿迭代法计算方程根?E=1e-4;for i=1:nfor j=1:50f=x(i)-(6*x(i)5-45*x(i)2+20)/(30*x(i)4-90
27、*x(i);b=abs(f-x(i);if bE %判断误差是否满足要求breakendx(i)=f; %求得的根重新赋值给 xend20enddisp(计算结果为 :,num2str(x);d) 计算结果与分析（1） 根的大致分布区间曲线图如下-8 -6 -4 -2 0 2 4 6 8-2-1.5-1-0.500.511.52x 105XY（2） 根的求解结果：（3） 分析总结由计算结果可知三个根分别为：-0.65455, 0.68118, 1.870，满足误差精度要求，在计算求解该题过程中，利用已知多项式函数的图像可以初步观察出根的大致分布区间，然后通过循环得出根区间的迭代初值，从而有效缩
28、减了求根区间。为防止迭代进入死循环，设置了最高迭代次数为 50 次，而结果显示方程的三个根的迭代次数均在 4 次以内，可以看出牛顿迭代法的迭代求解多项式根的速度非常快，精度也高。5. 线性方程组求解。（1）编写程序实现大规模方程组的高斯消去法程序，并对所附的方程组进行求解。所附方程组的类型为对角占优的带状方程组。21（2）针对本专业中所碰到的实际问题，提炼一个使用方程组进行求解的例子，并对求解过程进行分析、求解。(一)a) 解题思想高斯消去法的思想很简单，它的第一步是将 n 元方程组的 n-1 个方程实施“消元” ，形成一个与原方程等价的新方程组，而被消元的 n-1 个方程实际上形成了一个 n
29、-1 元方程组，逐步消元后得到一个上三角形状的方程组，它的从上到下的各个方程逐个减少一个未知量，最后一个方程是一元一次方程，然后，从最后一个方程解出一个未知量开始，逐次回代求得方程组的全部解。列主元消去法是当高斯消元到第 k 步时，从 k 列的 以下（包括 ）的kaka各元素中选出绝对值最大的，然后通过行交换将其交换到 的位置上，这样就可以有效的避免大数除以小数而造成的上溢现象。根据教材提供的算法，编写高斯列主元消去法程序，编写列主元消去法的子函数与适应于超大规模超出系统内存的方程组的改编程序。b) 算法结构1.数据文件的文件名后缀为.dat,形式为：文件名.dat；2.数据文件中的数据均为二
30、进制记录结构，因此必须使用二进制方式进行读取；3.数据文件的结构，分为以下四个部分：(1)文件标示部分,该部分用于存放数据文件的描述信息结构如下(用 C 语言格式进行描述)：typedef struct FileInfo long int id; / 数据文件标示long int ver; / 数据文件版本号long int id1; / 备用标志 FILEINFO;其中：id:为该数据文件的标识，值为 0xF1E1D1A0 即为：十六进制的 F1E1D1A022ver:为数据文件的版本号，值为 16 进制数据,版本号 说明0x101 系数矩阵为非压缩格式稀疏矩阵0x102 系数矩阵为非压缩格
31、式带状矩阵0x201 系数矩阵为压缩格式稀疏矩阵0x202 系数矩阵为压缩格式带状矩阵id1:为备用标志字段，暂时未用;(2)矩阵描述部分:此部分中包括矩阵的阶数和上下带宽,如果是稀疏矩阵，则上下带宽值为 0。结构如下：typedef struct HeadInfo long int n; / 方程组的阶数long int q; / 带状矩阵上带宽long int p; / 带状矩阵下带宽 HEADINFO;(3) 系数矩阵数据部分：该部分存放方程组系数矩阵中的所有元素如存贮格式为非压缩格式,则按行方式顺序存贮系数矩阵中的每一个元素，元素总个数为 n*n,每个元素的类型为 float 型；如果
32、存贮格式是压缩方式,则同样是按行方式进行存贮,每行中只放上下带宽内的非零元素,即每行中存贮的元素都为 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
33、(2); %这两句是进行进制转换，读取 id 与 ver5 A 的阶数 n236 For 1,23,kn6.1 找满足 的下标|max|kikik6.2 For ,j6.2.1 kkjj6.3 k6.4 For 1,2,in6.4.1 iki6.4.2 For 1,2,jn6.4.2.1 ijikjij6.4.3 iii/nnxFor 1,2,k1()/nkkjkjxc) 计算源程序clear; %清除工作空间变量clc; %清除命令窗口命令%读取系数矩阵f,p=uigetfile(*.dat,选择数据文件); %读取数据文件num=5; %输入系数矩阵文件头的个数name=strcat(p,
34、f);file=fopen(name,r);head=fread(file,num,uint); %读取二进制头文件24id=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); %把句柄值转向第六个元素
35、开头处A,count=fread(file,inf,float); %读取二进制文件，获取系数矩阵fclose(file); %关闭二进制头文件%对非压缩带状矩阵进行求解if ver=102a=zeros(n,n);for i=1:nfor j=1:na(i,j)=A(i-1)*n+j); %求系数矩阵 a(i,j)endend25b=zeros(n,1);for i=1:nb(i)=A(n*n+i);endfor k=1:n-1 %列主元高斯消去法m=k;for i=k+1:n, %寻找主元if abs(a(m,k)abs(a(i,k)m=i;endendif a(m,k)=0 %遇到条件终
36、止disp(there is a mistake！)returnendfor 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;endfor i=k+1:n %计算 l(i,k)并将其放到 a(i,k)中a(i,k)=a(i,k)/a(k,k);for j=k+1:na(i,j)=a(i,j)-a(i,k)*a(k,j);26endb(i)=b(i)-a(i,k)*b(k);endendx=zeros(n,1); %回代过程x(n)=b(n)/a(n,n);for k=n-1:-1:1x(k)=(b(
37、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:nfor j=1:1:ma(i,j)=A(i-1)*m+j); %求 a(i,j)endendb=zeros(n,1);for i=1:1:nb(i)=A(n*m+i); %求 b(i)endfor k=1:1:(n-1) %开始消去过程if a(k,(p+1)=0disp(there is a mistake！);break;end27st1=n;if (k+p)nst1=k+p;endfor
38、i=(k+1):1:st1a(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);endb(i)=b(i)-a(i,k+p-i+1)*b(k);endendx=zeros(n,1); %回代过程x(n)=b(n)/a(n,p+1);sum=0;for k=(n-1):-1:1sum=b(k);st2=n;if (k+q)nst2=k+q;endfor j=(k+1):1:st2sum=sum-a(k,j+p-k+1)*x(j);end
39、x(k)=sum/a(k,p+1);sum=0;endenddisp(方程组的的解为：) %输出方程组的解28disp(x)d) 计算结果与分析（1） dat51.dat 求解结果：文件标识符为 id =F1E1D1A0文件版本号为 ver =102读取矩阵的阶数 n =15读取矩阵的上带宽 q =3读取矩阵的下带宽 p =3方程组的的解为：1.00001.00001.00001.00001.00001.00001.00001.00001.00001.00001.00001.00001.00001.00001.0000（2） dat52 求解结果文件标识符为 id =F1E1D1A0文件版本号为 ver =202读取矩阵的阶数 n =20读取矩阵的上带宽 q =529读取矩阵的下带宽 p =5方程组的的解为：1.00001.00001.00001.00001.00001.00001.00001.00001.00001.00001.00001.00001.00001.00001.00001.00001.00001.00001.00001.0000（3） dat53 求解结果文件标识符为 id =F1E1D1A0文件版本号为 ver =102读取矩阵的阶数 n =2160读取矩阵的上带宽 q =5读取矩阵的下带宽 p =5方程组的的解为：
