1.2选主元三角分解
数值分析6(选主元高斯消元法)
2 1 2 2 3 3 2 1 2 0 5 0 4 6 19 9
14/27
三对角矩阵分解
f u over , with f
f u( x ) 0 x 1 f (0) a, f (1) b
20:22
( k = 2, 3, · · · ,n)
20:22
18/27
下三角方程组 LY = f
1 y1 f1 y f 1 2 2 2 n 1 y n f n
2 3 / 2 1 2
4 4 2 3 12 6 4 1 2 2 1 1
4 4 2 3 / 2 3 6 1 0 1 2 1 2
20:22
4 4 2 2 3 / 2 3 6 3 1 0 5 0 2 19 / 5 9 2
三角分解:
A = LU
三对角矩阵
=
单位下三角阵 上三角阵
AX=F LU X = F ①
20:22
L Y=F
②UX=Y
16/27
LU分解
2 2 1 1 / 2 1 3 2 2 4 2 3 5 1 2 2 1 / 2 1 / 2 5 / 2 2 4 / 5 12 / 5 2 3 5 2 1 1 / 2 1 U L 4/5 1 5 / 4 1 20:22
10 0 0
20:22
7 0 x1 7 x 6.1 0 6 2 2.5 5 2.5 x3
Guass列选主元消去法和三角分解法
Guass列选主元消去法和三⾓分解法 最近数值计算学了Guass列主消元法和三⾓分解法解线性⽅程组，具体原理如下：1、Guass列选主元消去法对于AX =B1）、消元过程：将（A|B）进⾏变换为，其中是上三⾓矩阵。
即：k从1到n-1a、列选主元选取第k列中绝对值最⼤元素作为主元。
b、换⾏c、归⼀化d、消元2）、回代过程：由解出。
2、三⾓分解法（Doolittle分解）将A分解为如下形式由矩阵乘法原理a、计算U的第⼀⾏，再计算L的第⼀列b、设已求出U的1⾄r-1⾏，L的1⾄r-1列。
先计算U的第r⾏，再计算L的第r列。
a）计算U的r⾏b）计算L的r列C#代码： 代码说明：Guass列主消元法部分将计算出来的根仍然储存在增⼴矩阵的最后⼀列，⽽Doolittle分解，将分解后的结果也储存⾄原来的数组中，这样可以节约空间。
using System;using System.Windows.Forms;namespace Test{public partial class Form1 : Form{public Form1(){InitializeComponent();}private void Cannel_Button_Click(object sender, EventArgs e){this.textBox1.Clear();this.textBox2.Clear();this.textBox3.Clear();boBox1.SelectedIndex = -1;}public double[,] GetNum(string str, int n){string[] strnum = str.Split(' ');double[,] a = new double[n, n + 1];int k = 0;for (int i = 0; i < n; i++){for (int j = 0; j < strnum.Length / n; j++){a[i, j] = double.Parse((strnum[k]).ToString());k++;}}return a;}public void Gauss(double[,] a, int n){int i, j;SelectColE(a, n);for (i = n - 1; i >= 0; i--){for (j = i + 1; j < n; j++)a[i, n] -= a[i, j] * a[j, n];a[i, n] /= a[i, i];}}//选择列主元并进⾏消元public void SelectColE(double[,] a, int n){int i, j, k, maxRowE;double temp; //⽤于记录消元时的因数for (j = 0; j < n; j++){maxRowE = j;for (i = j; i < n; i++)if (System.Math.Abs(a[i, j]) > System.Math.Abs(a[maxRowE, j]))maxRowE = i;if (maxRowE != j)swapRow(a, j, maxRowE, n); //与最⼤主元所在⾏交换//消元for (i = j + 1; i < n; i++){temp = a[i, j] / a[j, j];for (k = j; k < n + 1; k++)a[i, k] -= a[j, k] * temp;}}return;}public void swapRow(double[,] a, int m, int maxRowE, int n){int k;double temp;for (k = m; k < n + 1; k++){temp = a[m, k];a[m, k] = a[maxRowE, k];a[maxRowE, k] = temp;}}public void Doolittle(double[,] a, int n){for (int i = 0; i < n; i++){if (i == 0){for (int j = i + 1; j < n; j++)a[j, 0] = a[j, 0] / a[0, 0];}else{double temp = 0, s = 0;for (int j = i; j < n; j++){for (int k = 0; k < i; k++){temp = temp + a[i, k] * a[k, j];}a[i, j] = a[i, j] - temp;}for (int j = i + 1; j < n; j++){for (int k = 0; k < i; k++){s = s + a[j, k] * a[k, i];}a[j, i] = (a[j, i] - s) / a[i, i];}}}}private void Exit_Button_Click(object sender, EventArgs e){this.Close();}private void Confirm_Button_Click(object sender, EventArgs e){if (this.textBox2.Text.Trim().ToString().Length == 0){this.textBox2.Text = this.textBox1.Text.Trim();}else{this.textBox2.Text = this.textBox2.Text + "\r\n" + this.textBox1.Text.Trim();}this.textBox1.Clear();}private void Calculate_Button_Click(object sender, EventArgs e){string str = this.textBox2.Text.Trim().ToString();string myString = str.Replace("\n", " ").Replace("\r", string.Empty);double[,] a = new double[this.textBox2.Lines.GetUpperBound(0) + 1, this.textBox2.Lines.GetUpperBound(0) + 2];a = GetNum(myString, this.textBox2.Lines.GetUpperBound(0) + 1);if (boBox1.Text == "Guass列主消元法"){Gauss(a, this.textBox2.Lines.GetUpperBound(0) + 1);for (int i = 0; i < this.textBox2.Lines.GetUpperBound(0) + 1; i++){this.textBox3.Text = this.textBox3.Text + "\r\nX" + (i + 1) + "=" + a[i, this.textBox2.Lines.GetUpperBound(0) + 1]; }}else if (boBox1.Text == "Doolittle三⾓分解法"){this.textBox3.Enabled = true;Doolittle(a, this.textBox2.Lines.GetUpperBound(0) + 1);bel3.Text = "分解后的结果：";this.textBox3.Clear();this.textBox3.Text += "L矩阵：\r\n";for (int i = 0; i < this.textBox2.Lines.GetUpperBound(0) + 1; i++) {for (int j = 0; j < this.textBox2.Lines.GetUpperBound(0) + 1; j++) {if (j < i){this.textBox3.Text += a[i, j].ToString() + "\t";}else if (i == j){this.textBox3.Text += "1\t";}else{this.textBox3.Text += "0\t";}}this.textBox3.Text += "\r\n";}this.textBox3.Text += "\r\nU矩阵：\r\n";for (int i = 0; i < this.textBox2.Lines.GetUpperBound(0) + 1; i++) {for (int j = 0; j < this.textBox2.Lines.GetUpperBound(0) + 1; j++) {if (j >= i){this.textBox3.Text += a[i, j].ToString() + "\t";}else{this.textBox3.Text += "0\t";}}this.textBox3.Text += "\r\n";}}}private void textBox1_KeyDown(object sender, KeyEventArgs e){if (e.KeyCode == Keys.Enter){if (this.textBox1.Text.Trim().ToString().Length == 0){Calculate_Button_Click(sender, e);}else{Confirm_Button_Click(sender, e);}}}private void button1_Click(object sender, EventArgs e){this.textBox2.Enabled = true;}}} 运⾏截图： ⾄此完毕。
CH1-2选主元三角分解
选主元三角分解的思想
三角分解过程中存在的问题
Gauss消元法完成的条件是矩阵的各阶顺序过程中的除法运算要求分母不 能太小,否则
将可能产生不稳定情况. 选主元的目的就是为了完成消元且避免不稳定情况的发生
例1.2.1 ：在3位10进制的浮点数系下解方程组 0.001 1.00 x1 1.00 x 1.00 2.00 2 3.00
(k )
1 1 (r )
( k 1)
1 k
是所有元素模均小于1的 ( n k ) k 阶矩阵,
Ink 表示
n k 阶单位矩阵.
( k 1) L11 0 (k ) ( k 1) 1 L Pk L Pk Lk ( k 1) 1 % % Ln k L21 ( k 1) ( k 1) % p k 1 其中 L 是由 交换了第 1 行和 行得到的 , L 21 21
(0) (0) a % % a L 11 12 (1) 0 a L (1) (0) 22 % A L1 A M M O 0 a (1) L n2 (0) % a 1n (1) (1) a2 n A11 M 0 (1) ann (1) A12 , (1) A22
全主元消去法的具体过程
(0) (0) (1) 记A(0) = A, 设 a p q max{ aij :1 i, j n}.
0 0
若 a(0) p0 q0 0 , 则结束.
否则, 交换A(0) 的第1行与第p0行, 第1列与第q0列, (0) % 记交换后的矩阵为 A ,则 (0) (0) % A I1 p0 A(0) I1q0 PA Q1. 1 (0) % 对A 计算Gauss变换L1使得
用直接三角分解法解线性方程组
三角阵。等式左边是单位下三角阵，右边是上三角阵，要使等式
成则立L ， L只1 ,能U等于U单1，位即矩此阵三I。角于分是解L唯1一L1。
UU
1 1
I,
1 2 1
1 2 1
例7 解：
设 A 3 7
1
，试将A进行三角分解。
1 1 3
由高斯消去法得到
m21
3 1
3，m31
1 1
1
m 32
L1
1 1
0 0 1 2
例：
求
0 1 2
0 1 0
3 0 1
103的PLU分 解 。
解：用1，2， ，n的排列表示n阶置换阵P，其中排列的第i个元素
j，表示P的i行非零元素位于j列。则分解过程如下：
1 0 0 1 2
3 1 1 0 1
3 1 1 0
2 0
43
1 2
0 1 0
3 0 1
0 1 3
Ux j y j
Ly
j
bj
n1
k
n(n 1)
n2
n 次乘法
k 1
2
22
Ux j y j n k n(n 1) n2 n 次乘除法
k 1
2
22
即共需n 2 次乘除法运算。
n 2 次 乘 除 法
三角分解法的存放元素的方法：
以A (a ij )33 为例，
a11 A a21
1 mk1,k 1
k，
Lk1
1 mk1,k 1
k
mnk
1
mn,k
1
A ( L11 L21
L1 n1
)U
LU，1
a (1) 11
三角分解法
ai j =
min( i , j ) k =1
∑l
ik
uk j
ai j =
min( i , j )
一般采用列主元 对换, 将 i ,j 对换,对 j = i, i+1, …, n 有 一般采用列主元 ii a ji = ∑ l jk uki + l ji u k= k =1 法增强稳定性. 法增强稳定性.但注意 v i 1 b 也必须做相应的 l ji = ( a ji ∑ l jk uki ) / uii b 行交换. 行交换. k =1
=I
Upper-triangular
Lower-triangular With diagonal entries 1
注: L 为一般下三角阵而 U 为单位上三角阵的分解称为 单位上三角阵的分解称为 Crout 分解. 分解. ~~ 分解, 实际上只要考虑 A* 的 LU 分解,即A* = L U ,则 ~ ~ A= U * L* 即是 A 的 Crout 分解. 分解. =
(
)
n1 Step 6 Set l = 运算量为2 O(n3/6), 比普通 ann ∑ k =1 lnk ; , 比普通LU nn Step 7 Output ( lij for j = 1, …, i and i = 1, …, n );A = LDLT 分解少一半, 次开方. 分解少一半,但有 n 次开方.用
mn1
v A b
( 2) (2)
(1 (1 ( a11) a12) ... a11) n Step n 1: (2 ( v a22) ... a22) n Ln1Ln2 ... L1 A b = ... . . . (n ann)
数值线性代数第二版徐树方高立张平文上机知识题第一章实验报告
1.先用你所熟悉的的计算机语言将不选主元和列主元Gauss 消去法编写成通用的子程序；然后用你编写的程序求解84阶方程组；最后将你的计算结果与方程的精确解进行比较，并就此谈谈你对Gauss 消去法的看法。
Sol ：（1）先用matlab 将不选主元和列主元Gauss 消去法编写成通用的子程序，得到P U L ,,： 不选主元Gauss 消去法：[])(,A GaussLA U L =得到U L ,满足LU A = 列主元Gauss 消去法：[])(,,A GaussCol P U L =得到P U L ,,满足LU PA = （2）用前代法解()Pb or b Ly =，得y用回代法解y Ux =，得x求解程序为()P U L b A Gauss x ,,,,=（P 可缺省，缺省时默认为单位矩阵） （3）计算脚本为ex1_1 代码%算法1.1.3（计算三角分解：Gauss 消去法） function [L,U]=GaussLA(A) n=length(A); for k=1:n-1A(k+1:n,k)=A(k+1:n,k)/A(k,k);A(k+1:n,k+1:n)=A(k+1:n,k+1:n)-A(k+1:n,k)*A(k,k+1:n); endL=tril(A);L=L-diag(diag(L))+diag(ones(1,n));end%算法1.2.2(计算列主元三角分解：列主元Gauss消去法）function [L,U,P]=GaussCol(A)n=length(A);for k=1:n-1[s,t]=max(abs(A(k:n,k)));p=t+k-1;temp=A(k,1:n);A(k,1:n)=A(p,1:n);A(p,1:n)=temp;u(k)=p;if A(k,k)~=0A(k+1:n,k)=A(k+1:n,k)/A(k,k);A(k+1:n,k+1:n)=A(k+1:n,k+1:n)-A(k+1:n,k)*A(k,k+1:n); elsebreak;endendL=tril(A);U=triu(A);L=L-diag(diag(L))+diag(ones(1,n)); P=eye(n);for i=1:n-1temp=P(i,:);P(i,:)=P(u(i),:);P(u(i),:)=temp;endend%高斯消去法解线性方程组function x=Gauss(A,b,L,U,P)if nargin<5P=eye(length(A));endn=length(A);b=P*b;for j=1:n-1b(j)=b(j)/L(j,j);b(j+1:n)=b(j+1:n)-b(j)*L(j+1:n,j);endb(n)=b(n)/L(n,n);y=b;for j=n:-1:2y(j)=y(j)/U(j,j);y(1:j-1)=y(1:j-1)-y(j)*U(1:j-1,j);endy(1)=y(1)/U(1,1);x=y;endex1_1clc;clear;%第一题A=6*eye(84)+diag(8*ones(1,83),-1)+diag(ones(1,83),1); b=[7;15*ones(82,1);14];%不选主元Gauss消去法[L,U]=GaussLA(A);x1_1=Gauss(A,b,L,U);%列主元Gauss消去法[L,U,P]=GaussCol(A);x1_2=Gauss(A,b,L,U,P);%解的比较subplot(1,3,1);plot(1:84,x1_1,'o-');title('Gauss'); subplot(1,3,2);plot(1:84,x1_2,'.-');title('PGauss'); subplot(1,3,3);plot(1:84,ones(1,84),'*-');title('精确解');结果为（其中Gauss 表示不选主元的Gauss 消去法，PGauss 表示列主元Gauss消去法，精确解为[]'⨯8411,,1 ）：由图，显然列主元消去法与精确解更为接近，不选主元的Gauss 消去法误差比列主元消去法大，且不如列主元消去法稳定。
选主元的三角分解法
2012-2013（1）专业课程实践论文选主元（列）的三角分解法范俊，0818180124，R数学08-1班毛勇，0818180117，R数学08-1班一、算法理论从直接三角分解公式可看出，当时计算将中断或者当绝对值很小时，按分解公式计算可能要求舍入误差的累积，但如果非奇异，就可通过交换的行实现矩阵的分解，因此可采用与列主元素消去法类似的方法（可以证明下述方法与列主元素消去法等价），将直接三角分解法修改为（部分）选主元的三角分解法。
设第步分解已完成，这时有第步分解需要用到式及式为了避免用小的数作除数，引进量于是有,,用作为，交换的行与行元素(将位置的新元素仍记作及)，于是有。
由此在进行第步分解计算。
该程序即利用选主元的三角分解法解方程求方程的根。
先选择列主元，再构造LU矩阵，通过求解LY = B和UX = Y得出需要的解。
注：方程维数在程序中需按题意自定义。
否源代码源代码：LU_Decomposition.cpp#include<iostream.h>#include<math.h>#define N 4 //矩阵维数，可自定义static double A[N][N]; //系数矩阵static double B[N]; //右端项static double Y[N]; //中间项static double X[N]; //输出static double S[N]; //选取列主元的比较器int i,j,k; //计数器void main(){cout << "请输入线性方程组（ai1，ai2，ai3......ain, yi）："<<endl;for ( i = 0; i < N ;i++){for (int j = 0; j< N ;j++ )cin >> A[i][j];cin >>B[i];}for (k = 0 ; k < N; k++){//选列主元int index =k;for (i = k ; i < N; i++){double temp = 0;for (int m = 0 ; m < k ;m++){temp = temp + A[i][m]*A[m][k];}S[i] = A[i][k] - temp;if(S[index] < S[i]){index = i;}} //交换行double temp;for( i = k ; i < N ;i++ ){temp = A[index][i];A[index][i] = A[k][i];A[k][i] = temp;}temp = B[index];B[index] = B[k];B[k] = temp; // 构造L、U矩阵for (j = k ; j < N; j++){double temp = 0;for (int m = 0 ;m < k; m++ ){temp = temp + A[k][m] * A[m][j];}A[k][j] = A[k][j] - temp; //先构造U一行的向量}for( i = k+1; i < N; i++){double temp = 0;for (int m =0 ; m < k; m++ ){temp = temp + A[i][m] * A[m][k];}A[i][k] = (A[i][k] - temp)/A[k][k]; //再构造L一列的向量}}//求解L Y = BY[0] = B[0];for (i = 1; i < N ; i++){double temp = 0;for (int j =0 ; j < i; j++ ){temp = temp + A[i][j] * Y[j];}Y[i] = B[i] - temp;}//求解UX = YX[N-1] = Y[N-1]/A[N-1][N-1];for (i = N-2 ; i >= 0 ; i-- ){double temp = 0;for (int j =i+1 ; j < N; j++ ){temp = temp + A[i][j] * X[j];}X[i] = (Y[i] - temp)/A[i][i];}//打印Xcout << "线性方程组的解（X1，X2，X3......Xn）为："<<endl; for( i = 0; i < N ;i++){cout << X[i] <<" ";}}四、算法实现例1. 用选主元（列）的三角分解法解解：利用软件的解答方法为：（1）先将程序中的维数N改为3。
列主元高斯消去法和列主元三角分解法解线性方程
计算方法实验报告1课题名称用列主元高斯消去法和列主元三角分解法解线性方程目的和意义高斯消去法是一个古老的求解线性方程组的方法,但由它改进得到的选主元的高斯消去法则是目前计算机上常用的解低阶稠密矩阵方程组的有效方法;用高斯消去法解线性方程组的基本思想时用矩阵行的初等变换将系数矩阵A 约化为具有简单形式的矩阵上三角矩阵、单位矩阵等,而三角形方程组则可以直接回带求解 用高斯消去法解线性方程组b Ax =其中A ∈Rn ×n 的计算量为：乘除法运算步骤为32(1)(1)(21)(1)(1)262233n n n n n n n n n n nMD n ----+=+++=+-,加减运算步骤为(1)(21)(1)(1)(1)(25)6226n n n n n n n n n n AS -----+=++=;相比之下,传统的克莱姆法则则较为繁琐,如求解20阶线性方程组,克莱姆法则大约要19510⨯次乘法,而用高斯消去法只需要3060次乘除法;在高斯消去法运算的过程中,如果出现absAi,i 等于零或过小的情况,则会导致矩阵元素数量级严重增长和舍入误差的扩散,使得最后的计算结果不可靠,所以目前计算机上常用的解低阶稠密矩阵方程的快速有效的方法时列主元高斯消去法,从而使计算结果更加精确; 2、列主元三角分解法高斯消去法的消去过程,实质上是将A 分解为两个三角矩阵的乘积A=LU,并求解Ly=b 的过程;回带过程就是求解上三角方程组Ux=y;所以在实际的运算中,矩阵L 和U 可以直接计算出,而不需要任何中间步骤,从而在计算过程中将高斯消去法的步骤进行了进一步的简略,大大提高了运算速度,这就是三角分解法采用选主元的方式与列主元高斯消去法一样,也是为了避免除数过小,从而保证了计算的精确度计算公式1、 列主元高斯消去法设有线性方程组Ax=b,其中设A 为非奇异矩阵;方程组的增广矩阵为第1步k=1：首先在A 的第一列中选取绝对值最大的元素1l a ,作为第一步的主元素:111211212222112[,]n n n l n nn n a a a a b a a a b a a a b ⎡⎤⎢⎥⎢⎥⎢⎥=⎢⎥⎢⎥⎢⎥⎢⎥⎢⎥⎣⎦a b然后交换A,b 的第1行与第l 行元素,再进行消元计算;设列主元素消去法已经完成第1步到第k -1步的按列选主元,交换两行,消元计算得到与原方程组等价的方程组 Akx=bk第k 步计算如下：对于k=1,2,…,n -11按列选主元：即确定t 使 2如果t ≠k,则交换A,b 第t 行与第k 行元素; 3消元计算消元乘数mik 满足:4回代求解2、 列主元三角分解法 对方程组的增广矩阵 经过k -1步分解后,可变成如下形式：111max 0l i i n a a ≤≤=≠(1)(1)(1)(1)(1)1112111(2)(2)(2)(2)22222()(()1)()()()()()1,1()(,)()[,][,] k k k k nk k nk n k k k k k kk kn k k k k n k k k n nn a a a a b a a a b a a b a b b a a a +++⎡⎤⎢⎥⎢⎥⎢⎥⎢⎥→=⎢⎥⎢⎥⎢⎥⎢⎥⎢⎥⎢⎥⎣⎦A b A b ()()max 0k k tk ik k i na a ≤≤=≠,(1,,)ik ik ik kka a m i k n a ←=-=+, (,1,,), (1,,)ij ij ik kji i ik k a a m a i j k n b b m b i k n ←+=+⎧⎨←+=+⎩⎪⎪⎩⎪⎪⎨⎧--=-←←∑+=)1,,2,1(,)(1n n i a x a b x a b x ii n i j j ij i i nnn n [,]A A b =11121,11111222,122221,11,1,1,211,11,2121,112,112,1k k k k k k k j n k k j n k k k i i i k n n kk kj kn k ik ij in i nknjk k k j k n n nnk k n a a a b A a u u u u u u y l l l l l l ll l l l u u u u u y u u u u y a a b a a b l a -------------⎡→⎣⎤⎢⎥⎢⎥⎢⎥⎢⎥⎢⎥⎢⎥⎢⎥⎢⎥⎢⎥⎢⎥⎢⎥⎢⎥⎦第k 步分解,为了避免用绝对值很小的数kku 作除数,引进量1111 (,1,,;1,2,,) ()/ (1,2,,;1,2,,)k kj kj km mj m k ik ik im mk kkm u a l u j k k n k n l a l u u i k k n k n -=-=⎧=-=+=⎪⎪⎨⎪=-=++=⎪⎩∑∑11(,1,,)k i ik im mk m s a l u i k k n -==-=+∑,于是有kk u =ks ;如果 ,则将矩阵的第t 行与第k 行元素互换,将i,j 位置的新元素仍记为jjl 或jja ,然后再做第k 步分解,这时列主元高斯消去法程序流程图max t ik i n s s ≤≤= ()/ 1,2,,)1 (1,2,,),kk k k t iki k ik u s s s l s s i k k n l i k k n ===++≤=++即交换前的，（且列主元高斯消去法Matlab主程序function x=gauss1A,b,c %列主元法高斯消去法解线性方程Ax=bif lengthA~=lengthb %判断输入的方程组是否有误disp'输入方程有误'return;enddisp'原方程为AX=b：' %显示方程组Abdisp'------------------------'n=lengthA;for k=1:n-1 %找列主元p,q=maxabsAk:n,k; %找出第k列中的最大值,其下标为p,qq=q+k-1; %q在Ak:n,k中的行号转换为在A中的行号if absp<cdisp'列元素太小,detA≈0';break;elseif q>ktemp1=Ak,:; %列主元所在行不是当前行,将当前行与列主Ak,:=Aq,:; 元所在行交换包括bAq,:=temp1;temp2=bk,:;bk,:=bq,:;bq,:=temp2;end%消元for i=k+1:nmi,k=Ai,k/Ak,k; %Ak,k将Ai,k消为0所乘系数Ai,k:n=Ai,k:n-mi,kAk,k:n; %第i行消元处理bi=bi-mi,kbk; %b消元处理endenddisp'消元后所得到的上三角阵是'A %显示消元后的系数矩阵bn=bn/An,n; %回代求解for i=n-1:-1:1bi=bi-sumAi,i+1:nbi+1:n/Ai,i;endclear x;disp'AX=b的解x是' x=b;调用函数解题列主元三角分解法程序流程图列主元三角分解法Matlab主程序①自己编的程序：function x=PLUA,b,eps %定义函数列主元三角分解法函数if lengthA~=lengthb %判断输入的方程组是否有误disp'输入方程有误'return;enddisp'原方程为AX=b：' %显示方程组Abdisp'------------------------'n=lengthA;A=A b; %将A与b合并,得到增广矩阵for r=1:nif r==1for i=1:nc d=maxabsA:,1; %选取最大列向量,并做行交换if c<=eps %最大值小于e,主元太小,程序结束break;elseendd=d+1-1;p=A1,:;A1,:=Ad,:;Ad,:=p;A1,i=A1,i;endA1,2:n=A1,2:n;A2:n,1=A2:n,1/A1,1; %求u1,ielseur,r=Ar,r-Ar,1:r-1A1:r-1,r; %按照方程求取ur,iif absur,r<=eps %如果ur,r小于e,则交换行p=Ar,:;Ar,:=Ar+1,:;Ar+1,:=p;elseendfor i=r:nAr,i=Ar,i-Ar,1:r-1A1:r-1,i; %根据公式求解,并把结果存在矩阵A中endfor i=r+1:nAi,r=Ai,r-Ai,1:r-1A1:r-1,r/Ar,r; %根据公式求解,并把结果存在矩阵A中endendendy1=A1,n+1;for i=2:nh=0;for k=1:i-1h=h+Ai,kyk;endyi=Ai,n+1-h; %根据公式求解yiendxn=yn/An,n;for i=n-1:-1:1h=0;for k=i+1:nh=h+Ai,kxk;endxi=yi-h/Ai,i; %根据公式求解xiendAdisp'AX=b的解x是'x=x'; %输出方程的解②可直接得到P,L,U并解出方程解的的程序查阅资料得子函数PLU1,其作用是将矩阵A分解成L乘以U的形式;PLU2为调用PLU1解题的程序,是自己编的Ⅰ.function l,u,p=PLU1A %定义子函数,其功能为列主元三角分解系数矩阵A m,n=sizeA; %判断系数矩阵是否为方阵if m~=nerror'矩阵不是方阵'returnendif detA==0 %判断系数矩阵能否被三角分解error'矩阵不能被三角分解'endu=A;p=eyem;l=eyem; %将系数矩阵三角分解,分别求出P,L,Ufor i=1:mfor j=i:mtj=uj,i;for k=1:i-1tj=tj-uj,kuk,i;endenda=i;b=absti;for j=i+1:mif b<abstjb=abstj;a=j;endendif a~=ifor j=1:mc=ui,j;ui,j=ua,j;ua,j=c;endfor j=1:mc=pi,j;pi,j=pa,j;pa,j=c;endc=ta;ta=ti;ti=c;endui,i=ti;for j=i+1:muj,i=tj/ti;endfor j=i+1:mfor k=1:i-1ui,j=ui,j-ui,kuk,j;endendendl=trilu,-1+eyem;u=triuu,0Ⅱ.function x=PLU2A,b %定义列主元三角分解法的函数l,u,p=PLU1A %调用PLU分解系数矩阵A m=lengthA; %由于A左乘p,故b也要左乘p v=b;for q=1:mbq=sumpq,1:mv1:m,1;endb1=b1 %求解方程Ly=b for i=2:1:mbi=bi-sumli,1:i-1b1:i-1;endbm=bm/um,m; %求解方程Ux=y for i=m-1:-1:1bi=bi-sumui,i+1:mbi+1:m/ui,i;endclear x;disp'AX=b的解x是' x=b;调用函数解题①②编程疑难这是第一次用matlab编程,对matlab的语句还不是非常熟悉,因此在编程过程中,出现了许多错误提示;并且此次编程的两种方法对矩阵的运算也比较复杂;问题主要集中在循环控制中,循环次数多了一次或者缺少了一次,导致数据错误,一些基本的编程语句在语法上也会由于生疏而产生许多问题,但是语句的错误由于系统会提示,比较容易进行修改,数据计算过程中的一些逻辑错误,比如循环变量的控制,这些系统不会提示错误,需要我们细心去发现错误,不断修正,调试;。
列主元三角分解法步骤
列主元三角分解法步骤嘿，朋友们！今天咱来唠唠列主元三角分解法那些事儿。
你说这列主元三角分解法啊，就好比是搭积木，咱得一层一层稳稳地搭起来。
咱先说说啥是列主元。
这就好像是在一群数字小伙伴里，挑出那个最厉害、最能带头的数字来！它就是那个主元啦。
然后呢，咱就开始分解啦。
这就像是把一个大难题，一点点地拆解成小步骤，让它变得不再那么可怕。
你看啊，我们要通过巧妙的变换，让矩阵变得有秩序，就像整理房间一样，把东西都归置得整整齐齐。
比如说，咱有个矩阵，里面的数字就像一群调皮的小孩子，到处乱跑。
那我们就得想办法让它们听话呀！找到那个主元，把它放在合适的位置，然后其他数字就会跟着它的节奏来。
这过程中，可不能马虎哦！就像走路一样，一步一步都得走稳了。
要是有一步走错了，那可能就前功尽弃啦。
你想想，如果我们不仔细地找主元，随便乱来，那结果能对吗？那肯定不行呀！而且哦，这列主元三角分解法用处可大啦！在解决很多数学问题的时候，它就像是一把神奇的钥匙，能打开那扇困住我们的门。
它能让复杂的计算变得简单起来，就像给我们加了一双翅膀，能在数学的天空中自由翱翔。
咱再打个比方，它就像是一个好向导，带着我们在数学的迷宫里找到正确的路。
这可不是随便说说的呀！很多厉害的数学家都在用它呢！他们用它解决了一个又一个难题，创造了一个又一个奇迹。
咱虽然不是数学家，但咱也可以学会呀！只要我们有耐心，肯钻研，就一定能掌握这神奇的列主元三角分解法。
所以啊，朋友们，别害怕数学，别害怕这列主元三角分解法。
它其实没那么难，只要我们用心去学，去理解，就一定能学好它。
让我们一起加油，在数学的海洋里畅游吧！就这么定啦！。
列主元高斯消去法和列主元三角分解法解线性方程
b(k,:)=b(q，:)；
b(q,:）=ｔemｐ2;
ｅｎｄ
%消元
forｉ=k＋1：n
m(i,ｋ)=A(i,k)/A（ｋ,k);%Ａ（k,k)将Ａ(i,k)消为0所乘系数
A(i，k:n)=A(ｉ,k:ｎ)-m（i,k)*A（k,k:n);％第i行消元处理
b(i)=b(i)－m(ｉ，k）＊ｂ(k);%ｂ消元处理
用高斯消去法解线性方程组的基本思想时用矩阵行的初等变换将系数矩阵A约化为具有简单形式的矩阵(上三角矩阵、单位矩阵等)，而三角形方程组则可以直接回带求解
用高斯消去法解线性方程组 （其中A∈Rn×n)的计算量为：乘除法运算步骤为 ,加减运算步骤为 。相比之下，传统的克莱姆法则则较为繁琐，如求解20阶线性方程组,克莱姆法则大约要 次乘法，而用高斯消去法只需要３０６０次乘除法。
对方程组的增广矩阵经过k-1步分解后，可变成如下形式:
第k步分解，为了避免用绝对值很小的数 作除数,引进量
，于是有 = 。如果,则将矩阵的第t行与第k行元素互换,将（i,j）位置的新元素仍记为 或 ，然后再做第ｋ步分解,这时
【列主元高斯消去法程序流程图】
【列主元高斯消去法Mａtｌab主程序】
ｆｕnction x=gauss1(Ａ，b,c)%列主元法高斯消去法解线性方程Ax=ｂ
【列主元三角分解法程序流程图】
【列主元三角分解法Ｍatlab主程序】
①自己编的程序：
function x=PＬU（A,ｂ,eｐｓ)%定义函数列主元三角分解法函数
if（length(A)～=length(b))％判断输入的方程组是否有误
ｄisp('输入方程有误!＇)
retｕｒn;
end
选主元三角分解
§1.2 选主元三角分解1.2.1 问题的提出大家知道，对于方程组来说，只要A非奇异方程组就存在唯一的解。
然而，A非奇异并不能保证其顺序主子阵均非奇异。
因此，A非奇异并不能保证Gauss消去过程能够进行到底。
这样，我们的问题自然便是，怎样修改算法1·1·3才能使其适应于非奇异矩阵呢？此外，在算法1·1·3中计算时，位于分母上的主元虽不为零但很小时，是否会对算法产生不良影响呢？若有影响，该如何解决？下面来看一个例子。
例1·2·1假定我们是在3位10进制数的浮点数系下来解方程组用算法1·1·3得，，从而得该方程组的计算解为，这与精确解相差甚远。
上例中的问题是由小主元引起的。
当然，如果用更高精度的计算机来计算，可使计算解的精度提高，然而仅以提高计算机精度的方法去解决这个问题是不明智的，因为计算机的精度毕竟是有限的。
事实上，我们可以用下面的方法来避免小主元的出现。
如果交换例1·2·1的第一个与第二个方程的位置，原方程组变为再用算法1·1·3在同样的数系下进行计算，可得进而得原方程组的计算解为，这已与方程组的精确解相当接近了。
1.2.2 全主元三角分解这种交换方程顺序的方法其实并不是解决小主元总是的唯一方法，当然亦可通过交换未知向量x在方程组中的顺序来解决。
这样，如果出现小的主元，我们就可以选择一个合适的元素，交换A的行和列，将此元素换到主元位置上。
例如在第k步中，若太小，并且选择作为主元，则我们需要先交换第k行，再交换第k列和第q列，从而将的移至（k,k）位置上，消去过程用新的主元继续进行。
为了不打乱在消去过程中已经引入的零元素的分布，所选的的位置应该满足。
为了下面叙述简单起见，我们初等置换矩阵，它是单位矩阵I的第p列与第q列交换所得到的矩阵，即用左乘矩阵A便交换了A的第p行与第q行；用右乘A便交换了A的第p列与第q列。
列主元三角分解法分解三阶矩阵
列主元三角分解法分解三阶矩阵1.引言1.1 概述列主元三角分解法是一种经典的数值计算方法，用于将一个三阶矩阵分解为一个下三角矩阵和一个上三角矩阵的乘积形式。
这种分解方法可以帮助我们解决线性方程组和求逆矩阵等数值计算问题。
在实际问题中，我们经常会遇到需要求解线性方程组的情况。
而列主元三角分解法的主要作用就是将线性方程组的求解转化为两个步骤：矩阵分解和回代求解。
通过将一个矩阵分解为一个下三角矩阵和一个上三角矩阵的乘积形式，我们可以简化线性方程组的求解过程，提高计算效率。
列主元三角分解法的步骤包括：选取列主元、消元和回代。
其中，选取列主元的过程是为了减小计算误差，保证数值计算的稳定性。
消元过程则是通过逐行操作，将原始矩阵逐步转化为下三角和上三角矩阵的乘积形式。
回代过程是求解三角方程组，得到线性方程组的解。
在本篇文章中，我们将详细介绍列主元三角分解法的原理和步骤。
我们将首先讲解列主元三角分解法的原理，包括选取列主元的方法和消元过程的具体操作。
然后，我们将详细介绍回代过程，以及列主元三角分解法的优点和应用。
通过本文的学习，读者将能够了解到列主元三角分解法的基本原理和操作步骤，掌握如何应用列主元三角分解法求解线性方程组和求逆矩阵。
同时，读者还能够了解到列主元三角分解法在实际问题中的重要性和广泛应用，为进一步深入学习数值计算提供基础知识和思路。
1.2文章结构1.2 文章结构本文将按照以下结构来进行阐述列主元三角分解法分解三阶矩阵的原理、步骤以及应用。
第一部分，引言，将对列主元三角分解法进行概述。
首先介绍三阶矩阵的基本概念和性质，然后引出列主元三角分解法的出发点和主要思想。
通过对该方法的简要介绍，读者将能够掌握本文所要介绍的内容。
第二部分，正文，将详细介绍列主元三角分解法的原理和步骤。
首先，我们将解释列主元三角分解法的原理，包括如何选择主元元素和使用主元消去的思想。
接着，我们将逐步阐述列主元三角分解法的具体步骤，包括将矩阵转化为上三角矩阵和求解最终的解向量。
三角分解
列主元素高斯消去法相当于先进行一系列行交换后再对 PAX Pb 应用顺序高斯消去法.
定理8(列主元三角分解) 若A为非奇异矩阵, 则存在排列 矩阵P使得 PA LU 其中L为单位下三角矩阵,U为上三角矩阵. 说明: L, U, Ip的存贮.
§7.4
高斯消去法的变形
设有线性方程组：AX=b
a11 a12 a1n x1 b1 a x b a22 a2 n 21 , X 2 , b 2 . A 如何简单的实现三 a x 角分解？ bn an 2 ann n1 n
本节主要内容
1、回顾高斯消去法与三角分解的关系 2、三角分解的条件、方法与应用
下面用矩阵描述列主元消去法
L1I1,i1 A(1) A( 2) , L1I1,i1b(1) b( 2) ,, Lk I k ,ik A( k ) A( k 1) , Lk I k ,ik b( k ) b( k 1) .
其中I k ,ik 为初等置换阵.
于是 Ln 1I n 1,in1 L2 I 2,i2 L1I1,i1 A A( n ) U . ~ ~ P为排列矩阵 即 P A U , P b b( n ) . ~ L为单位下三角矩阵 下面就n 4考察P .
UA
( 4)
L3 I3,i3 L2 I 2,i2 L1I1,i1 A
PA 选主元三角分解算法： A, 整型Ip(n)记录主行, x b.
1. 对r 1,2,, n, r 1 (1) 计算si : air si air lik ukr (i r ,, n). k 1 (2) 选主元：取ir 使得 sir max si ，Ip(r ) ir
全主元的三角分解定理
全主元的三角分解定理
全主元的三角分解定理是一种矩阵分解方法，用于将一个矩阵分解成一个上三角矩阵和一个下三角矩阵的乘积。
全主元的三角分解定理可以表示为：如果一个n\times n的矩阵A是非奇异的，那么存在一个主元矩阵P和一个上三角矩阵U，以及一个下三角矩阵L，使得A=PLU。
其中，主元矩阵P是一个置换矩阵，它的每行和每列都只有一个$1$，其他元素都是$0$。
上三角矩阵U的主对角线元素都是1，其他元素都是0。
下三角矩阵$L$的主对角线元素都是1，其他元素都是0。
全主元的三角分解定理的证明可以通过高斯消元法来完成。
具体来说，我们可以通过一系列的消元操作，将矩阵$A$转化为一个上三角矩阵U，然后通过回代操作，将U 转化为一个下三角矩阵L。
在这个过程中，我们需要使用主元来保证消元操作的正确性。
全主元的三角分解定理在数值计算中有很多应用，例如求解线性方程组、计算矩阵的逆、计算特征值和特征向量等。
它的优点是可以避免在计算过程中出现除以零的情况，从而提高计算的稳定性和精度。
需要注意的是，全主元的三角分解定理的计算量比较大，因此在实际应用中，我们通常使用部分主元的三角分解定理来减少计算量。
线性代数方程组数值解法及MATLAB实现综述
线性代数方程组数值解法及MATLAB 实现综述廖淑芳 20122090 数计学院 12计算机科学与技术1班（职教本科） 一、分析课题随着科学技术的发展，提出了大量复杂的数值计算问题，在建立电子计算机成为数值计算的主要工具以后，它以数字计算机求解数学问题的理论和方法为研究对象。
其数值计算中线性代数方程的求解问题就广泛应用于各种工程技术方面。
因此在各种数据处理中，线性代数方程组的求解是最常见的问题之一。
关于线性代数方程组的数值解法一般分为两大类：直接法和迭代法。
直接法就是经过有限步算术运算，可求的线性方程组精确解的方法（若计算过程没有舍入误差），但实际犹如舍入误差的存在和影响，这种方法也只能求得近似解，这类方法是解低阶稠密矩阵方程组级某些大型稀疏矩阵方程组的有效方法。
直接法包括高斯消元法，矩阵三角分解法、追赶法、平方根法。
迭代法就是利用某种极限过程去逐步逼近线性方程组精确解的方法。
迭代法具有需要计算机的存储单元少，程序设计简单，原始系数矩阵在计算过程始终不变等优点，但存在收敛性级收敛速度问题。
迭代法是解大型稀疏矩阵方程组（尤其是微分方程离散后得到的大型方程组）的重要方法。
迭代法包括Jacobi 法SOR 法、SSOR 法等多种方法。
二、研究课题-线性代数方程组数值解法 一、 直接法 1、 Gauss 消元法通过一系列的加减消元运算，也就是代数中的加减消去法，以使A 对角线以下的元素化为零，将方程组化为上三角矩阵；然后，再逐一回代求解出x 向量。
1.1消元过程1. 高斯消元法(加减消元)：首先将A 化为上三角阵，再回代求解。
11121121222212n n n n nn n a a a b a a a b a a a b ⎛⎫ ⎪ ⎪ ⎪ ⎪⎝⎭L L M M O M M L (1)(1)(1)(1)(1)11121311(2)(2)(2)(2)222322(3)(3)(3)3333()()000000nn n n n nn n a a a a b a a a b a a b a b ⎛⎫ ⎪ ⎪ ⎪ ⎪⎪ ⎪⎝⎭L L L M M M OM M L 步骤如下：第一步：1111,2,,i a i i n a -⨯+=L 第行第行 11121121222212n nn n nn n a a a b a a a b a a a b ⎛⎫⎪ ⎪ ⎪⎪⎝⎭L L M M O M M L111211(2)(2)(2)2222(2)(2)(2)200n nn nn n a a a b a a b a a b ⎛⎫⎪⎪ ⎪ ⎪⎝⎭LL M M O M M L第二步：(2)2(2)222,3,,i a i i n a -⨯+=L 第行第行111211(2)(2)(2)2222(2)(2)(2)200nnn nn n a a a b a a b a a b ⎛⎫ ⎪⎪ ⎪ ⎪⎝⎭L L M M O M M L11121311(2)(2)(2)(2)222322(3)(3)(3)3333(3)(3)(3)300000n n nn nn n a a a a b a a a b a a b a a b ⎛⎫⎪ ⎪ ⎪⎪⎪ ⎪⎝⎭LL LM M M O M M L 类似的做下去，我们有：第k 步：()()k ,1,,k ikk kka i i k n a -⨯+=+L 第行第行。
方程组直接三角分解法
如果U的第1至k-1列和L的第1至k-1列已经算出，则由
ak1 lk1 , k 2,3, u11
k r 1
, n.
（4.2.2）
a l u ,j k , k 1 , , n , kj kr rj
可得U的第k行元素
ukj =akj 同理，由
k
k 1 r 1
l kr u rj
用向后回代的方法即可求得x。设x=（x1 ，x2， · · · xn) T, y=（y1, y2, · · · yn) T,b= （b1 ，b2， · · · bn) T, 则有计算公式
y b 1 1 i 1 （4.2.5） y b l y 1 , 2 ,..., n i i ir r ,i r 1
(4.2.8)
利用(4.2.7)和(4.2.7)可得
u1 b 1 n li ai / ui1, i 2,3,... u b l c , i 2,3,... n i i i1 i
(4.2.9)
由此可求得L和U的所有元素.。解原方程组Ax=b可分为两步Ly=d 和Ux=y,计算公式为
由于方车程组的右端参与了消元计算，所以Ly=Pb的解为y=b（3）= （20，14/3，216/39） T 。解Ux=y得x=（1，2，3） T
4.2.2
三对角方程组的追赶法
b a1 A c1 b c2 a n1 bn1 an c n1 bn
(k ) akk
uk 1,k 1 l k , k 1 l n , k 1
(k ) ank
u1n u2 n uk 1,n (k ) akn (k ) ann
数值分析第4讲矩阵的三角分解法
ln2
u1n u2n unn
数值分析
例：用直角三角分解法解
1 2 3 x1 14
252x2源自183 1 5 x3 20
解：用分解计算公式得
1 0 0 1 2 3
A 2
1
0 0 1
4
LU .
3 5 1 0 0 24
求解 Ly (14,18, 20)T y (14, 10, 72)T ,
n
n3/3
GaussJordan elimination
n3/2
Square root
/improved square root
n3/6
追赶法
5n-4
Row/column/ complete Pivoting
Only eliminating elements in a
column below the diagonal one
k 1 n
4. bn bn / unn , bi (bi uikbk ) / uii k i 1
(i n 1, ,1).
A1 U 1L1P
数值分析
数值分析
数值分析
三对角阵的追赶法（A的前n-1个顺序主子式非零）
在数值求解常微分方程边值问题、热传导方程和建 立三次样条函数时，都会要解三对角方程组：AX=b
r 1 及 urr arr lrkukr .
k 1
希望urr为大的数！
数值分析
为了避免用小的数urr作除数，引进量
r 1 si air likukr (i r, , n).
k 1
当不选主元时，
urr sr , lir si / sr (i r 1, , n).
当选主元时，取ir使得
列主元高斯消去法和列主元三角分解法解线性方程组
列主元高斯消去法和列主元三角分解法解线性方程组一、课题名称列主元高斯消去法和列主元三角分解法解线性方程组二、班级:周二7—10，周四3、4三、目的和意义1、进一步熟悉matlab的编程环境。
2、体会列主元高斯消去法和列主元三角分解法解线性方程组在matlab环境下的编程实现。
3、列主元高斯消去法是解线性方程组的一般方法，列主元三角分解法的计算量比列主元高斯消去法少，更适用于解大规模矩阵方程。
四、流程图1、列主元高斯消去法:输入A，Bn=A的阶A，B是否匹配停止K=1,2,…，n-1按列选主元，设k列第max行最大将A，B的第k行与第max行交换输出”detA=0”A(k,k)=0 停止消元l=k+1,k+2,…,nA(l,k)=A(l,k)/A(k,k)l,m=k+1,k+2,…,nA(l,m)=A(l,m)-A(k,m)*A(l,k) B(l)=B(l)-B(k)*A(l,k)回代求解X(n)=B(n)/A(n,n)k=n-1,n-2,…,1sum=0l=k+1,k+2,…,nsum=sum+A(k,l)*X(l)X(k)=(B(k)-sum)/A(k,k)输出X停止2、列主元三角分解法输入A，Bn=A的阶A，B是否匹配停止是第一行是否要交换交换否计算U的第一行元素计算L的第一列元素r=2,3,…,n计算s(i)(i=r,r+1,…,n)是第r行是否要交换交换否计算U的第r行元素计算L的第r列元素回代求解Ly=B回代求解Ux=y输出X停止五、matlab执行结果六、结果讨论和分析Matlab编写列主元高斯消去法和列主元三角分解法解线性方程组的程序有很多循环，必须先画流程图，清晰地理清每一层循环的含义，才能得到正确的结果。
三角分解法
三角分解法
三角分解法亦称因子分解法，由消元法演变而来的解线性方程组的一类方法。
设方程组的矩阵形式为Ax=b，三角分解法就是将系数矩阵A分解为一个下三角矩阵L和一个上三角矩阵U之积：A=LU，然后依次解两个三角形方程组Ly=b 和Ux=y，而得到原方程组的解，例如，杜利特尔分解法、乔莱斯基分解法等就是三角分解法。
【基本介绍】
若能通过正交变换，将系数矩阵A分解为A=LU，其中L是单位下三角矩阵(主对角线元素为1的下三角矩阵)，而U是上三角矩阵，则线性方程组Ax=b 变为LUx=b，若令Ux=y，则线性方程组Ax=b的求解分为两个三角方程组的求解：
(1)求解Ly=b，得y；
(2)再求解Ux=y，即得方程组的解x；
因此三角分解法的关键问题在于系数矩阵A的LU分解
【矩阵LU能分解的充分条件】
一般地，任给一个矩阵不一定有LU分解，下面给出一个矩阵能LU分解的充分条件。
定理1 对任意矩阵AϵR n×n（n≥2），若A的各阶顺序主子式均不为零，则A有唯一的Doolittle分解(或Crout分解)。
定理2 若矩阵AϵR n×n（n≥2）非奇异，且其LU分解存在，则A的LU 分解是唯一的。
- 1、下载文档前请自行甄别文档内容的完整性，平台不提供额外的编辑、内容补充、找答案等附加服务。
- 2、"仅部分预览"的文档,不可在线预览部分如存在完整性等问题,可反馈申请退款(可完整预览的文档不适用该条件!)。
- 3、如文档侵犯您的权益，请联系客服反馈,我们会尽快为您处理(人工客服工作时间：9:00-18:30)。
(k A22 −1)
设 a (pk −1) q 若a
k −1 k −1
( = max{ aijk −1) : k ≤ i, j ≤ n}.
选
( k −1) pk −1qk −1
= 0 , 则消去过程结束.
}
元
主
否则, 交换A(k-1) 的第k行与第pk-1行, 第k列与第qk-1列,
% 记交换后的矩阵为 A( k −1) , 则
选 主 元
% A(1)
Gauss
A
(2)
换L2
% = L2 A(1)
(2) A11 = 0 (2) A12 , (2) A22
(2) A11 为2
阵, 则A(2) = L2P2L1P1AQ1Q2.
Gauss 消 去 法
(3) 假定已求出k－1个Gauss变换L1 , L2 , ··· , Lk－1∈Rn×n 和 2(k－1)个初等置换矩阵P1 , P2 , ··· , Pk－1 和 Q1 , Q2 , ··· , Qk－1 ∈Rn×n 使得
1 × 1.00 − 0.001 × 2.00 = 0.100 × 10 − 0.0002 × 10 = 1.00 0 y1 3.00 y1 = 3.00 1 y = ⇒ y2 = 1.00 − 0.003 0.001 1 2 1.00 = 0.100 × 10 − 0.0003 × 10 = 1.00 x1 = 1.00 1.00 2.00 x1 3.00 x = ⇒ x = 1.00 0 1.00 的第p列与第q列交换所得到的矩阵. Ipq = (e1 ··· ep-1 eq ep+1 ··· eq-1 ep eq+1 en), (p≤q) IpqA : 交换A的第p行与第q行 AIpq : 交换A的第p列与第q列
全主元消去法的具体过程
(0) (0) (1) 记A(0) = A, 设 a p q = max{ aij :1 ≤ i, j ≤ n}.
A( k −1) = Lk −1 Pk −1 L L1 P AQ1 L Qk −1 1
( 其中A11k −1) 为k－1阶上三角阵, 且
(k (k akk −1) L akn −1) = M M . (k (k ank −1) L ann −1)
( A11k −1) = 0 ( A12k −1) , ( k −1) A22
则A(1) = L1P1AQ1.
Gauss 消 去 法
(1) (2) 设 a (1)q = max{ aij : 2 ≤ i, j ≤ n}. p
1 1
(1) 若 a p1q1 = 0 , 则消去过程结束. 否则, 交换A(1) 的第2行与第p1行, 第2列与第q1列, (1) % (1) % 记交换后的矩阵为 A(1) = A11 A12 , 则 % (1) 0 A22 % A(1) = I 2 p1 A(1) I 2 q1 = P2 A(1) Q2 .
因此计算解为x = (1.00, 1.00)T. 与精确解接近了 与精确解接近了. ②交换未知向量 的分量在方程组中的顺序 即 交换未知向量x的分量在方程组中的顺序 交换未知向量 的分量在方程组中的顺序,
1.00 0.001 x2 1.00 x = 2.00 1.00 1 3.00 1.00 0.001 1 0 ˆ 1 0 ˆ −1 ∴ L= , L = L = , U = LA = 0 1.00 −2 1 2 1
设A∈Rn×n 非奇异, 则利用列主元消去法求解线性方 程组Ax = b的步骤: (1) 用算法1.2.2计算A的列主元LU分解: PA = LU; (2) 用算法1.1.1解下三角形方程组Ly = Pb; (3) 用算法1.1.2解上三角形方程组Ux = y. 小结: 小结 列主元消去法与全主元消去法的数值稳定性差 不多, 但运算量大为减少, 是目前求解中小型线性方程组 的常用方法.
0 ˆ−1 a21 1.00 1 用算法1.1.3得l21 = = = 1000, L = =L a11 0.001 −1000 1 0 0.001 1.00 0.001 1.00 ˆ 1 所以 LA = = =U −1000 −1000 1 1.00 2.00 0
−1000 × 1.00 + 1 × 2.00 = −0.100 × 104 + 0.200 × 10 = −0.100 × 104 + 0.0002 × 104 = −0.100 × 104 = −1000
因此计算解为x = (0, 1)T. 与精确解相差甚远 与精确解相差甚远.
在例1.2.1中, 用Gauss消去法求出的结果与精确解相差 甚远, 这是由小主元引起的 解决方法: 这是由小主元引起的. (1) 用更高精度的计算机; (2) 避免出现小主元. 比如可以①交换方程组顺序 即 交换方程组顺序, 交换方程组顺序
则 PAQ = LU， (1.2.1) 其中P, Q为初等置换矩阵, 且可以证明L是单位下三角阵, 并且| lij |≤1. （书P24） (1.2.1)称为A的全主元三角分解 全主元Gauss消去法 全主元三角分解(全主元 消去法) 全主元三角分解 全主元 消去法
定理1.2.1 设A∈Rn×n, 则存在排列矩阵P, Q∈Rn×n 以及 定理 单位下三角阵L∈Rn×n 和U∈Rn×n, 使得 PAQ = LU 且| lij |≤1, U的非零对角元的个数正好等于矩阵A的秩. 算法1.2.1(全主元 全主元Gauss消去法 见P25. 消去法) 算法 全主元 消去法 一些解释： 一些解释 (1) 因需要出所有的Pi , Qi , 全主元三角分解相当于对A做 好所有的行列交换得到PAQ, 再对PAQ应用不选主元 的Gauss消去法进行三角分解. 故不能在实际中进行. (2) 全主元Gauss消去法需要大量的比较, 十分费时. 为此 有列主元 列主元Gauss消去法 即仅在列中找模最大的, 只进 消去法, 列主元 消去法 行行交换.
因此计算解为x = (1.00, 1.00)T. 与精确解接近了 与精确解接近了. 1.2.2 全主元消去法
y1 1.00 1 0 y1 1.00 y = ⇒ y = 2 1 2 3.00 1.00 2 x1 = 1.00 1.00 0.001 x2 1.00 x = ⇒ x2 = 1.00 − 0.001 × 1.00 1.00 1 1.00 0 = 0.100 × 10 − 0.0001 × 10 = 1.00
例1.2.1 在3位10进制的浮点数系下解方程组
0.001 1.00 x1 1.00 x = 1.00 2.00 2 3.00 −1 x1 0.001 1.00 1.00 1.002... 解: 精确解 = = x2 1.00 2.00 3.00 0.998...
0 0
若a
(0) p0 q0
= 0 , 则结束.
否则, 的第1行与第p0行, 第1列与第q0列, % 记交换后的矩阵为 A(0) , 则
% A(0) = I1 p0 A(0) I1q0 = P A(0) Q1 . 1
% A(0)
(1)
交换A(0)
选 主 元
Gauss
% = L1 A
(0)
换L1
% (0) %n a12 L a1(0) (1) (1) (1) a22 L a2 n A11 = M O M 0 (1) (1) an 2 L ann
(1) A12 , (1) A22
A
% (0) a11 0 = M 0
Gauss
换Lk
Gauss Gauss
(4) 设全主元Gauss消去过程进行到第r步终止, 则有 LrPr···L1P1AQ1···Qr = U 其中U为上三角阵, P1 , P2 , ··· , Pr和Q1 , Q2 , ··· , Qr为初等 置换矩阵, L1 , L2 , ··· , Lr为单位下三角阵, 且| lij |≤1. 令 Q = Q1Q2···Qr P = PrPr-1···P1 L = P(LrPr···L2P2L1P1)－1,
§1.2 选主元三角分解
1.2.1 选主元的重要性 一些疑问： 一些疑问： (1) A非奇异可保证Ax = b有唯一解，但A非奇异并不能 保证其顺序主子阵均非奇异，因此不能保证Gauss消去 过程能够进行到底, 怎样修正算法1.1.3？ (2) Gauss消去法中计算lik时位于分母上的主元虽不为0 但很小时, 是否会对算法产生不良影响？怎样修正该算 法？
1.00 2.00 x1 3.00 x = 0.001 1.00 2 1.00 a21 0.001 用算法1.1.3得 l21 = = = 0.001, a11 1.00 0 ˆ 0 1 1 −1 L= , L = L = −0.001 1 0.001 1 0 1.00 2.00 1.00 2.00 ˆ = LA = 1 U = −0.001 1 0.001 1.00 0 1.00
% A( k −1) = I 2 pk A( k −1) I 2 qk
% A( k −1)
( A11k −1) = Pk A( k −1) Qk = 0
% A1(2k −1) . ( k −1) % A22
选 主 元