线性方程组求解Matlab程序(精.选)
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是解向量。
求线性方程组AX=b通解的Matlab实现程序
［ ｎ １ 】 ［ ２ ｎ ２ 】
的通解加 上其 自身 的一个 特解 。在理 论基础 上 ，我们 利用下 面 的例子说 明Ｍａ ｔ １ ａ ｂ 实现 程序 。
解 ：Ｍａ ｔ ｌ ａ ｂ 实现程序如下：
＞ ＞ Ａ ＝ ［ １ １ ０— ３ － １ ； １— １ ２— １
一
７ ］ ；
＞ ＞ ｆ ｏ ｒ ｍａ ｔ ｒ ａ ｔ
＝
下面是 其简化形 式
［一 ｎ ｌ ＋７ ｎ ２ ］
［ｎ １ ＋５ ｎ ２ ］
Ｃ１ １ ＋Ｃ２ （ ｚ ２＋… ＋Ｃ
一
，
０ 【
一
， ，
在这里 Ｃ ， Ｃ ２ ， …， Ｃ ～ 设为任 意 的常数 。 利用Ｍａ ｔ ｌ ａ ｂ ，我 们要 求零 空 间 ，可 以调用 函数ｎ ｕ ｌ ｌ ，从
％限定 输 出格式 为有理 式
Ｉ ２ 十 Ｘ ２ 一 ｘ ３ ＋ ｘ ４ ＝ １ Ｉ ３ 一 ２ ＋ ２ 一 ３ ＝ ２
＞ ＞Ｐ＝ｎ ｕ ｌ ｌ （ Ａ ， ＇ ｒ ’ ） ％ 求线性方程组的解空间的有理形式
的基
＞ ＞ ｓ ｙ ｍｓ ｎ ｌ ，ｎ ２
例 ２ 求 方 程 组 ： １ ５ ＋ 一 Ｘ ３ ＋ ２ Ｘ ４ ： 一 １ 的 通 解。
关键 词 ：齐次线性 方程 组 ；非齐次线性 方程组 ；通解
引 言
求 解 线性 方 程 组 的 问题 是 数 值线 性 代 数 的三 大 问 题之
数值分析(hilbert矩阵)病态线性方程组的求解matlab程序
（Hilbert 矩阵）病态线性方程组的求解理论分析表明，数值求解病态线性方程组很困难。
考虑求解如下的线性方程组的求解Hx = b ，期中H 是Hilbert 矩阵，()ij n n Hh ，11ij h i j ，i ，j = 1,2，…，n 1.估计矩阵的2条件数和阶数的关系2.对不同的n ，取(1,1,,1)nx K ?，分别用Gauss 消去，Jacobi 迭代，Gauss-seidel 迭代，SOR 迭代和共轭梯度法求解，比较结果。
3.结合计算结果，试讨论病态线性方程组的求解。
第1小题：condition.m %第1小题程序t1=20;%阶数n=20x1=1:t1;y1=1:t1;for i=1:t1H=hilb(i);y1(i)=log(cond(H));endplot(x1,y1);xlabel('阶数n');ylabel('2-条件数的对数（log(cond(H)）');title('2-条件数的对数（log(cond(H)）与阶数n 的关系图');t2=200;%阶数n=200x2=1:t2;y2=1:t2;for i=1:t2H=hilb(i);y2(i)=log(cond(H));endplot(x2,y2);xlabel('阶数n');ylabel('2-条件数的对数（log(cond(H)）');title('2-条件数的对数（log(cond(H)）与阶数n 的关系图');画出Hilbert 矩阵2-条件数的对数和阶数的关系n=200时n=20时从图中可以看出，1）在n小于等于13之前，图像近似直线log（cond（H））~1.519n-1.8332）在n大于13之后，图像趋于平缓，并在一定范围内上下波动，同时随着n的增加稍有上升的趋势第2小题：solve.m%m第2小题主程序N=4000;xGauss=zeros(N,1);xJacobi=zeros(N,1);xnJ=zeros(N,1);xGS=zeros(N,1);xnGS=zeros(N,1);xSOR=zeros(N,1);xnSOR=zeros(N,1);xCG=zeros(N,1);xnCG=zeros(N,1);for n=1:N;x=ones(n,1);t=1.1;%初始值偏差x0=t*x;%迭代初始值e=1.0e-8;%给定的误差A=hilb(n);b=A*x;max=100000000000;%可能最大的迭代次数w=0.5;%SOR迭代的松弛因子G=Gauss(A,b);[J,nJ]=Jacobi(A,b,x0,e,max);[GS,nGS]=G_S(A,b,x0,e,max);[S_R,nS_R]=SOR(A,b,x0,e,max,w);[C_G,nC_G]=CG(A,b,x0,e,max);normG=norm(G'-x);xGauss(n)=normG;normJ=norm(J-x);nJ;xJacobi(n)=normJ;xnJ(n)=nJ;normGS=norm(GS-x);nGS;xGS(n)=normGS;xnGS(n)=nGS;normS_R=norm(S_R-x);nS_R;xSOR(n)=normS_R;xnSOR(n)=nS_R;normC_G=norm(C_G-x);nC_G;xCG(n)=normC_G;xnCG(n)=nC_G;endGauss.m%Gauss消去法function x=Gauss(A,b)n=length(b);l=zeros(n,n);x=zeros(1,n);%消去过程for i=1:n-1for j=i+1:nl(j,i)=A(j,i)/A(i,i);for k=i:nA(j,k)=A(j,k)-l(j,i)*A(i,k);endb(j)=b(j)-l(j,i)*b(i);endend%回代过程x(n)=b(n)/A(n,n);for i=n-1:-1:1c=A(i,:).*x;x(i)=(b(i)-sum(c(i+1:n)))/A(i,i);endJacobi.m%Jacobi迭代,x0表示迭代初值,e表示允许误差(迭代停止条件),n表示迭代次数,m 可能最大的迭代次数function [x,n]=Jacobi(A,b,x0,e,m)n=length(A);D=diag(diag(A));U=-triu(A,1);L=-tril(A,-1);B=D\(L+U);f=D\b;x=B*x0+f;n=1;while norm(x-x0)>ex0=x;x=B*x0+f;n=n+1;if n>mdisp('Jacobi迭代次数过多，迭代可能不收敛');break;endendG_S.m%Gauss-Seidel迭代,x0表示迭代初值,e表示允许误差(迭代停止条件),n表示迭代次数,m可能最大的迭代次数function [x,n]=G_S(A,b,x0,e,m)n=length(A);D=diag(diag(A));U=-triu(A,1);L=-tril(A,-1);B=(D-L)\U;f=(D-L)\b;x=B*x0+f;n=1;while norm(x-x0)>ex0=x;x=B*x0+f;n=n+1;if n>mdisp('Gauss-Seidel迭代次数过多，迭代可能不收敛');break;endendSOR.m%SOR超松弛迭代,x0表示迭代初值,e表示允许误差(迭代停止条件),n表示迭代次数,m可能最大的迭代次数,w松弛因子function [x,n]=SOR(A,b,x0,e,m,w)n=length(A);D=diag(diag(A));U=-triu(A,1);L=-tril(A,-1);B=(D-w*L)\((1-w)*D+w*U);f=(D-w*L)\b*w;x=B*x0+f;n=1;while norm(x-x0)>ex0=x;x=B*x0+f;n=n+1;if n>mdisp('SOR超松弛迭代次数过多，迭代可能不收敛');break;endendCG.m%CG共轭梯度法,x0表示迭代初值,e表示允许误差(迭代停止条件),n表示迭代次数,m可能最大的迭代次数function [x,n]=CG(A,b,x0,e,m)r=b-A*x0;p=r;alpha=(r'*r)/(p'*(A*p));x=x0+alpha*p;r1=b-A*x;n=1;while norm(r1)>ebelta=(r1'*r1)/(r'*r);p=r1+belta*p;r=r1;x0=x;alpha=(r'*r)/(p'*(A*p));x=x0+alpha*p;r1=b-A*x;n=n+1;if n>mdisp('CG共轭梯度法迭代次数过多，迭代可能不收敛');break;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．“：”号的用法。
用“：”号生成行向量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.求解微分方程。
matlab矩阵分解与线性方程组求解
格式
[Q, R] = rsf2csf(q, r) 例4-7
A=[1 1 1 3;1 2 1 1;1 1 3 1;-2 1 1 4]; [q, r]=schur (A) [Q, R]=rsf2csf(q, r)
4.2 秩与线性相关性
4.2.1
汪远征
矩阵和向量组的秩与向量组的线性相关性
矩阵 A 的秩是指矩阵 A 中最高阶非零子式的阶数，或
是矩阵线性无关的行数与列数；向量组的秩通常由该
向量组构成的矩阵来计算。 k = rank(A) 返回矩阵A的行(或列)向量中线性无关个数 k = rank(A,tol) tol为给定误差
在 MATLAB 中，求矩阵秩的函数是 rank 。其格式为：
4.2 秩与线性相关性
4.2.1
汪远征
矩阵和向量组的秩与向量组的线性相关性
4.2 秩与线性相关性
4.2.2
汪远征
求行阶梯矩阵及向量组的基
Matlab 将矩阵化成行最简形的命令是 rref或 rrefmovie 。
其格式为：
R = rref(A) R 是A的行最简行矩阵 [R,jb] = rref(A) jb 是一个向量，其含义为： r = length(jb) 为 A 的秩； A(:, jb)为A的列向量基；jb中元素表示基向量所在的 列。
阵。
4.1 矩阵分解
4.1.2
汪远征
Cholesky分解
例4-2
A=pascal(4) %产生4阶pascal矩阵 [R,p]=chol(A)
4.1 矩阵分解
4.1.3
汪远征QBiblioteka 分解将矩阵A分解成一个正交矩阵Q与一个上三角矩阵R的
乘积A=QR，称为QR分解。
MATLAB-平方根法和改进平方根法求解线性方程组例题与程序演示教学
M A T L A B-平方根法和改进平方根法求解线性方程组例题与程序（2）设对称正定阵系数阵线方程组12345678424024000221213206411418356200216143323218122410394334411142202531011421500633421945x x x x x x x x -⎡⎤⎡⎤⎡⎤⎢⎥⎢⎥⎢⎥---⎢⎥⎢⎥⎢⎥⎢⎥⎢⎥⎢⎥----⎢⎥⎢⎥⎢⎥----⎢⎥⎢⎥⎢⎥=⎢⎥⎢⎥⎢⎥----⎢⎥⎢⎥⎢⎥----⎢⎥⎢⎥⎢⎢⎥⎢⎥⎢---⎢⎥⎢⎥⎢--⎢⎥⎢⎢⎥⎣⎦⎣⎦⎣⎦⎥⎥⎥⎥ (1,1,0,2,1,1,0,2)T x *=--二、数学原理 1、平方根法解n 阶线性方程组Ax=b 的choleskly 方法也叫做平方根法，这里对系数矩阵A 是有要求的，需要A 是对称正定矩阵，根据数值分析的相关理论，如果A 对称正定，那么系数矩阵就可以被分解为的T A=L L •形式，其中L 是下三角矩阵，将其代入Ax=b 中，可得：T LL x=b 进行如下分解：T L xL by y ⎧=⎨=⎩ 那么就可先计算y,再计算x ，由于L 是下三角矩阵，是T L 上三角矩阵，这样的计算比直接使用A 计算简便，同时你应该也发现了工作量就转移到了矩阵的分解上面，那么对于对称正定矩阵A 进行Cholesky 分解，我再描述一下过程吧： 如果你对原理很清楚那么这一段可以直接跳过的。
设T A=L L •，即1112111112112122221222221212....................................n n n n n n nn n n nn nn a a a l l l l aa a l l l l a a a l l l l ⎡⎤⎡⎤⎡⎤⎢⎥⎢⎥⎢⎥⎢⎥⎢⎥⎢⎥=⎢⎥⎢⎥⎢⎥⎢⎥⎢⎥⎢⎥⎣⎦⎣⎦⎣⎦其中,,1,2,...,ij ji a a i j n ==第1步，由矩阵乘法，211111111,i i a l a l l ==g 故求得111111,2,3,...i i a l l i n a === 一般的，设矩阵L 的前k-1列元素已经求出 第k 步，由矩阵乘法得112211k k kk kmkkik im km ik kkm m a l l a l l l l --===+=+∑∑, 于是11(2,3,...,n)1(),1,2,...kk k ik ik im km m kk l k l a l l i k k n l -=⎧=⎪⎪=⎨⎪=-=++⎪⎩∑ 2、改进平方根法在平方根的基础上，为了避免开方运算,所以用TLDL A =计算；其中，11122.........n d D D D d ⎤⎤⎡⎤⎥⎥⎢⎥⎢⎥⎢⎥⎢⎥⎢⎥⎢⎥⎢⎥===⎢⎥⎢⎥⎢⎥⎢⎥⎢⎥⎢⎥⎢⎢⎢⎥⎣⎦⎣⎣；得1121212212111111n n n n n d l l l d l A l l d ⎡⎤⎡⎤⎡⎤⎢⎥⎢⎥⎢⎥⎢⎥⎢⎥⎢⎥=⎢⎥⎢⎥⎢⎥⎢⎥⎢⎥⎢⎥⎣⎦⎣⎦⎣⎦L L M MO O O M L按行计算的L 元素及D 对元素公式 对于n i ,,2,1Λ=11(1,21)j ij ij ik jk k t a t l j i -==-=-∑…，./(1,2,)ij ij j l t d j ==…，i-1.11i i ii ik ikk d a t l -==-∑计算出LD T =的第i 行元素(1,2,i-1)ij t j =…，后，存放在A 的第i 行相置，然后再计算L 的第i 行元素，存放在A 的第i 行.D 的对角元素存放在A 的相应位置.对称正定矩阵A 按T LDL 分解和按T LL 分解计算量差不多，但T LDL 分解不需要开放计算。
利用matlab解线性方程组
数值计算实验——解线性方程组西南交通大学2012级茅7班20123257 陈鼎摘要本报告主要介绍了基于求解线性方程组的高斯消元法和列主消元法两种数值分析方法的算法原理及实现方法。
运用matlab数学软件辅助求解。
实验内容1．编写用高斯消元法解线性方程组的MATLAB程序，并求解下面的线性方程组，然后用逆矩阵解方程组的方法验证。
2．编写用列主消元法解线性方程组的MATLAB程序，并求解下面的线性方程组，然后用逆矩阵解方程组的方法验证。
给定方程组如下：①0.325x1+2.564x2+3.888x3+5x4=1.521②-1.548x1+3.648x2+4.214x3-4.214x4=2.614③-2.154x1+1.647x2+5.364x3+x4=3.978④0x1+2.141x2-2.354x3-2x4=4.214A.高斯消元法一、算法介绍高斯消元法是一种规则化的加减消元法。
基本思想是通过逐次消元计算把需要求解的线性方程组转化成为上三角方程组，即把现形方程组的系数矩阵转化为上三角矩阵，从而使一般线性方程组的求解转化为等价的上三角方程组的求解。
二、matlab程序function [RA,RB,n,X]=gaus(A,b)B=[A b]; n=length(b); RA=rank(A);RB=rank(B);zhica=RB-RA;if zhica>0,disp(‘因为RA~=RB，所以此方程组无解.')returnendif RA==RBif RA==ndisp(‘因为RA=RB=n，所以此方程组有唯一解.')X=zeros(n,1); C=zeros(1,n+1);for p= 1:n-1for k=p+1:nm= B(k,p)/ B(p,p); B(k,p:n+1)= B(k,p:n+1)-m* B(p,p:n+1);endendb=B(1:n,n+1);A=B(1:n,1:n); X(n)=b(n)/A(n,n);for q=n-1:-1:1X(q)=(b(q)-sum(A(q,q+1:n)*X(q+1:n)))/A(q,q);endelsedisp(‘因为RA=RB<n，所以此方程组有无穷多解.')endend三、实验过程与结果输入的量：系数矩阵A和常系数向量b；输出的量：系数矩阵A和增广矩阵B的秩RA、RB，方程中未知量的个数n和有关方程组解X及其解的信息。
- 1、下载文档前请自行甄别文档内容的完整性，平台不提供额外的编辑、内容补充、找答案等附加服务。
- 2、"仅部分预览"的文档,不可在线预览部分如存在完整性等问题,可反馈申请退款(可完整预览的文档不适用该条件!)。
- 3、如文档侵犯您的权益，请联系客服反馈,我们会尽快为您处理(人工客服工作时间：9:00-18:30)。
线性方程组求解1.直接法Gauss消元法：function x=DelGauss(a,b)% Gauss消去法[n,m]=size(a);nb=length(b);det=1;%存储行列式值x=zeros(n,1);for k=1:n-1for i=k+1:nif a(k,k)==0returnendm=a(i,k)/a(k,k);for j=k+1:na(i,j)=a(i,j)-m*a(k,j);endb(i)=b(i)-m*b(k);enddet=det*a(k,k);enddet=det*a(n,n);for k=n:-1:1 %回代for j=k+1:nb(k)=b(k)-a(k,j)*x(j);endx(k)=b(k)/a(k,k);endExample:>> A=[1.0170 -0.0092 0.0095;-0.0092 0.9903 0.0136;0.0095 0.0136 0.9898]; >> b=[1 0 1]';>> x=DelGauss(A,b)x =0.9739-0.00471.0010列主元Gauss消去法：function x=detGauss(a,b)% Gauss列主元消去法[n,m]=size(a);nb=length(b);det=1;%存储行列式值x=zeros(n,1);for k=1:n-1amax=0;% 选主元for i=k:nif abs(a(i,k))>amaxamax=abs(a(i,k));r=i;endendif amax<1e-10return;endif r>k %交换两行for j=k:nz=a(k,j);a(k,j)=a(r,j);a(r,j)=z;endz=b(k);b(k)=b(r);b(r)=z;det=-det;endfor i=k+1:n %进行消元m=a(i,k)/a(k,k);for j=k+1:na(i,j)=a(i,j)-m*a(k,j);endb(i)=b(i)-m*b(k);enddet=det*a(k,k);enddet=det*a(n,n);for k=n:-1:1 %回代for j=k+1:nb(k)=b(k)-a(k,j)*x(j);endx(k)=b(k)/a(k,k);endExample:>> x=detGauss(A,b)x =0.9739-0.00471.0010Gauss-Jordan消去法：function x=GaussJacobi(a,b)% Gauss-Jacobi消去法[n,m]=size(a);nb=length(b);x=zeros(n,1);for k=1:namax=0;% 选主元for i=k:nif abs(a(i,k))>amaxamax=abs(a(i,k));r=i;endendif amax<1e-10return;endif r>k %交换两行for j=k:nz=a(k,j);a(k,j)=a(r,j);a(r,j)=z;endz=b(k);b(k)=b(r);b(r)=z;end%进行消元b(k)=b(k)/a(k,k);for j=k+1:na(k,j)=a(k,j)/a(k,k);endfor i=1:nif i~=kfor j=k+1:na(i,j)=a(i,j)-a(i,k)*a(k,j);endb(i)=b(i)-a(i,k)*b(k);endendendfor i=1:nx(i)=b(i);endExample:>> x=GaussJacobi(A,b) x =0.9739-0.00471.0010LU分解法：function [l,u]=lu(a)%LU分解n=length(a);l=eye(n);u=zeros(n);for i=1:nu(1,i)=a(1,i);endfor i=2:nl(i,1)=a(i,1)/u(1,1);endfor r=2:n%%%%for i=r:nuu=0;for k=1:r-1uu=uu+l(r,k)*u(k,i);endu(r,i)=a(r,i)-uu;end%%%%for i=r+1:nll=0;for k=1:r-1ll=ll+l(i,k)*u(k,r);endl(i,r)=(a(i,r)-ll)/u(r,r);end%%%%Endfunction x=lusolv(a,b)%LU分解求解线性方程组aX=b if length(a)~=length(b)error('Error in inputing!')return;endn=length(a);[l,u]=lu(a);y(1)=b(1);for i=2:nz=0;for k=1:i-1z=z+l(i,k)*y(k);endy(i)=b(i)-z;endx(n)=y(n)/u(n,n);for i=n-1:-1:1z=0;for k=i+1:nz=z+u(i,k)*x(k);endx(i)=(y(i)-z)/u(i,i);endExample:>> x=lusolv(A,b)x =0.9739 -0.0047 1.0010对称正定矩阵之Cholesky分解法：function L=Cholesky(A)%对对称正定矩阵A进行Cholesky分解n=length(A);L=zeros(n);for k=1:ndelta=A(k,k);for j=1:k-1delta=delta-L(k,j)^2;endif delta<1e-10return;endL(k,k)=sqrt(delta);for i=k+1:nL(i,k)=A(i,k);for j=1:k-1L(i,k)=L(i,k)-L(i,j)*L(k,j);endL(i,k)=L(i,k)/L(k,k);endendfunction x=Chol_Solve(A,b)%利用对称正定矩阵之Cholesky分解求解线性方程组Ax=b n=length(b);l=Cholesky(A);x=ones(1,n);y=ones(1,n);for i=1:nz=0;for k=1:i-1z=z+l(i,k)*y(k);endy(i)=(b(i)-z)/l(i,i);endfor i=n:-1:1for k=i+1:nz=z+l(k,i)*x(k);endx(i)=(y(i)-z)/l(i,i);endExample:>> a=[9 -36 30 ;-36 192 -180;30 -180 180]; >> b=[1 1 1]';>> x=Chol_Solve(a,b)x =1.8333 1.0833 0.7833对称正定矩阵之LDL’分解法：function [L,D]=LDL_Factor(A)%对称正定矩阵A进行LDL'分解n=length(A);L=eye(n);D=zeros(n);d=zeros(1,n);for k=1:nd(k)=A(k,k);for j=1:k-1d(k)=d(k)-L(k,j)*T(k,j);endif abs(d(k))<1e-10return;endfor i=k+1:nT(i,k)=A(i,k);for j=1:k-1T(i,k)=T(i,k)-T(i,j)*L(k,j);endL(i,k)=T(i,k)/d(k);endendD=diag(d);function x=LDL_Solve(A,b)%利用对对称正定矩阵A进行LDL'分解法求解线性方程组Ax=b n=length(b);[l,d]=LDL_Factor(A);y(1)=b(1);for i=2:nz=0;for k=1:i-1z=z+l(i,k)*y(k);endy(i)=b(i)-z;endx(n)=y(n)/d(n,n);for i=n-1:-1:1z=0;for k=i+1:nz=z+l(k,i)*x(k);endx(i)=y(i)/d(i,i)-z; endExample:>> x=LDL_Solve(a,b)x =1.8333 1.0833 0.78332.迭代法Richardson迭代法：function [x,n]=richason(A,b,x0,eps,M) %Richardson法求解线性方程组Ax=b %方程组系数矩阵:A%方程组之常数向量:b%迭代初始向量:X0%e解的精度控制:eps%迭代步数控制：M%返回值线性方程组的解：x%返回值迭代步数：nif(nargin == 3)eps = 1.0e-6;M = 200;elseif(nargin == 4)M = 200;endI =eye(size(A));x1=x0;x=(I-A)*x0+b;n=1;while(norm(x-x1)>eps)x1=x;x=(I-A)*x1+b;n = n + 1;if(n>=M)disp('Warning: 迭代次数太多，现在退出！');return;endendExample:>> A=[1.0170 -0.0092 0.0095;-0.0092 0.9903 0.0136;0.0095 0.0136 0.9898]; >> b=[1 0 1]';x0=[0 0 0]';>> [x,n]=richason(A,b,x0)x =0.9739-0.00471.0010n =5Jacobi迭代法：function [x,n]=jacobi(A,b,x0,eps,varargin) if nargin==3eps= 1.0e-6;M = 200;elseif nargin<3errorreturnelseif nargin ==5M = varargin{1};endD=diag(diag(A)); %求A的对角矩阵L=-tril(A,-1); %求A的下三角阵U=-triu(A,1); %求A的上三角阵B=D\(L+U);f=D\b;x=B*x0+f;n=1; %迭代次数while norm(x-x0)>=epsx0=x;x=B*x0+f;n=n+1;if(n>=M)disp('Warning: 迭代次数太多，可能不收敛！');return;endendExample:>> [x,n]=Jacobi(A,b,x0)x =0.97391.0010n =5Gauss-Seidel迭代法：function [x,n]=gauseidel(A,b,x0,eps,M) if nargin==3eps= 1.0e-6;M = 200;elseif nargin == 4M = 200;elseif nargin<3errorreturn;endD=diag(diag(A)); %求A的对角矩阵L=-tril(A,-1); %求A的下三角阵U=-triu(A,1); %求A的上三角阵f=(D-L)\b;x=G*x0+f;n=1; %迭代次数while norm(x-x0)>=epsx0=x;x=G*x0+f;n=n+1;if(n>=M)disp('Warning: 迭代次数太多，可能不收敛！');return;endendExample:>> [x,n]=gauseidel(A,b,x0)x =0.9739-0.00471.0010n =4超松驰迭代法：function [x,n]=SOR(A,b,x0,w,eps,M) if nargin==4eps= 1.0e-6;M = 200;elseif nargin<4errorreturnelseif nargin ==5M = 200;endif(w<=0 || w>=2)error;return;endD=diag(diag(A)); %求A的对角矩阵L=-tril(A,-1); %求A的下三角阵U=-triu(A,1); %求A的上三角阵B=inv(D-L*w)*((1-w)*D+w*U);f=w*inv((D-L*w))*b;x=B*x0+f;n=1; %迭代次数while norm(x-x0)>=epsx0=x;x =B*x0+f;n=n+1;if(n>=M)disp('Warning: 迭代次数太多，可能不收敛！');return;endendExample:>> [x,n]=SOR(A,b,x0,1)x =0.9739-0.00471.0010n =4对称逐次超松驰迭代法：function [x,n]=SSOR(A,b,x0,w,eps,M) if nargin==4eps= 1.0e-6;M = 200;elseif nargin<4errorreturnelseif nargin ==5M = 200;endif(w<=0 || w>=2)error;return;endD=diag(diag(A)); %求A的对角矩阵L=-tril(A,-1); %求A的下三角阵U=-triu(A,1); %求A的上三角阵B1=inv(D-L*w)*((1-w)*D+w*U);B2=inv(D-U*w)*((1-w)*D+w*L);f1=w*inv((D-L*w))*b;f2=w*inv((D-U*w))*b;x12=B1*x0+f1;x =B2*x12+f2;n=1; %迭代次数while norm(x-x0)>=epsx0=x;x12=B1*x0+f1;x =B2*x12+f2;n=n+1;if(n>=M)disp('Warning: 迭代次数太多，可能不收敛！');return;endendExample:>> [x,n]=SSOR(A,b,x0,1)x =0.9739-0.00471.0010n =3两步迭代法：function [x,n]=twostep(A,b,x0,eps,varargin)if nargin==3eps= 1.0e-6;M = 200;elseif nargin<3errorreturnelseif nargin ==5M = varargin{1};endD=diag(diag(A)); %求A的对角矩阵L=-tril(A,-1); %求A的下三角阵U=-triu(A,1); %求A的上三角阵B1=(D-L)\U;B2=(D-U)\L;f1=(D-L)\b;f2=(D-U)\b;x12=B1*x0+f1;x =B2*x12+f2;n=1; %迭代次数while norm(x-x0)>=epsx0 =x;x12=B1*x0+f1;x =B2*x12+f2;n=n+1;if(n>=M)disp('Warning: 迭代次数太多，可能不收敛！');return;endendExample:>> [x,n]=twostep(A,b,x0)x =0.9739-0.00471.0010n =3最速下降法：function [x,n]=fastdown(A,b,x0,eps) if(nargin == 3)eps = 1.0e-6;endr = b-A*x0;d = dot(r,r)/dot(A*r,r);x = x0+d*r;n=1;while(norm(x-x0)>eps)x0 = x;r = b-A*x0;d = dot(r,r)/dot(A*r,r);x = x0+d*r;n = n + 1;endExample:>> [x,n]=fastdown(A,b,x0)x =0.9739-0.00471.0010n =5共轭梯度法：function [x,n]=conjgrad(A,b,x0) if(nargin == 3)eps = 1.0e-6;endr1 = b-A*x0;p1 = r1;d = dot(r1,r1)/dot(p1,A*p1);x = x0+d*p1;r2 = r1-d*A*p1;f = dot(r2,r2)/dot(r1,r1);p2 = r2+f*p1;n = 1;for(i=1:(rank(A)-1))x0 = x;p1 = p2;r1 = r2;d = dot(r1,r1)/dot(p1,A*p1);x = x0+d*p1;r2 = r1-d*A*p1;f = dot(r2,r2)/dot(r1,r1);p2 = r2+f*p1;n = n + 1;endd = dot(r2,r2)/dot(p2,A*p2);x = x+d*p2;n = n + 1;Example:>> [x,n]=conjgrad(A,b,x0)x =0.9739-0.00471.0010n =4预处理的共轭梯度法：当AX=B为病态方程组时，共轭梯度法收敛很慢。