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

1 i n
消元：用a11消a i 1为0 :
ai 1 第1行 a 加到第i行, 第i行第j个元素成为 11 ai 1 a ij a1 j a a ij 11 （i 2,3, , n；j 1,2, , n 1）
………………………
…….
(k)
1求u的第k行：用L的第k行 u的第j列
(j k, k 1, , n) ( l k 1 , l k 2 , , l kk， 0 ,0) ( u1 j , u2 j , , u jj， 0 ,0)' a kj
m 1
l
k 1
km
由矩阵乘法 (1) 1) l11 a11 l11
（三角因子分解）
定义3.1 A LU 叫 A 的三角（因子）分解，其中 L是
下三角, U是上三角。 定义3.2 若 L为单位下三角阵（对角元全为1），
L 若 L 是下三角， U 是单位上三角，则称 A LU
U 为上三角阵，则称 A LU 为Doolittle分解;
为Crwenku.baidu.comut分解。
为什么要讨论三角分解？若在消元法进行前能实 现三角分解 A LU ， 则
Ax b ( LU ) x b
Ly b (下三角方程组） Ux y (上三角方程组）
从而容易回代求解。
1．直接三角分解法(Doolittle分解为例）
a11 a21 a n1 a12 a1n a22 a2 n an 2 a nn
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 )
直到(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)
% 回代 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消去法---计算结果可靠
(1)找行号r1 使 ar1 1 max ai 1 ，对调1 r1行：
3 3
Gauss 列主元消去法：
优点 ------ 计算结果更可靠； 缺点 ------ 挑主元花机时更多, 有变动，程序复杂。 次序
x1 ,, xn
用Matlab实现选列主元Gauss消去法解线性方程组 在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 % 消元 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);
程序运行
结果:
2．平方根法
定理3.1 设A对称正定，则有非奇异下三角阵L，使
A LLT 分解方法：设
a11 a12 a1n l11 a l a a l 21 22 2 n 21 22 an1 an 2 ann l n1 l n 2 l 其中aij a ji l11 l 21 l n1 l l 22 n2 l nn nn
例3．1 2 1 2 6 2 1 2 6 5 4 18 2 3 0 6 4 6 －3 5 5 3 －2 －1 －1 所以
1 x3 1, x2 (6 0 x3 ) / 3 2 1 x1 (6 2 x3 1 x2 ) / 2 1
( 2) 找r2，使 a r2 2 max a i 2 ，
2 i n
对调2 r2 行 .
消元：用a 22把a i 2消为0 ( i 3,4, , n) : ai 2 第2行 a 第i行，则 22 ai 2 a ij a 2 j a a ij 22 （i 3,4， , n；j 2,3， , n 1 ）
消 元： 用a11将ai 1 ( i 2,, n)化为零；
ai 1 把 a 第1行，加到第i 行。 11
( 3.1)
问 题：a11 0或 a11 0 ?以后各步类似。
用Matlab实现顺序Gauss消去法
在Matlab程序编辑器中输入：
function x=nagauss(a,b,flag) %解线形方程组ax=b， a为系数矩阵，b为右端列向量，flag若为0，则显示中间 过程，否则不显示，默认为0，x为解向量 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); if flag==0,a,end
u11 u12 u22 u1n u2 n unn
1 l 21 1 l l 1 n1 n 2
由矩阵乘法
(1) 1) 求u的第1行：用L的第一行 u的第j列 1. u1 j a1 j u1 j a1 j (j 1,2, , n) 2) 求L的第1列：用L的第i行 u的第1列 ai1 l i 1 u11 a i 1 l i 1 u11 ( i 2,3, , n)
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
（3.3)

系数矩阵为对称正定阵或严格对角占优 阵的方程组按高斯消去法计算是数值稳 定的，因而不必选主元。 严格对角占优阵：至少有一个主对角线 元素的绝对值严格大于此行或此列其他 元素的绝对值之和。

说明： (1)也可采用无回代的列主元消去法(叫Gauss--Jordan消去法)，但比有回代的列主元消 去法的乘除运算次数多。 (2)有回代的列主元消去法所进行的乘除运算 次数为 1 n 3 n 2 1 n，量很小。
( A | b) ( u | g ) (| 1) Lk Lk 1 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 （下三角 上三角）
第三章线性方程组数值解法
解线性方程组 Ax b
迭代法：格式，无穷序 列 解向量 x
直接法：理论，无舍入 误差，有限步，精确解
§3.1
直接法
一、
Gauss
消去法
a11 x1 a12 x2 a1n xn a1,n1 a21 x1 a22 x2 a2 n xn a2,n1 设有 ai1 x1 ai 2 x2 ain xn ai ,n1 an1 x1 an 2 x2 ann xn an ,n1
到此原方程组化为
a11 x1 a12 x2 a1n xn a1,n1 a22 x2 a2 n xn a2 ,n1

ai 2 x2 ain xn ai ,n1

an 2 x2 ann xn an ,n1
用Matlab实现LU分解
在Matlab程序编辑器中输入： function [L,U]=nalu(a) % a为可逆方阵；L返回单位下 三角矩阵；U返回上三角矩阵 n=length(a); U=zeros(n,n);L=eye(n,n); U(1,:)=a(1,:);L(2:n,1)=a(2:n,1)/U(1,1); for k=2:n U(k,k:n)=a(k,k:n)-L(k,1:k-1)*U(1:k-1,k:n); L(k+1:n,k)=(a(k+1:n,k)-L(k+1:n,1:k-1)*U(1:k1,k))/U(k,k); end
( 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
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消元，初等行变换，化原方程组为上三 角型。