MATLAB计算方法3解线性方程组计算解法
MATLAB计算方法3解线性方程组计算解法
直到(n-1) 原方程组化为
a11 x1 a12 x2 a1n xn a1,n1 a22 x2 a2 n xn a2 ,n1
ann xn an ,n1
(上三角方程组) （3.2) 以上为消元过程。
(n) 回代求解公式
a n ,n1 xn a nn n x k 1 [a k ,n1 a kj x j ] a kk j k 1 ( k n 1, n 2,...,1)
由矩阵乘法 (1) 1) l11 a11 l11
umj 1 ukj a kj ukj a kj l km umj
m 1
k 1
2 求L的第k列：用L的第i行 u的第k列
(i k 1, , n)，即 ( l i 1 , , l ik , l kk , 0 0) ( u1k , u2 k , , ukk , 0 0)' a ik
( 2) 1)求u的第2行：用L的第2行 u的第j列 （j 2, , n) l 21 u1 j 1 u2 j a 2 j u2 j a 2 j l 21u1 j 2)求L的第2列：用L的第i行 u的第2列 （i 3,4, , n) l i 1 u12 l i 2 u22 a i 2 l i 2 (a i 2 l i 1 u12 ) / u22
m 1
l
k 1
im
umk l ik ukk a ik
k 1
l ik a ik l im umk ukk m 1
LU分解式： u1 j a1 j ( j 1,2, n) l i 1 a i 1 u11 ( i 2,3, , n) k 1 ukj a kj l km umj a kj m 1 ( j k , k 1, , n) k 1 l ik a ik l im umk ukk a ik m 1 ( i k 1, , n) ( k 2, 3, , n )
MATLAB计算方法3解线性方程组计算解法
MATLAB计算方法3解线性方程组计算解法线性方程组是数学中的一个重要问题，解线性方程组是计算数学中的一个基本计算，有着广泛的应用。
MATLAB是一种功能强大的数学软件，提供了多种解线性方程组的计算方法。
本文将介绍MATLAB中的三种解线性方程组的计算方法。
第一种方法是用MATLAB函数“linsolve”解线性方程组。
该函数使用高斯消元法和LU分解法求解线性方程组，可以处理单个方程组以及多个方程组的情况。
使用该函数的语法如下：X = linsolve(A, B)其中A是系数矩阵，B是常数向量，X是解向量。
该函数会根据A的形式自动选择求解方法，返回解向量X。
下面是一个使用“linsolve”函数解线性方程组的例子：A=[12;34];B=[5;6];X = linsolve(A, B);上述代码中，A是一个2×2的系数矩阵，B是一个2×1的常数向量，X是一个2×1的解向量。
运行代码后，X的值为[-4.0000;4.5000]。
第二种方法是用MATLAB函数“inv”求解逆矩阵来解线性方程组。
当系数矩阵A非奇异（可逆）时，可以使用逆矩阵求解线性方程组。
使用“inv”函数的语法如下：X = inv(A) * B其中A是系数矩阵，B是常数向量，X是解向量。
该方法先计算A的逆矩阵，然后将逆矩阵与B相乘得到解向量X。
下面是一个使用“inv”函数解线性方程组的例子：A=[12;34];B=[5;6];X = inv(A) * B;上述代码中，A是一个2×2的系数矩阵，B是一个2×1的常数向量，X是一个2×1的解向量。
运行代码后，X的值为[-4.0000;4.5000]。
第三种方法是用MATLAB函数“mldivide”（或“\”）求解线性方程组。
该函数使用最小二乘法求解非方阵的线性方程组。
使用“mldivide”函数的语法如下：X=A\B其中A是系数矩阵，B是常数向量，X是解向量。
A03 线性代数部分之MATLAB解方程
解方程组 矩阵的特征值与特征向量
郑丰华
成都信息工程学院数学学院 郑丰华
求解代数方程
一般代数方程包括线性（Linear）、非 一般代数方程包括线性（Linear）、非 线性（Nonlinear）和超越方程 线性（Nonlinear）和超越方程 （Transcedental equation）等。求解 equation）等。求解 的指令是solve。 的指令是solve。
成都信息工程学院数学学院 郑丰华
求矩阵的秩
rank rref 讨论线性方程组的解
– 求特解：linsolve(),pinv(A)*b 求特解：linsolve(),pinv(A)*b – 求通解：null() 求通解：null()
成都信息工程学院数学学院 郑丰华
特征值与特征向量
表3-3 eig命令 命令 d=eig(A) [V,D]=eig(A) 功能 求特征值及特征向量
成都信息工程学院数学学院 郑丰华
表3-1 solve命令 命令 S=solve('eq1','eq2',… S=solve('eq1','eq2',…,'eqn', 'v1','v2',… 'v1','v2',…,'vn') 功能 求方程或方程组关于 指定变量的解
说明： 1.'eq1','eq2',…,'eqn'或是字符串表达的方程，或是字 符串表达式；'v1','v2',…,'vn'是字符串表达的求解 变量名； 2.如果'eq1','eq2',…,'eqn'是不含有“等号”的表达式 ，则指令是对eq1=0,eq2=0,…,eqn=0的求解； 3.求解结果S是一个构架数组。如果要显示求解结果，必 须采用S.v1,S.v2,…,S.vn的援引方式； 4.指令solve在缺省规则下，还有一些形式更为简单的调 用方法，但推荐大家使用字符串格式； 5.在得不到“封闭型解析解”时，solve会给出数值解。
matlab线性方程组的矩阵解法
function x=lupdsv(A,b) n=length(b); [LU,p]=lupd(A); y(1)=b(p(1)); for i=2:n y(i)=b(p(i))-LU(i,1:i-1)*y(1:i-1)'; end x(n)=y(n)/LU(n,n); for i=(n-1):-1:1 x(i)=(y(i)-LU(i,i+1:n)*x(i+1:n)')/LU(i,i); end
lupdsv.m %功能：调用列主元三角分解函数 [LU,p]=lupd(A) % 求解线性方程组Ax=b。 。 求解线性方程组
%解法：PA=LU, Ax=b←→PAx=Pb 解法： 解法 % % LUx=Pb, Ly=f=Pb, y=Ux f(i)=b(p(i))
%输入：方阵A，右端项 （行或列向量均可） 输入：方阵 ，右端项b（行或列向量均可） 输入 %输出：解x（行向量） 输出： 输出 （行向量）
Ax = d 用矩阵表示 应用追赶法求解三对角线性方程组。追赶法仍然 追赶法求解三对角线性方程组 应用追赶法求解三对角线性方程组。追赶法仍然 保持LU分解特性,它是一种特殊的LU分解。 LU分解特性 LU分解 保持LU分解特性,它是一种特殊的LU分解。充分利用 了系数矩阵的特点，而且使之分解更简单, 了系数矩阵的特点，而且使之分解更简单,得到对三对 角线性方程组的快速解法。 角线性方程组的快速解法。
基于matlab平台的三种迭代法求解矩阵方程
数值分析第二次作业学院：电子工程学院基于matlab平台的三种迭代法求解矩阵方程组求解系数矩阵由16阶Hilbert方程组构成的线性方程组的解，其中右端项为[2877/851,3491/1431,816/409,2035/1187,2155/1423,538/395,1587/1279,573/502,947 /895,1669/1691,1589/1717,414/475,337/409,905/1158,1272/1711,173/244].要求：1）Gauss_Sedel迭代法；2）最速下降法；3）共轭梯度法；4）将结果进行分析对比。
解：根据题目要求，编写了对应算法的matlab程序，求解结果如下：（求解精度为10e-4，最大迭代次数1000）1、方程的解：如下图1所示图1 三种方法求解的结果对比图2 Gause_Sedel算法收敛特性图3 最速下降法收敛特性图3 共轭梯度法收敛特性从图中可以看到，在相同的最大迭代次数和预设求解精度条件下，共轭梯度算法仅需要4次迭代便可求出方程组的解，耗时0.000454秒，而且求出解的精度最高；Gauss_Sedel方法需要465次迭代，耗时0.006779秒，求解精度最差；最速下降法需要398次迭代，耗时0.007595秒，求解精度与共轭梯度算法差不多，因此两者求出的解也几乎相同。
从中可以得出结论，共轭梯度算法无论从求解精度还是求解速度上都优于其他两种，最速下降法在求解精度上几乎与共轭梯度算法持平，但求解速度更慢。
Gauss_Sedel方法在求解精度和速度两方面都最差。
具体的解为:Gauss_Sedel迭代法:(共需465次迭代，求解精度达到9.97e-5) X=[0.995328360833192 1.01431732497804 1.052861239300110.934006974137998 0.931493373808838 0.9665081384030661.00661848511341 1.03799789809258 1.051806903036541.06215849948572 1.04857676431223 1.028561990411131.01999170162638 0.971831831519515 0.9525261666348130.916996019179182].最速下降法:(共需398次迭代，求解精度达到9.94e-5)X=[0.998835379744322 1.01507463472900 0.9825890937201850.980191460759243 0.991245169713628 1.003780222253291.01350884374478 1.01928337905816 1.020859096651941.01930314197028 1.01444777381651 1.007040589892970.998384452250809 0.987399404644377 0.9757678149709120.963209150871750].共轭梯度法:(共需4次迭代，求解精度达到3.98e-5)X=[0.996472751179456 1.02707840189049 0.9776233734098530.973206695321590 0.986133032967607 1.001289025642341.01322158496914 1.02047386502293 1.023009050605651.02163015083975 1.01678089454399 1.009203108638740.999772406055155 0.988443827498859 0.9760941924969490.962844741655005].Matlab程序主程序:clc;clear;%% 本程序用于计算第二次数值分析作业，关于希尔伯特矩阵方程的解，用三种方法，分析并比较，也可推广至任意n维的矩阵方程%%A=hilb(16); %生成希尔伯特系数矩阵b=[2877/851;3491/1431;816/409;2035/1187;2155/1423;538/395;1587/1279;573/502;947/895;166 9/1691;1589/1717;414/475;337/409;905/1158;1272/1711;173/244]; %右端向量M=1000; %最大迭代次数err=1.0e-4; %求解精度[x,n,xx,cc,jingdu]=yakebi_diedai(A,b,err,M); % 雅克比算法求解tic;[x1,n1,xx1,cc1,jingdu1]=gauss_seidel(A,b,err,M); % gauss_seidel算法求解toc;tic;[x2,n2,xx2,jingdu2]=zuisuxiajiangfa(A,b,err,M); % 最速下降法求解toc;tic;[x3,flag,jingdu3,n3]=bicg(A,b,err); % matlab内置双共轭梯度算法求解toc;tic;[x4,xx4,n4,jingdu4]=con_grad(A,b,err,M); % 教材共轭梯度算法求解toc;%% 计算相应结果，用于作图%%num=[1:16]';jie=[num,x1,x2,x4]; % 三者的解对比% 三者的收敛情况对比num1=[1:n1]';fit1=[num1,jingdu1'];num2=[1:n2]';fit2=[num2,jingdu2'];num4=[1:n4]';fit4=[num4,jingdu4'];子函数1(Gause_Sedel算法):function [x,n,xx,cc,jingdu] = gauss_seidel(A,b,err,M)% 利用迭代方法求解矩阵方程这里是高斯赛尔得迭代方法% A 为系数矩阵b 为右端向量err为精度大小返回求解所得向量x及迭代次数% M 为最大迭代次数cc 迭代矩阵普半径jingdu 求解过程的精度n 所需迭代次数xx 存储求解过程中每次迭代产生的解for ii=1:length(b)if A(ii,ii)==0x='error';break;endendD=diag(diag(A));L=-tril(A,-1);U=-triu(A,1);B=(D-L)\U;cc=vrho(B); %迭代矩阵普半径FG=(D-L)\b;x0=zeros(length(b),1);x=B*x0+FG;k=0;xx(:,1)=x;while norm(A*x-b)>errx0=x;x=B*x0+FG;k=k+1;xx(:,k+1)=x;if k>=Mdisp('迭代次数太多可能不收敛！');break;endjingdu(k)=norm(A*x-b);endend子函数2(最速下降算法):function [x,n,xx,jingdu]=zuisuxiajiangfa(A,b,eps,M)% 利用迭代方法求解矩阵方程这里是最速下降迭代方法% A 为系数矩阵b 为右端向量err为精度大小返回求解所得向量x及迭代次数% % M 为最大迭代次数jingdu 求解过程的精度n 所需迭代次数xx 存储求解过程中每次迭代产生的解x0=zeros(length(b),1);r0=b-A*x0;t0=r0'*r0/(r0'*A*r0);x=x0+t0*r0;r=b-A*x;xx(:,1)=x;k=0;while norm(r)>epsr=r;x=x;t=r'*r/(r'*A*r);x=x+t*r;r=b-A*x;k=k+1;xx(:,k+1)=x;if k>=Mdisp('迭代次数太多可能不收敛！');break;endn=k;jingdu(k)=norm(r);endend子函31(共轭梯度法):function [x,xx,n,jingdu]=con_grad(A,b,eps,M)% 利用迭代方法求解矩阵方程这里是共轭梯度迭代方法% A 为系数矩阵b 为右端向量err为精度大小返回求解所得向量x及迭代次数% M 为最大迭代次数jingdu 求解过程的精度n 所需迭代次数xx 存储求解过程中每次迭代产生的解x0=zeros(length(b),1);r0=b-A*x0;p0=r0;% t0=r0'*r0/(r0'*A*r0);% x=x0+t0*r0;% xx(:,1)=x;k=0;x=x0;r=r0;p=p0;while norm(r)>epsx=x;r=r;p=p;afa=r'*r/(p'*A*p);x1=x+afa*p;r1=r-afa*A*p;beta=r1'*r1/(r'*r);p1=r1+beta*p;x=x1;r=r1;p=p1;k=k+1;xx(:,k)=x;if k>=Mdisp('迭代次数太多可能不收敛！');break;endn=k;jingdu(k)=norm(r);endend。
Matlab解方程(方程组)
Matlab 解方程这里系统的介绍一下关于使用Matlab求解方程的一系列问题，网络上关于Matlab求解方程的文章数不胜数，但是我大体浏览了一下，感觉很多文章都只是零散的介绍了一点，都只给出了一部分Matlab函数例子，以至于刚接触的人面对不同文章中的不同函数一脸茫然，都搞不清楚这些函数各自的用途，也不知道在什么样的情况下该选择哪个函数来求解方程，在使用Matlab解方程时会很纠结。
不知道读者是否有这样的感觉，反正我刚开始接触时就是这样的感觉，面对网络搜索到一系列函数都好想知道他们之间是个什么关系。
所谓的方程就是含有未知数的等式，解方程就是找出使得等式成立时的未知数的数值。
求方程的解可以转换成不同形式，比如求函数的零点、多项式的根。
方程分类很多，按照未知数个数分为一元、二元、多元方程；按照未知数组合形式分为线性方程和非线性方程；按照非零项次数是否一致分为齐次方程和非齐次方程。
线性方程就是方程中未知数次数是一次的，未知数之间不存在指、对、2及以上幂次的关系，线性方程又分为一元线性方程，也就是一元一次方程；多元线性方程，也就是多元一次方程，多以线性方程组的形式出现（包括齐次线性方程组和非齐次线性方程组）。
在Matlab中求解方程的函数主要有roots、solve、fzero、和fsolve函数等，接下来详细的介绍一下各个Matlab函数的使用方法和使用场合。
一、直接求解法(线性方程组)直接求解法不需要借助任何的Matlab函数，主要用于求解线性方程组，也就是未知数次数是一次的方程组，包括齐次线性方程组合非齐次线性方程组。
当然既然可以求解方程组自然也就可以求解单个方程。
主要针对A x=b形式的方程，其中A是未知数系数矩阵，x是未知数列向量，b是常数列向量，当b=0时就是齐次线性方程组，b ≠0时是非齐次线性方程组。
用左除法，x=A\b例：求解线性方程组的解12341242341234251357926640x x x x x x x x x x x x x x +-+=⎧⎪-+=-⎪⎨+-=⎪⎪+--=⎩解：即直接利用b 左除以A 。
matlab-解方程
1、解方程组问1：如何用matlab解方程组？这个问题其实在matlab中解方程组还是很方便的，例如，对于代数方程组Ax=b(A为系数矩阵，非奇异)的求解，MATLAB中有两种方法：(1)、x=inv(A)*b —采用求逆运算解方程组；(2)、x=A\ —采用左除运算解方程组。
例:x1+2x2=82x1+3x2=13>>A=[1,2;2,3];b=[8;13];>>x=inv(A)*bx =2.003.00>>x=A\bx =2.003.00；即二元一次方程组的解x1和x2分别是2和3。
问2：如何用matlab解多次的方程组？有符号解法，方法是：先解出符号解，然后用vpa(F,n)求出n位有效数字的数值解。
具体步骤如下：第一步:定义变量syms x y z ...;第二步:求解[x,y,z,...]=solve('eqn1','eqn2',...,'eqnN','var1','var2',...'varN');第三步:求出n位有效数字的数值解x=vpa(x,n);y=vpa(y,n);z=vpa(z,n);...。
如：解二(多)元二(高)次方程组：x^2+3*y+1=0y^2+4*x+1=0解法如下：>>syms x y;>>[x,y]=solve('x^2+3*y+1=0','y^2+4*x+1=0');>>x=vpa(x,4);>>y=vpa(y,4);结果是：x =1.635+3.029*i1.635-3.029*i-.283-2.987y =1.834-3.301*i1.834+3.301*i-.3600-3.307。
二元二次方程组，共4个实数根；问3，如何用matlab解高次方程组(非符号方程组)？举个例子好吗？解答如下：基本方法是：solve(s1,s2,…,sn,v1,v2,…,vn)，即求表达式s1,s2,…,sn组成的方程组，求解变量分别v1,v2,…,vn。
Matlab数组运算及线型方程组的求解
数组运算及线型方程组的求解1．“：”号的用法。
用“：”号生成行向量a=[1 2 3 4 5 6 7 89 10]、b=[5 3 1 -1 -3 -5]；用线性等分命令linspace重新生成上述的a和b向量。
另，在100和10000之间用对数等分命令logspace生成10维的向量c。
linspace(1,10,10) linspace(5,-5,6) ak=logspace(2,4,10)2. 已知多项式a(x)=x2+2x+3，b(x)=4x2+5x+6，求a，b积的微分。
a=[1 2 3];b=[4 5 6];c=polyder(a,b)c=poly2str(c,'x')3.用生成下列矩阵，取出方框内的数组元素a=[1:5;10:-1:6;11:15;16:20;21:25]q=a(2,2:3)m=a(2:4,4)n=a(4:5,1:3)4. 生成一个9×9维的魔方矩阵，提取其中心的3×3维子矩阵M，利用sum函数检验其各行和各列的和是否相等。
并且实现上述中心矩阵左旋90°或右旋90°，左右翻转，上下翻转a=magic(9)a =47 58 69 80 1 12 23 34 45 57 68 79 9 11 22 33 44 46 67 78 8 10 21 32 43 54 56 77 7 18 20 31 42 53 55 66 6 17 19 30 41 52 63 65 76 16 27 29 40 51 62 64 75 5 26 28 39 50 61 72 74 4 1536 38 49 60 71 73 3 14 2537 48 59 70 81 2 13 24 35>> m=a(4:6,4:6)m =20 31 42 30 41 52 40 51 62 >> c=rot90(m)c =42 52 62 31 41 51 20 30 40 >> c=rot90(m,-1)c =40 30 20 51 41 31 62 52 42 >> s=fliplr(m)s =42 31 20 52 41 30 62 51 40 >> w=flipud(m)w =40 51 62 30 41 52 20 31 425.已知a=[1 2 3:4 5 6;7 8 0]，求其特征多项式并求其根、特征值和特征多项式6. 计算二重不定积分7.求解微分方程。
- 1、下载文档前请自行甄别文档内容的完整性，平台不提供额外的编辑、内容补充、找答案等附加服务。
- 2、"仅部分预览"的文档,不可在线预览部分如存在完整性等问题,可反馈申请退款(可完整预览的文档不适用该条件!)。
- 3、如文档侵犯您的权益，请联系客服反馈,我们会尽快为您处理(人工客服工作时间：9:00-18:30)。
(3.1)
ai1x1 ai2 x2 ain xn ai,n1
an1x1 an2 x2 ann xn an,n1
消 元： 用a11将ai1(i 2, , n)化为零；
把
ai1 a11
第1行，加到第i
行。
--Jordan消去法)，但比有回代的列主元消
去法的乘除运算次数多。
(2)有回代的列主元消去法所进行的乘除运算
次数为 1 n3 n2 1 n，量很小。
3
3
Gauss 列主元消去法：
优点 ------ 计算结果更可靠；
缺点 ------ 挑主元花机时更多, 次序
有变动，程x序1 ,复,杂x。n
.
用Matlab实现选列主元Gauss消去法解线性方程组
问 题：a11 0或 a11 0?以后各步类似。
.
用Matlab实现顺序Gauss消去法 在Matlab程序编辑器中输入：
function x=nagauss(a,b,flag) %解线形方程组ax=b，
a为系数矩阵，b为右端列向量，flag若为0，则显示中间
过程，否则不显示，默认为0，x为解向量
if flag==0,a,end
end
.
% 回代 x=zeros(n,1); x(n)=a(n,n+1)/a(n,n); for k=n-1:-1:1
x(k,:)=(a(k,n+1)a(k,(k+1):n)*x((k+1):n))/a(k,k); end
.
.
程序运行 结果：
.
二 列主元素Gauss消去法---计算结果可靠
n
akj x j ]
j k 1
(k n 1, n 2,...,1)
（3.3)
.
系数矩阵为对称正定阵或严格对角占优 阵的方程组按高斯消去法计算是数值稳 定的，因而不必选主元。
严格对角占优阵：至少有一个主对角线 元素的绝对值严格大于此行或此列其他 元素的绝对值之和。
.
说明：
(1)也可采用无回代的列主元消去法(叫Gauss-
if nargin<3,flag=0;end
n=length(b); a=[a,b];
% 消元
for k=1:(n-1)
a((k+1):n,(k+1):(n+1))=a((k+1):n,(k+1):(n+1))-
a((k+1):n,k)/a(k,k)*a(k,(k+1):(n+1));
a((k+1):n,k)=zeros(n-k,1);
(1)找行号r1
使 ar11
max
1 i n
ai1 ，对调1
r1行：
消元：用a11消ai1为0 :
第1行
ai1 a11
加到第i行, 第i行第j个元素成为
aij
a1 j
ai1 a11
aij
（i 2,3, , n；j 1,2, , n 1）
.
到此原方程组化为
a11 x1 a12 x2 a1n xn a1,n1 a22 x2 a2n xn a2,n1
第三章线性方程组数值解法
解线性方程组 Ax b
直接法：理论，无舍入误差，有限步，精确解 迭代法：格式，无穷序列 解向量 x
§3.1 直接法
.
一、 Gauss 消去法
a11x1 a12 x2 a1n xn a1,n1 设 有 a21x1 a22 x2 a2n xn a2,n1
.
直到(n-1) 原方程组化为
a11 x1 a12 x2 a1n xn a1,n1 a22 x2 a2n xn a2,n1
ann xn an,n1
(上三角方程组) （ 3.2) 以上为消元过程。
.
(n) 回代求解公式
xn
an,n1
ann
xk
1 akk
[ak ,n1
为Crout分解。
.
为什么要讨论三角分解？若在消元法进行前能实
现三角分解 A LU ， 则
Ax b (LU)x b
Ly b (下三角方程组） Ux y (上三角方程组）
L ( LK L1 )1 ，则
A LU （下三角 上三角）
（三角因子分解）
.
Hale Waihona Puke 定义3.1 A LU 叫 A 的三角（因子）分解，其中 L是
下三角, U是上三角。
定义3.2 若 L为单位下三角阵（对角元全为1），
U 为上三角阵，则称 A LU 为Doolittle分解;
若L是L下三角，U 是单位上三角，则称 A LU
在Matlab程序编辑器中输入：
function x=nagauss2(a,b,flag) %a为系数矩阵；b为右 端列向量；flag若为0，则显示中间过程，否则不显示 if nargin<3,flag=0;end n=length(b); a=[a,b]; % 选主元 for k=1:(n-1) [ap,p]=max(abs(a(k:n,k)));p=p+k-1; if p>k,t=a(k,:); a(k,:)=a(p,:); a(p,:)=t; end % 消元
x(k,:)=(a(k,n+1)a(k,(k+1):n)*x((k+1):n))/a(k,k); end
.
.
程序运行 结果:
.
三．矩阵三角分解法
Gauss消元，初等行变换，化原方程组为上三 角型。
( A | b) (u | g) (| 1) Lk Lk1 L2 L1 ( A | b) (u | g) LK LK 1 L2 L1 A u A (LK LK 1 L2 L1 )1 u
ai 2 x2 ain xn ai,n1
an2 x2 ann xn an,n1
.
(2) 找r2，使 ar2 2
max
2 i n
ai2 ，
对调2 r2行.
消
元
：
用a22把ai
消
2
为0
(i
3,4,
, n) :
第2行
ai2 a22
第i行，则
aij
a2
j
ai2 a22
aij
（i 3,4， , n；j 2,3， , n 1）
a((k+1):n,(k+1):(n+1))=a((k+1):n,(k+1):(n+1))a((k+1):n,k)/a(k,k)*a(k,(k+1):(n+1));
a((k+1):n,k)=zeros(n-k,1);
.
if flag==0,a,end end %回代 x=zeros(n,1); x(n)=a(n,n+1)/a(n,n); for k=n-1:-1:1