5.4.2 Matlab运行结果分析
对于任意一个方程组，所有根都相等，因此Matlab求解结果如下表所示：
方程 阶数 根 Dat51.dat 15 1.0000 Dat53.dat 2160 3.1400 Dat52.dat 20 1.0000 Dat54.dat
43240
5.2100
43
由上面结果可以看出，所编写程序可以解决正常系数矩阵方程组，同时能解决带状系数矩阵方程组，并且在解决该类系数矩阵时具有很高的效率，说明了列主元消去法在这一类方程组的有效性。
如果要继续提高计算的精度，可以采用全主元消去法。但是此时的工作量要大的多。
5.5 实际问题
在实际工程技术中所遇到的电路，大多数是很复杂的，这些电路是由电器元件按照一定方式互相连接而构成的网络。在电路中，含有元件的导线称为支路，而三条或三条以上的支路的会合点称为节点。本题我们以济南教育学院学报上的一篇论文为例，利用Gauss消去法进行复杂电路的计算。
其中，R1?R2?10?,R3?4?,R4?R5?8?,R6?2?,US3?10V,US6?40V 问题的解答方法如下: 1、问题分析：
对于这类问题的计算，通常采用基尔霍夫（Kirchhoff）定律来解决。设各节点的电流如图所示，则由基尔霍夫第一定律（简记为KCL）回路电压的关系，任何一个闭合回路的电流服从希尔霍夫电压定律：沿某个方向环绕回路一周的所有电压降U的代数和等于沿同一方向环绕该回路一周的电源电压的代数和,本题用支路电流法求解各支路电流。
2、解：假设 各支路电流方向如图所示，选三个网孔作为回路，回路绕行方向均为顺时针方向。
列支路电流方程如下： 据KCL可得：
?i1?i2?i6?0 ?i2?i3?i4?0
?i4?i5?i6?0 44
据KVL可得：
10i1?10i2?4i3??20 ?10i2?8i4?2i6??40
?4i3?8i4?8i5?20上面的线性方程组可表示为下列矩阵形式：
10??1?0?11??000? ?4?1010?0?100??0?4?0
3、Matlab求解：
0?1?i1???i??10?0??2??i3??11???1?????i4???00??0i5???80??2????88?0??i6???0?0?0??0? ?20?40?20?利用列主元高斯消去法函数：
%列主元高斯消去法
function gauss(A,b) [m,n]=size(A); %增广阵 Ab=[A b];
for k=1:n-1 max=k; %寻找主元 for i=k+1:n
if Ab(i,k)>=Ab(max,k) max=i; end end
%交换主元所对应的行 temp=Ab(k,k:n+1);
Ab(k,k:n+1)=Ab(max,k:n+1); Ab(max,k:n+1)=temp; %高斯消元
for i=k+1:n
Ab(i,k)=-Ab(i,k)/Ab(k,k);
Ab(i,k+1:n+1)=Ab(i,k+1:n+1)+Ab(i,k)*Ab(k,k+1:n+1); end end Ab
%回代解方程
45
x=zeros(n,1);
x(n,1)=Ab(n,n+1)/Ab(n,n); for i=n-1:-1:1 sum=0; for j=i+1:n
sum=sum+x(j,1)*Ab(i,j); end
x(i,1)=(Ab(i,n+1)-sum)/Ab(i,i); end
%显示方程组的解
disp('方程组的解：') for xxp=1:n
fprintf('x[%d]=%f\\n',xxp,x(xxp)) end
求解的主程序： zyt.m clear; clc;
A=[-1,1,0,0,0,1; 0,-1,1,1,0,0; 0,0,0,-1,1,-1; 10,10,4,0,0,0; 0,-10,0,-8,0,2; 0,0,-4,8,8,0]; b=[0;0;0;-20;-40;20]; gauss(A,b) 结果截图：
46
Matlab在电气工程中有着重要的应用。例如：可用于直流电路分析、交流电路分析及动态电路计算，也能达到令人满意的效果。在高斯主元消去法中，尤其缩小了很大的运算量。
5.6
分析总结
Gauss 消去法要能顺利的进行，必须保证消去过程中的每一步都满足主元不为零，否则 Gauss 消去法将无法进行下去；另一方面，即使在整个消去过程中，但若在某一步，对角线上元的绝对值比其下方的某些元素的绝对值小得 将是一个绝对值很大的数，这可能引起上溢而中断计算，或带来大的舍入误差。所以在 Gauss 消去过程中，适当交换方程的顺序对保证消去过程能顺利进行及计算解的精确度都是必要的。列主元消去法交换方程的原则是使绝对值最大的成为第 k 步消去的主元，这样可以最大可能的避免这些现象，最大限度保证计算结果的精确性。
47
计算方法B上机作业报告
姓名：姚威 班级：硕5081班 学号：3115312022 学科专业：电气工程
任课教师：苏剑老师
2015年12月
目录
1. 上机题一 ................................................... 1 1.1 解题源程序.............................................. 1 1.2 计算结果................................................ 3 1.3 实验结果分析及总结...................................... 3 2. 上机题二..................................................... 3
2.1 三次样条插值拟合数据点.................................. 4
2.1.1算法组织 .......................................... 4 2.1.2 算法.............................................. 5 2.1.3算法的Matlab实现 ................................. 6 2.1.4 实验结果及分析.................................... 9 2.2 复化Simpson公式近似计算光缆长度....................... 12
2.2.1算法组织 ......................................... 12 2.2.2算法的Matlab实现 ................................ 13 2.2.3实验结果及分析 ................................... 14
3. 上机题三.................................................... 14
3.1 算法思想............................................... 14 3.2 算法组织............................................... 14 3.3 算法的Matlab实现...................................... 15 3.4 实验结果及分析......................................... 17
3.4.1法方程最小二乘法输出结果 ......................... 17 3.4.2 基于最小二乘法ployfit函数法................... 19 3.5 QR分解 ................................................ 24
3.5.1 结果分析......................................... 29
4. 上机题四.................................................... 30
4.1 二分法算法组织......................................... 30
4.1.1算法思想及依据 ................................... 30 4.1.2算法结构 ......................................... 30 4.2 二分法算法的Matlab实现................................ 30
4.2.1利用Matlab校验根的区间 .......................... 30 4.2.2使用二分法Matlab程序 ............................ 31 4.3 二分法实验结果及分析................................... 32 4.4 割线法算法组织......................................... 33
4.4.1算法思想及依据 ................................... 33 4.4.2算法结构 ......................................... 33 4.5 割线法算法的Matlab实现................................ 33
4.5.1割线法实验结果及分析 ............................. 34
5. 上机题五.................................................... 34
5.1 算法思想............................................... 35 5.2 算法组织............................................... 35
5.2.1算法思想及依据 ................................... 35 5.2.2算法结构 ......................................... 35 5.3 算法的Matlab实现...................................... 36
5.3.1非压缩带状对角方程组 ............................. 36 5.3.2压缩带状对角方程组 ............................... 38 5.4 实验结果及分析......................................... 40
5.4.1 Matlab运行结果 ................................. 40 5.4.2 Matlab运行结果分析 ............................. 43 5.5 实际问题............................................... 44 5.6 分析总结............................................... 47 6. 总结及致谢.................................................. 48
1. 上机题一
对以下和式计算：S??1?4211?????，要求： n?168n?18n?48n?58n?6??n?0? (1)若只需保留11个有效数字，该如何进行计算； (2)若要保留30个有效数字，则又将如何进行计算；
计算方法
任何一个β进制、具有t位有效数字的实数x总可以表示成x??(其中，dj是满足
d1d2d3dt???2??3?...??l，)??t1?d1??0?dj??的整数。
有效数字是指一个近似数的“有意义”的数字的数位，通常在十进制中讨论，设
~?10l~x??0.????l12t，若x??(0.?1?2??t?)?10?R,其中?i??0,1,??,9?1,?，0近似数
??x?x1~~?10l?t，则称x是x的具有t位有效数字的近似数，或称x具有t位有效数字。根据有2效数字的要求可以得到计算公式的精度要求。再利用后验误差估计式可以根据精度要求得到和式累加的项数。首先满足题目要求的情况下，利用后验误差估计计算出算法中的最大运行次数。即以一下条件作为迭代终止的准则：S(k?1)?S(k)??，其中?为误差数值。该实验第一问利用后验误差估计截取的项数，使误差小于1e-10；第二问利用后验误差估计截取的项数，使误差小于1e-30。得出运行次数，作为累加的终止条件。特别地，为了减小舍入误差在计算S时所采用的方法是逆序相加，其依据是：S中各项的绝对值为递减的正数，而两个数量级相差较大的数字相加减时，较小数的有效数字会丢失，从而导致最后的运算结果失真。为避免“大数吃小数”现象的发生，采用逆序相加；同时为了避免相近数相减和减小计算步骤，我们将
4211960n2?1208n?376???化简为进行计算。另外，要限定计算8n?18n?48n?5n8?6(8n?1)(8n?4)(8n?5)(8n?6)结果的有效数字的个数，在MATLAB软件中可直接调用函数vpa(变量名,位数)来控制计算结果的位数，最终完成精度要求较高的计算。
1.1 解题源程序
function solve1
clc; clear; sum1=0; sum2=0;
1
n=0; N=0;
%第一小问
%根据精度要求利用后验误差式，估计要使误差小于1?101?11?1?10?10，需要累加的项数n，
22所以要保证计算精度第n+1项要小于1?10?11，才能保证计算的精度
2while abs((960*n^2+1208*n+376)/((8*n+1)*(8*n+4)*(8*n+5)*(8*n+6))/(16.^n))>10^(-11)
%为了减小误差，我们对计算公式进行变形，一是防止相近数相减，二是为了减少计算步骤；
n=n+1; end n;
%累加进行的项数
disp('累加的次数:') n
%为减少舍入误差的影响，按绝对值递增的顺序求和。因此按n按从小到大的顺序相加,步距取为-1
for i=n:-1:0
sum1=sum1+((960*i.^2+1208*i+376)/((8*i+1)*(8*i+4)*(8*i+5)*(8*i+6)))/(16.^i); end
disp('保留11位有效数字时的近似值:') Sum1=vpa(sum1,11) %第二小问
%根据精度要求利用后验误差式，估计要使误差小于0.5e-30需要累加的项数N While
abs((960*N^2+1208*N+376)/((8*N+1)*(8*N+4)*(8*N+5)*(8*N+6))/(16.^N))>10^(-30) N=N+1; end N;
%累加进行的项数
disp('累加的次数:') N
%为减少舍入误差的影响，按绝对值递增的顺序求和。因此按N按从小到大的顺序相加,步距取为-1
for j=N:-1:0
sum2=sum2+((960*j.^2+1208*j+376)/((8*j+1)*(8*j+4)*(8*j+5)*(8*j+6)))/(16.^j); end
disp('保留30位有效数字时的近似值:') Sum2=vpa(sum2,30) End
2
if A(Rmax,k)==0 break; end
%交换当前行和主元所在行 for j=1:n
temp=A(k,j);
A(k,j)=A(Rmax,j); A(Rmax,j)=temp; end
temp=B(k); B(k)=B(Rmax); B(Rmax)=temp; for 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
B(i)=B(i)-A(i,k)*B(k); end end A
%--------回代过程---------------- x=zeros(1,n); x(n)=B(n)/A(n,n); for k=n-1:-1:1 S=B(k); for j=k+1:n
S=S-A(k,j)*x(j); end
x(k)=S/A(k,k); end
5.3.2压缩带状对角方程组 （1）%压缩格式Dat52.dat ； clear;
clc;
format short %读出数据
fid=fopen('C:\\Documents and Settings\\Administrator\\桌面\\最近学习计划\\计算方法\\上机\\2015上机题目\\Dat52.dat','r');
38
D=fread(fid,5,'long',0); n=D(3)%矩阵维数 p=D(5)%上带宽
q=D(4)%下带宽
A=fread(fid,[p+q+1,n],'float'); B=fread(fid,[n,1],'float'); C = A'; D = B';
%压缩格式Dat54.dat; clear;
clc;
format short %读出数据
fid=fopen('C:\\Documents and Settings\\Administrator\\桌面\\最近学习计划\\计算方法\\上机\\2015上机题目\\Dat54.dat','r');
D=fread(fid,5,'long',0); n=D(3)%矩阵维数
p=D(5)%上带宽
q=D(4)%下带宽
A=fread(fid,[p+q+1,n],'float'); B=fread(fid,[n,1],'float'); C = A'; D = B';
%采用求解压缩型的列主元Gauss消去法 x=GAUSSPP1(C,D,n,p,q) （2）针对压缩格式的带状矩阵，尤其是大型矩阵的列主元Gauss消去法子程序GAUSSPP1如下：
%针对压缩格式的带状方程组的列主元Gauss消去法 function x= GAUSSPP1(A, B, n, p, q) %-------消去过程--------- for k=1:n-1
pmax = k + p; if pmax > n pmax = n; end
for i=k+1:pmax
A(i,k-i+p+1) = A(i,k-i+p+1)/A(k,p+1); qmax = k+q;
39
if qmax > n qmax = n; end
for j=k+1:qmax
A(i,j-i+p+1) = A(i,j-i+p+1) - A(i,k-i+p+1) * A(k,j-k+p+1); end
B(i) = B(i) - B(k) * A(i,k-i+p+1); end end
%--------回代过程---------------- x(n) = B(n)/A(n,p+1); for k=n-1:-1:1 sum = 0; qmax = k+q; if qmax > n qmax = n; end
for j=k+1:qmax
sum = sum + A(k,j-k+p+1)*x(j); end
x(k) = (B(k) - sum)/A(k,p+1); end end
5.4 实验结果及分析 5.4.1 Matlab运行结果
1、Dat51.dat运行结果如下：
40
2、Dat53.dat运行结果如下：41
3、Dat52.dat运行结果如下：
4、Dat54.dat运行结果如下：
42
x(n)=GG(n,n+1)/GG(n,n); for i=n-1:-1:1 gx=0;
for j=i+1:n
gx=gx+GG(i,j)*x(j); end
x(i)=(GG(i,n+1)-gx)/GG(i,i); end
5、误差函数
%误差
function [rou] = wucha_wlg(GG) [m,n]=size(GG); n=n-1; rou=0;
for i=n+1:m
rou=rou+GG(i,n+1)^2; end %误差
3.5.1 QR分解输出结果
拟合二次函数如下及相应的系数 p2 =-0.0938 2.6064 8.3063
Px2 =(5869189440802495*x)/2251799813685248 - (1262*x^2)/13455 + 24296/2925 估算误差e2 =280.3395
拟合三次函数如下及相应的系数
p2 =-0.0084 0.2075 -0.2273 13.3880
Px3 =- (2412224987398535*x^3)/288230376151711744 (4095195769401595*x)/18014398509481984 + 7832/585
估算误差e3 =131.0618
+ (6142*x^2)/29601 -
拟合四次函数如下及相应的系数
p4 =0.0009 -0.0532 0.8909 -3.7050 16.7939 Px4 =(8622444092144931*x^4)/9223372036854775808 - (7672945562807963*x^3)/144115188075855872 + (8024210619243375*x^2)/9007199254740992 - (8342830140799029*x)/2251799813685248 + 21916/1305
估算误差e4 =59.0412 拟合五次函数如下及相应的系数
p5 =0.0001 -0.0045 0.0623 -0.1202 -0.4997 14.9504 Px5 =(3360357177204419*x^5)/36893488147419103232 - (2611432097870069*x^4)/576460752303423488 + (8975699045606965*x^3)/144115188075855872 - (4328904268112707*x^2)/36028797018963968 - (9001816270954031*x)/18014398509481984 + 4208151280463277/281474976710656
28
估算误差e5 =33.1446 Matlab输出图形如图所示
QR分解二阶函数拟合曲线3530252520201510515QR分解三阶函数拟合曲线35300510152025100510152025QR分解四阶函数拟合曲线3535QR分解五阶函数拟合曲线3030252520201515100510152025100510152025
3.5.1 结果分析
由结果可以看出，随着拟合阶数的增多误差有明显减少的趋势，随着拟合次数的增多正交分解法拟合曲线很接近，但是并不是阶数越多越好，阶数过大的话，局部会产生震荡。法方向最小二乘法随着次数的增多拟合误差减小，八阶至十五阶函数拟合较好，误差最小。在同样拟合次数的情况下，正交分解法比法方向最小二乘法的拟合误差要小。故有时在不知道函数形态的情况下，有时候如果选择正交多项式作为拟合的基函数效果往往会很不错。从拟合曲线可以大致看出一天中温度的变化规律是温度从零点时到大约三点是呈现持续降温趋势，之后到14点的过程中呈现持续升温趋势且14点时温度达到最高，之后到24点的过程中又呈现持续降温趋势，当日平均气温为21.2摄氏度。十六阶拟合时误差只有3.6941。
29
4. 上机题四
设计算法，求出非线性方程6x5?45x2?20?0的所有根，并使误差不超过10。
?44.1 二分法算法组织
4.1.1算法思想及依据
当函数f(x)连续时，设x?[a,b]，若f(a)f(b)?0，则方程f(x)?0在x?[a,b]至少有一个根。对此区间二分，找到有根区间，再二分，直到有根区间小到误差所允许的范围内，二分过程结束，近似根为此区间的中点。显然，每次迭代使区间长度减少一半，故二分法总是收敛的。
由f(x)?6x5?45x2?20?3x2(2x3?15)?20可知，当x?2时，f(x)必定永远大于0，当x??1时，f(x)必定永远小于0，所以6x5?45x2?20?0在???,?1???2,??的区间内必定没有根。而由
52f(?1)??31?0，f(0)?20?0，f(1)??19?0和f(2)?32?0可知，6x?45x?20?0在区间
52??1,0?，?0,1?和?1,2?上分别有一个根。并且利用Matlab校验非线性方程6x?45x?20?0在
x?[?1,2]上的个数。在每一个区间上分别采用二分法，即可求出此线性方程的三个根。
4.1.2算法结构
（1） 确定初始区间I(0)?[x(0),x(1)]，满足f(x(0))f(x(1))?0,k?0； （2） 计算区间I(k)?[x(k),x(k?1)]的中点x?1(x(k)?x(k?1))；
2（3） 若f(x(k))f(x)?0，则令x(k?1)?x，取I(k?1)?[xk(x；否则，令x(k)?x，取,)]k?(1)；I(k?1)?[x,x]
（4） 计算区间长度，若转（2）。
x(k)?x(k?1)??，则停止，令x*?1(x(k)?x(k?1))，输出；否则，令k?k?1，
24.2 二分法算法的Matlab实现
4.2.1利用Matlab校验根的区间
在[-1,2]上做出函数的图形，由图形可知函数的根所在的区间。之所以在求解函数零点之前，需要绘制函数图形，一方面是为了能够在后面的步骤中，在二分法、割线法中更好的选取初值；另一方面是为了校验算法的可靠性，从而由函数图像判断根的准确性。Matlab程序如下：
clear
x=[-1:0.01:2];
y=6.*x.^5-45.*x.^2+20;
30
%绘制函数图形
plot(x,y,'r','Linewidth',1.5) hold on
%添加水平线
h=line([-1,2],[0,0]);
%设置直线的宽度和颜色 set(h,'Linewidth',1.5) set(h,'color','k') %设置坐标轴刻度
set(gca,'Xtick',[-1:0.1:2])
%添加图形标题和坐标轴名称 title('The zero of function') grid
xlabel('x') ylabel('f(x)')
生成图形如下图所示，可见上述二分法所选区间符合条件。即分别取区间[-1,0],[0,1]和[1,2]为二分法的区间，分别在每个区间上应用二分法Matlab程序。即可得到函数f(x)?0在每个区间上所得到的根的数值解。
函数图像40302010f(x)0-10-20-30-40-1-0.9-0.8-0.7-0.6-0.5-0.4-0.3-0.2-0.100.10.20.30.40.5x0.60.70.80.911.11.21.31.41.51.61.71.81.92
4.2.2使用二分法Matlab程序 clc clear
31
f=inline('6*x^5-45*x^2+20'); e=1e-4;
%求解在区间[-1,0]上的根
[r1,e1]=erfen(f,e,-1,0); %调用二分法子程序 disp('在区间[-1,0]上的根为'); r1
disp('误差='); e1
%求解在区间[0,1]上的根 [r2,e2]=erfen(f,e,0,1);
disp('在区间[0,1]上的根为'); r2
disp('误差为'); e2
%求解在区间[1,2]上的根 [r3,e3]=erfen(f,e,1,2);
disp('在区间[1,2]上的根为'); r3
disp('误差为'); e3
编制二分法程序：erfen.m function [r,a]=erfen(f,e,x0,x1) Q=1;
while Q %设置循环变量标识 xm=(x0+x1)/2; %求取中点值 if f(xm)*f(x0)<0
x1=xm; %取中点值为右端值 else x0=xm; %否则取中点值为左端值 end
a=x1-x0;
if a r=x1; %输出非线性方程根的值 end 4.3 二分法实验结果及分析 方程6x5?45x2?20?0的三个根分别位于三个区间??1,0?，32 ?0,1?和?1,2?上，程序计算结果 为： 在区间[-1,0]上的根为：r1 = -0.6545；误差为 e1 = 6.1035e-005； 在区间[0,1]上的根为： r2 = 0.6812；误差为 e2 = 6.1035e-005； 在区间[1,2]上的根为： r3 = 1.8708；误差为 e3 = 6.1035e-005。 可以看出，通过用二分法对方程组根的求解可以发现这种方法总是收敛的，但是用这种方法的运算过程较慢而且运算量较大，因为他没有充分利用函数的性质而只是用到函数的符号，且每次只将区间缩小一半，所以其精度也不高，可以用其他更好的方法如牛顿法或割线法来求解。 下面采用割线法求解： 4.4 割线法算法组织 4.4.1算法思想及依据 在Newton法中，涉及函数求导，从而增加了计算量。可以用迭代前后两次的割线来替代曲线，把割线与x轴的交点作为根的近似，如此反复，这样就可以求解f(x)?0的计算方法，而不需要计算一次导数。 割线的方程为 y?f(xk)?f(xk)?f(xk?1)(x?xk) xk?xk?1f(xk)(xk?xk?1) f(xk)?f(xk?1)由该直线方程与x轴交点横坐标容易算出 x?xk?我们就把这个x作为下一次迭代的初值，那么就得到割线法的迭代公式 xk?1?xk?4.4.2算法结构 f(xk)(xk?xk?1) f(xk)?f(xk?1)（1）设定前两步的迭代值，利用迭代公式进行迭代。 （2）设置循环变量，当两次迭代的误差小于1.0e-4时，迭代终止，输出计算值及迭代次数。 4.5 割线法算法的Matlab实现 编制割线法程序：ger.m function y=ger(x0,x1) x2=x1-fc(x1)*(x1-x0)/(fc(x1)-fc(x0)); n=1; 33 while(abs(x1-x0)>1.0e-4)&(n<100000000) x0=x1;x1=x2； x2=x1-fc(x1)*(x1-x0)/(fc(x1)-fc(x0)); n=n+1; end x2 n end 编制函数文件：fc.m function y=fc(x) y=6*x^5-45*x^2+20; end 在Matlab命令窗口中分别计算 ger(-0.5,-0.4) ger(0.7,0.6) ger(1.9,1.8) 4.5.1割线法实验结果及分析 526x?45x?20?0的三个根在上述程序计算结果为： 方程 在区间[-1,0]上的根为：ger(-0.5,-0.4) x1 =-0.7050 x1 =-0.6406 x1 =-0.6539 x1 =-0.6546 x1 =-0.6545 x2 =-0.6545 n =6 在区间[0,1]上的根为：ger(0.7,0.6) x1 =0.6804 x1 =0.6812 x1 =0.6812 x2 =0.6812 n =4 在区间[1,2]上的根为：ger(1.9,1.8) x1 =1.8670 x1 =1.8713 x1 =1.8708 x1 =1.8708 x2 =1.8708 n =5 从计算结果可以看出迭代次数相对而言还是比较少的，总在可接受的范围之内。当然迭代次数与函数本身有关，同时和初始值的选取也是有关的。从计算结果来看，计算结果是有效的。 一般而言，割线法比牛顿法收敛速度要慢一些，但是割线法每一步只需计算一次函数值，而牛顿法则需计算导数值。因此，割线法的收敛速度也是相当快的。割线法具有超线性的收敛性，接近于平方收敛，故实际应用比较多。 5. 上机题五 编写程序实现大规模方程组的列主元高斯消去法程序，并对所附的方程组进行求解。具体的要求参见所附的相关文档。针对本专业中所碰到的实际问题，提炼一个使用方程组进行求解 34 的例子，并对求解过程进行分析、求解。 5.1 算法思想 列主元高斯消去法与LU分解的算法思想类似，但为了保证消去时，不会发生上溢，Gauss消去过程中，要适当交换方程的顺序以保证消去过程能够顺利进行并且计算解的精确满足要求。 ?k?1??k?1?消去过程中产生的数akk称为第k步消去的主元。交换方程的原则是使aik??i?k,k?1,...,n?中，绝对值最大的一个换到(k,k)位置而成为第k步消去的主元，带有这种交换的Gauss消去法即为列主元Gauss消去法。 程序要求对稀疏矩阵、条状对角矩阵以及压缩格式的带状矩阵进行列主元高斯消去求解，列主元高斯消去程序参考《计算方法教程》P31页的GAUSS算法。 5.2 算法组织 5.2.1算法思想及依据 消去法的中心是“降维”，即求解n元方程组的问题转化为先解n-1元方程组，一旦这个n-1元方程组的解取得，则剩余的一个未知量自然可以求得。这样逐步减少未知量个数的方法，是求解多元方程组的一个重要思想。 Gauss消去过程中，适当交换方程的顺序对保证消去过程能够顺利进行及计算解的精确度都是必要的。消去过程中产生的数akk(k?1)称为第k步消去的主元。交换方程的原则是使 ?aikk?1???i?k,k?1,...,n?中，绝对值最大的一个换到(k,k)位置而成为第k步消去的主元，带有这种交 换的Gauss消去法为列主元Gauss消去法。 5.2.2算法结构 一般地，列主元消去包含两个主要步骤：消去和回代，算法实现结构具体过程如下： 列主元Gauss消去法GAUSSPP(A,b,n,x) 1. For k=1,2,…,n-1 1.1找满足 ??kk?max?iki?k的下标 ?k ??01.2If ?kk then 输出错误信息；stop 1.3For j=1,2,…,n ????kj1.3.1 kj ????k1.4 k 1.5 For i=k+1,k+2,…,n 1.5.1 ?ik/?kk??ik 1.5.2 For j=k+1,k+2,…,n 35 ???ik?kj??ij1.5.2.1 ij 1.5.3 ?i??ik?k??i 1. [回代过程] 1.1 ?n/?nn?xn 1.2 For k=n-1,n-2,…,1 n?????x?k?kjj?/?kk?xkj?k+1? ? 5.3 算法的Matlab实现 5.3.1非压缩带状对角方程组 （1）%非压缩格式Dat51.dat; clear; clc; format short %读出数据 fid=fopen('C:\\Documents and Settings\\Administrator\\桌面\\最近学习计划\\计算方法\\上机\\2015上机题目\\Dat51.dat','r') D=fread(fid,3,'long',0); n=D(3) %矩阵维数 offset=20; fseek(fid,offset,'bof'); A=zeros(n); B=zeros(1,n); C=fread(fid,n*(n+1),'float',0); for i=1:n for j=1:n A(i,j)=C((i-1)*n+j); end end A; for i=1:n B(i)=C(n*n+i); end B; x=GAUSSPP(A,B) %非压缩格式Dat53.dat; 36 clear; clc; format short %读出数据 fid=fopen('C:\\Documents and Settings\\Administrator\\桌面\\最近学习计划\\计算方法\\上机\\2015上机题目\\Dat53.dat','r') D=fread(fid,3,'long',0); n=D(3) %矩阵维数 offset=20; fseek(fid,offset,'bof'); A=zeros(n); B=zeros(1,n); C=fread(fid,n*(n+1),'float',0); for i=1:n for j=1:n A(i,j)=C((i-1)*n+j); end end A; for i=1:n B(i)=C(n*n+i); end B; x=GAUSSPP(A,B) （2）列主元高斯消去法Matlab函数程序： %列主元高斯消去法 function x=GAUSSPP(A,B) n=length(B); %--------消去过程------------ %找出列主元 for k=1:n-1 MAX=abs(A(k,k)); Rmax=k; for i=k+1:n if abs(A(i,k))>MAX MAX=abs(A(i,k)); Rmax=i; end end 37