MATLAB计算方法3解线性方程组计算解法

.
直到(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
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消去法---计算结果可靠
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
为Crout分解。
.
为什么要讨论三角分解？若在消元法进行前能实
现三角分解 A LU ， 则
Ax b (LU)x b
Ly b (下三角方程组） Ux y (上三角方程组）
n
akj x j ]
j k 1
(k n 1, n 2,...,1)
（3.3)
.
系数矩阵为对称正定阵或严格对角占优 阵的方程组按高斯消去法计算是数值稳 定的，因而不必选主元。
严格对角占优阵：至少有一个主对角线 元素的绝对值严格大于此行或此列其他 元素的绝对值之和。
.
说明：
(1)也可采用无回代的列主元消去法(叫Gauss-
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）
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
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
L ( LK L1 )1 ，则
A LU （下三角 上三角）
（三角因子分解）
.
定义3.1 A LU 叫 A 的三角（因子）分解，其中 L是
下三角, U是上三角。
定义3.2 若 L为单位下三角阵（对角元全为1），
U 为上三角阵，则称 A LU 为Doolittle分解;
若L是L下三角，U 是单位上三角，则称 A LU
--Jordan消去法)，但比有回代的列主元消
去法的乘除运算次数多。
(2)有回代的列主元消去法所进行的乘除运算
次数为 1 n3 n2 1 n，量很小。
3
3
Gauss 列主元消去法：
优点 ------ 计算结果更可靠；
缺点 ------ 挑主元花机时更多, 次序
有变动，程x序1 ,复,杂x。n
.
用Matlab实现选列主元Gauss消去法解线性方程组
(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
行。
问 题：a11 0或 a11 0?以后各步类似。
.
用Matlab实现顺序Gauss消去法 在Matlab程序编辑器中输入：
function x=nagauss(a,b,flag) %解线形方程组ax=b，
a为系数矩阵，b为右端列向量，flag若为0，则显示中间
过程，否则不显示，默认为0，x为解向量
在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 % 消元