三角分解法解线性方程组

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;}}} 运⾏截图： ⾄此完毕。

线性方程组AX=B 的数值计算方法实验一、 实验描述：随着科学技术的发展，线性代数作为高等数学的一个重要组成部分，在科学实践中得到广泛的应用。

本实验的通过C 语言的算法设计以及编程，来实现高斯消元法、三角分解法和解线性方程组的迭代法（雅可比迭代法和高斯-赛德尔迭代法），对指定方程组进行求解。

二、 实验原理：1、高斯消去法：运用高斯消去法解方程组，通常会用到初等变换，以此来得到与原系数矩阵等价的系数矩阵，达到消元的目的。

初等变换有三种：(a)、(交换变换)对调方程组两行；(b)、用非零常数乘以方程组的某一行；(c)、将方程组的某一行乘以一个非零常数，再加到另一行。

通常利用(c)，即用一个方程乘以一个常数，再减去另一个方程来置换另一个方程。

在方程组的增广矩阵中用类似的变换，可以化简系数矩阵，求出其中一个解，然后利用回代法，就可以解出所有的解。

2、选主元：若在解方程组过程中，系数矩阵上的对角元素为零的话，会导致解出的结果不正确。

所以在解方程组过程中要避免此种情况的出现，这就需要选择行的判定条件。

经过行变换，使矩阵对角元素均不为零。

这个过程称为选主元。

选主元分平凡选主元和偏序选主元两种。

平凡选主元：如果()0p pp a ≠，不交换行；如果()0p pp a =，寻找第p 行下满足()0p pp a ≠的第一行，设行数为k ，然后交换第k 行和第p 行。

这样新主元就是非零主元。

偏序选主元：为了减小误差的传播，偏序选主元策略首先检查位于主对角线或主对角线下方第p 列的所有元素，确定行k ，它的元素绝对值最大。

然后如果k p >，则交换第k 行和第p 行。

通常用偏序选主元，可以减小计算误差。

3、三角分解法：由于求解上三角或下三角线性方程组很容易所以在解线性方程组时，可将系数矩阵分解为下三角矩阵和上三角矩阵。

其中下三角矩阵的主对角线为1，上三角矩阵的对角线元素非零。

有如下定理：如果非奇异矩阵A 可表示为下三角矩阵L 和上三角矩阵U 的乘积： A LU = (1) 则A 存在一个三角分解。

三角线性方程组的数值计算1. 实验描述有回代的高斯消元法：如果A 是N ×N 非奇异矩阵，则存在线性方程组UX=Y 与线性方程组AX=B 等价，这里U 是上三角矩阵，并且u kk ≠0。

当构造出U 和Y 后，可以用回代法求解UX=Y ，并且得到方程组的解X 。

三角分解法：设线性方程组AX=B 的系数矩阵A 存在三角分解（如果非奇异矩阵A 可以表示为下三角矩阵L 和上三角矩阵U 的乘积：A=LU 则A 存在一个三角分解），则线性方程组可以表示为 LUX=B 而方程组的解可以通过定义Y=UX 并求解下面的俩个方程组得到。

首先对方程组LY=B 求解Y ，然后对方程组UX=Y 求解X 。

雅可比迭代：设矩阵A 具有严格对角优势，则AX=B 有唯一解X=P 利用迭代式（x j（k+1)=b j −a j 1x k −⋯−a jj −1x j −1 k −a jj +1x jj +1 k−⋯a jN x N(k )a jj(其中j=1，2，…，N)）可产生一个向量序列﹛P k ﹜，而且对于任意初始向量P 0，向量序列都将收敛到P 0。

高斯-赛德尔迭代：x j（k+1)=b j −a j 1x k +1 −⋯−a jj −1x j −1 k +1 −a jj +1x jj +1 k−⋯−a jN x N(k )a jj2. 实验内容设有如下三角线性方程组，而且系数矩阵具有严格对角优势：d 1x 1+c 1x 1 =b 1a 1+d 2x 2+c 2x 2=b 2a 2x 2+d 3x 3+c 3x 4 =b 3 . . . . . . . . . . . . a N-2x N-2+d N-1x N-1+c N-1x N =b N-1a N-1x N-1+d N x N =b N假设d1=d2=…=dn=12，b1=b2=…=bn=5;c1=c2=..=cn=a1=a2=…=an=-2; 1.用高斯消元法计算；2.采用LU 分解来计算；3.采用雅可比迭代法计算；4.采用高斯-赛德尔迭代法计算。

三角分解法解线性方程组#include<iostream.h>#include<iomanip.h>#include<stdlib.h>//----------------------------------------------全局变量定义区 const int Number=15; //方程最大个数doublea[Number][Number],b[Number],copy_a[Number][Number],copy_b[Number]; // 系数行列式int A_y[Number]; //a[][]中随着横坐标增加列坐标的排列顺序,如a[0][0],a[1][2],a[2][1]...则A_y[]={0,2,1...}; int lenth,copy_lenth;//方程的个数char * x; //未知量a,b,c的载体int i,j;//----------------------------------------------函数声明区 voidinput(); //输入方程组void print_menu(); //打印主菜单int Doolittle_check(double a[][Number],double b[Number]); //判断是否行列式>0,若是,调整为顺序主子式全>0void xiaoqu_u_l(); //将行列式Doolittle分解 void calculate_u_l(); //计算Doolittle结果 void exchange(int m,int i); //交换A_y[m],A_y[i] void exchange_lie(int j); //交换a[][j]与b[]; void exchange_hang(int m,int n);//分别交换a[][]和b[]中的m与n两行 void exchange_a_lie(int m,int n); //交换a[][]中的m和n列 void exchange_x(int m,int n); //交换x[]中的x[m]和x[n] //函数定义区void print_menu(){system("cls");cout<<"------------方程系数和常数矩阵表示如下:\n"; for(intj=0;j<lenth;j++)cout<<"系数"<<j+1<<" ";cout<<"\t常数";cout<<endl;for(int i=0;i<lenth;i++){for(j=0;j<lenth;j++)cout<<setw(8)<<setiosflags(ios::left)<<a[i][j];cout<<"\t"<<b[i]<<endl; }}void input(){int i,j;cout<<"方程的个数:";cin>>lenth;if(lenth>Number){cout<<"It is too big.\n";return;}x=new char[lenth];for(i=0;i<lenth;i++)x[i]='a'+i;//输入方程矩阵//提示如何输入cout<<"====================================================\n";cout<<"请在每个方程里输入"<<lenth<<"系数和一个常数:\n"; //输入每个方程for(i=0;i<lenth;i++){cout<<"输入方程"<<i+1<<":";for(j=0;j<lenth;j++)cin>>a[i][j];cin>>b[i];}}void Doolittle() //Doolittle消去法计算方程组 {double temp_a[Number][Number],temp_b[Number];int i,j,flag;for(i=0;i<lenth;i++)for(j=0;j<lenth;j++)temp_a[i][j]=a[i][j]; flag=Doolittle_check(temp_a,temp_b);if(flag==0) cout<<"\n行列式为零.无法用Doolittle求解."; xiaoqu_u_l();calculate_u_l();cout<<"用Doolittle方法求得结果如下:\n";for(i=0;i<lenth;i++) //输出结果{for(j=0;x[j]!='a'+i&&j<lenth;j++);cout<<x[j]<<"="<<b[j]<<endl;}}void calculate_u_l() //计算Doolittle结果{ int i,j;double sum_ax=0; for(i=0;i<lenth;i++){for(j=0,sum_ax=0;j<i;j++)sum_ax+=a[i][j]*b[j];b[i]=b[i]-sum_ax;}for(i=lenth-1;i>=0;i--){for(j=i+1,sum_ax=0;j<lenth;j++)sum_ax+=a[i][j]*b[j];b[i]=(b[i]-sum_ax)/a[i][i]; }}void xiaoqu_u_l() //将行列式按Doolittle分解{ int i,j,n,k;double temp; for(i=1,j=0;i<lenth;i++)a[i][j]=a[i][j]/a[0][0]; for(n=1;n<lenth;n++){ //求第n+1层的上三角矩阵部分即Ufor(j=n;j<lenth;j++){ for(k=0,temp=0;k<n;k++)temp+=a[n][k]*a[k][j];a[n][j]-=temp;}for(i=n+1;i<lenth;i++) //求第n+1层的下三角矩阵部分即L{ for(k=0,temp=0;k<n;k++)temp+=a[i][k]*a[k][n];a[i][n]=(a[i][n]-temp)/a[n][n];}}}int Doolittle_check(double temp_a[][Number],double temp_b[Number]) //若行列式不为零,将系数矩阵调整为顺序主子式大于零{int i,j,k,maxi;double lik,temp;for(k=0;k<lenth-1;k++){j=k;for(maxi=i=k;i<lenth;i++)if(temp_a[i][j]>temp_a[maxi][j]) maxi=i;if(maxi!=k){ exchange_hang(k,maxi);for(j=0;j<lenth;j++){ temp=temp_a[k][j];temp_a[k][j]=temp_a[maxi][j];temp_a[maxi][j]=temp;}}for(i=k+1;i<lenth;i++){lik=temp_a[i][k]/temp_a[k][k];for(j=k;j<lenth;j++)temp_a[i][j]=temp_a[i][j]-temp_a[k][j]*lik;temp_b[i]=temp_b[i]-temp_b[k]*lik;}}if(temp_a[lenth-1][lenth-1]==0) return 0;return 1;}void exchange_hang(int m,int n) //交换a[][]中和b[]两行 { int j; double temp;for(j=0;j<lenth;j++){ temp=a[m][j];a[m][j]=a[n][j];a[n][j]=temp;}temp=b[m];b[m]=b[n];b[n]=temp;}void exchange(int m,int i) //交换A_y[m],A_y[i] { int temp;temp=A_y[m];A_y[m]=A_y[i];A_y[i]=temp;}void exchange_lie(int j) //交换未知量b[]和第i列 { double temp;int i; for(i=0;i<lenth;i++){ temp=a[i][j];a[i][j]=b[i];b[i]=temp;}}void exchange_a_lie(int m,int n) //交换a[]中的两列{ double temp;int i; for(i=0;i<lenth;i++) { temp=a[i][m];a[i][m]=a[i][n];a[i][n]=temp;}}void exchange_x(int m,int n) //交换未知量x[m]与x[n] { char temp;temp=x[m];x[m]=x[n];x[n]=temp;}//主函数void main(){int flag=1;input(); //输入方程while(flag){print_menu(); //打印主菜单cout<<"用Doolittle方法求得结果如下:\n";for(i=0;i<lenth;i++) //输出结果{for(j=0;x[j]!='a'+i&&j<lenth;j++);cout<<x[j]<<"="<<b[j]<<endl;}}}。

计算方法实验报告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的语句还不是非常熟悉,因此在编程过程中,出现了许多错误提示;并且此次编程的两种方法对矩阵的运算也比较复杂;问题主要集中在循环控制中,循环次数多了一次或者缺少了一次,导致数据错误,一些基本的编程语句在语法上也会由于生疏而产生许多问题,但是语句的错误由于系统会提示,比较容易进行修改,数据计算过程中的一些逻辑错误,比如循环变量的控制,这些系统不会提示错误,需要我们细心去发现错误,不断修正,调试;。

Cholesky分解法Cholesky分解法又称三角分解法，或称因子化法设线性方程组Ax b=（1）式中A为对称、正定的矩阵。

对于对称、正定的矩阵A，可进行分解TA LDL=（2）式中L是下单位三角阵，D是对角线矩阵。

右端项列向量（列阵）也作相应的分解b LDb=（3）将式（2）和式（3）代入方程（1），得到上三角方程组TL x b=再按诸如高斯消元法的回代过程就可解出x下面，通过编写一个完整的C语言程序求1234121312235416155118341720xxxx⎛⎫⎛⎫⎛⎫⎪⎪ ⎪- ⎪⎪ ⎪=⎪⎪ ⎪--⎪⎪ ⎪-⎝⎭⎝⎭⎝⎭，完整的C语言程序代码如下：#include<>#include<>#include<>#define N 4void Cholesky(int n,double A[N][N],double x[N],double b[N]) {int i,j,k;double L[N][N],D[N][N],b2[N];i=1;D[i-1][i-1]=A[i-1][i-1];for(i=2;i<=n;i++){j=1;L[i-1][j-1]=A[i-1][j-1]/D[j-1][j-1];if(j==i-1){D[i-1][i-1]=0;for(k=1;k<=i-1;k++)D[i-1][i-1]+=pow(L[i-1][k-1],2)*D[k-1][k-1];D[i-1][i-1]=A[i-1][i-1]-D[i-1][i-1];continue;}else{for(j=2;j<=i-1;j++){L[i-1][j-1]=0;for(k=1;k<=j-1;k++)L[i-1][j-1]+=L[i-1][k-1]*L[j-1][k-1]*D[k-1][k-1];L[i-1][j-1]=(A[i-1][j-1]-L[i-1][j-1])/D[j-1][j-1];if(j==i-1){D[i-1][i-1]=0;for(k=1;k<=i-1;k++)D[i-1][i-1]+=pow(L[i-1][k-1],2)*D[k-1][k-1];D[i-1][i-1]=A[i-1][i-1]-D[i-1][i-1];continue;}}}}i=1;b2[i-1]=b[i-1]/D[i-1][i-1];for(i=2;i<=n;i++){b2[i-1]=0;for(k=1;k<=i-1;k++)b2[i-1]+=L[i-1][k-1]*D[k-1][k-1]*b2[k-1];b2[i-1]=(b[i-1]-b2[i-1])/D[i-1][i-1];}x[n-1]=b2[n-1] ;for(i=n-1;i>=1;i--){x[i-1]=;for(k=i+1;k<=n;k++)x[i-1]+=L[k-1][i-1]*x[k-1];x[i-1]=b2[i-1]-x[i-1];}for(i=0;i<=n-1;i++)printf("%f\n",x[i]);}void main(){double A[4][4]={{1,2,1,3},{2,3,-5,4},{1,-5,5,-1},{3,4,-1,7}}; double x[4]={0};double b[4]={12,16,18,20};Cholesky(4,A,x,b);}程序运行结果：。

第2章线性方程组的解法--------学习小结本章学习体会本章主要学习的是线性方程组的解法。

而我们则主要学习了高斯消去法、直接三角分解法以及迭代法三种方法。

这三种方法的优缺点以及适用范围各有不同。

高斯消去法中，我们又学习了顺序高斯消去法以及列主元素高斯消去法。

顺序高斯消去法可以得到方程组的精确解，但要求系数矩阵的主对角线元素不为零，而且该方法的数值稳定性没有保证。

但列主元素高斯消去法因为方程顺序的调整，其有较好的数值稳定性。

直接三角分解法中，我们主要学习了Doolitte分解法与Crout分解法。

其思想主要是：令系数矩阵A=UL，其中L为下三角矩阵，U是上三角矩阵，为求AX=b 的解，则引进Ly=b,Ux=y 两个方程，以求X得解向量。

这种方法计算量较小，但是条件苛刻，且不具有数值稳定性。

迭代法（逐次逼近法）是从一个初始向量出发，按照一定的计算格式，构造一个向量的无穷序列，其极限才是所求问题的精确解，只经过有限次运算得不到精确解。

该方法要求迭代收敛，而且只经过有限次迭代，减少了运算次数，但是该方法无法得到方程组的精确解。

二、本章知识梳理针对解线性方程组，求解线性方程组的方法可分为两大类：直接法和迭代法，直接法（精确法）：指在没有舍入误差的情况下经过有限次运算就能得到精确解。

迭代法（逐次逼近法）：从一个初始向量出发，按照一定的计算格式，构造一个向量的无穷序列，其极限才是所求问题的精确解，只经过有限次运算得不到精确解。

我们以前用的是克莱姆法则，对于计算机来说，这种方法运算量比较大，因此我们学习了几种减少运算次数的方法，有高斯消去法、直接三角分解法，同时针对病态方程组，也提出了几种不同的解法。

Gauss消去法Gauss消去法由消元和回代两个过程组成，消元过程是指针对方程组的增广矩阵，做有限次初等行变化，使它系数矩阵变为上三角矩阵。

顺序Gauss消去法消元过程：对于K=1,2,3…,n-1执行如果,则算法失效，停止计算；否则转（2）对于计算回代过程：综上：顺序Gauss消去法的数值稳定性是没有保证的。

线性代数方程组的解法关键词：线性代数方程组；高斯消元法；列主元消元法；三角分解法；杜立特尔分解法；迭代法；雅可比迭代法；高斯-赛德尔迭代法1引言目前，解线性代数方程组在计算机上常用的的方法大致把它分为两类：“直接法”与“迭代法”.在线性代数中曾指出阶线性代数方程组有唯一的解，并且可以用克拉默法则求方程组的解，初次看来问题已经解决，但从使用效果看并不是这样的.因为求阶线性代数方程组，如果用克拉默法则，需要计算个阶行列式，每个阶行列式为项之和，每项又是个元素的乘积，所以计算中仅乘法次数就高达次，当较大时，它的计算量是非常惊人的.因为现在所碰到的很多问题都需要很大的计算量，故需要好用的算法来求解.先来回顾一下回代过程和迭代过程.(1)是一个三角形方程组，当有唯一解时，可以用反推的方式求解，也就是先从第个方程解得， (2)然后代入第个方程，可得到， (3)如此继续下去，假设已得到，， , ,代进第个方程即得的计算， (4)上述求解的过程叫做回代过程.定义1[1] (向量的范数) 若向量的某个实值函数满足1.是非负的，即且的充要条件是 ;2.是齐次的，即 ;3.三角不等式，即对，总是有.那么上向量的范数（或模）就是 .下面给几个最常遇到的向量范数.向量的“1”范数：(5)向量的“2”范数：(6)向量的范数：(7)例1设求 , , .解由式（5）,（6）及（7）知.定义2若矩阵的某个实值函数满足1.是非负的，即且的充要条件是 ;2.是齐次的，即 ;3.三角不等式，即对总有;1.矩阵的乘法不等式，即对总有，那么称为上矩阵的范数（或模）.表 1是矩阵几个常用算子范数的定义与算式.表 1范数名称记号定义计算公式“1”范数（又名列模）“2”范数（又名谱模）“”范数（又名行模）的极限就是方程组的解向量，这时候在给定允许的误差内，只要适当的大,就可以作为方程组在满足精度要求条件下的近似解.这种求近似解的方法就是解线性方程组的一类基本的迭代解法，其中称为迭代矩阵，公式（9）称迭代公式（或迭代过程）,由迭代公式得到的序列叫做迭代序列.如果迭代的序列是收敛的，则称为迭代法收敛；如果迭代的序列是不收敛，则称它是迭代法发散.定理3设 .如果约化主元素，则可以利用高斯消元的方法把方程组约化成三角形方程组来求解，其计算公式如下：（1）消元计算：对依次计算（2）回代计算：3用高斯消元法与列主元消元法解线性代数方程组（重点）！3.1 高斯消元法解方程组用高斯消元的方法求线性代数方程组的解的整个计算过程可分为两个环节，也就是利用按照次序消去未知数的方法，把原来的方程组转化成跟它同解的三角形方程组（这个转化的过程叫消元过程），再通过回代过程求三角形方程组的解，最终得到原来方程组的解.其中按照方程的顺进行消元的高斯消元法，又叫顺序消元法.3.2列主元消元法解方程组列主元消元法实际上是一种行交换的消元法，它跟顺序消元法比较而言，主要特点是在进行第次消元前，不管的值是否等于零，都在子块的第一列中选择一个元，使，并将中的第行元与第行元互相变换（相当于交换同解方程组中的第个方程），然后再进行消元计算得到结果.注：列主元素法的精度虽然稍低于全主元素法[1]，但它计算简单，相对比全主元素法它的工作的量大大减少，并且从计算经验和理论分析都可以表明，它与全主元素法同样拥有很好的值稳定性，列主元素法是求解中小型浓密型方程组的最好的方法之一.4用三角分解法解线性代数方程组4.1 矩阵的三角分解定义4把一个阶矩阵分解成两个三角矩阵相乘的形式称为矩阵的三角分解.常见的矩阵三角分解是其中是下三角形的矩阵,是上三角形的矩阵.定理5[1]（矩阵三角分解基本定理）设 .若的顺序主子式，那么存在唯一的杜利特尔分解其中是单位下三角形矩阵，为非奇异的上三角形矩阵.如果是单位下三角形的矩阵，是上三角形的矩阵，那么把这种分解法称为杜利特尔分解法，其中杜利特尔分解法是这种三角分解的一种特例，下面主要介绍利用杜利特尔分解法来求方程组的解.4.2 用杜利特尔分解法解线性代数方程组用杜利特尔分解法解方程组的步骤可以把它归纳为（1）实现分解，也就是1.按算式（11）（12）依次计算的第一行元与的第一列元；1.对按算式（13）（14）依次计算的第行元与的第列元.（2）求解三角形方程组，即按算式依次计算 .（3）求解三角形方程组，即按算式依次计算.利用杜利特尔分解法解方程组与高斯消元法是相似的，它重要的优点是：在利用分解，解有相同的系数矩阵的方程组时,用杜利特尔分解法非常方便，只用两个式子就可以得到方程组的解.5用迭代法解线性代数方程组用迭代法求方程组的解，需要考虑迭代过程的收敛性，在下面的讨论中，都假设方程组的系数矩阵的对角阵是不为零的.5.1 用雅可比迭代法解方程组对于一般线性方程组，如果从第个方程解出，就可以把它转化成等价的方程组. （15）从而可以得到对应的迭代公式（16）这就是解一般方程组的分量形式的雅可比(Jacobi)迭代公式.如果把它改成（17）并把系数矩阵表示成（18）其中则可以看出式的左右两端分别是向量和的第个分量，故因为可逆，所以于是就可以得到是雅可比迭代的公式.其中（称为雅可比迭代矩阵）， .5.2 用高斯-赛德尔迭代法解方程组高斯-赛德尔迭代法也是常用的迭代法，设线性代数方程组为，则高斯-赛德尔迭代法的迭代公式为（19）其中迭代法（19）就称为高斯-赛德尔迭代法.通过雅可比迭代法类似的途径，就可以得到矩阵的表达式其中（称为高斯-赛德尔迭代矩阵）， .高斯-赛德尔迭代法与雅可比迭代法都有算式简单、容易在计算机上实现等优点，但是用计算机来计算时，雅可比迭代法需要两组工作单元用来寄存与的量，而高斯赛-德尔迭代法只需一组工作单元存放或的分量.对于给定的线性方程组，用这两种方法求解可能都收敛或者都不收敛，也可能一个收敛另一个不收敛，两种方法的收敛速度也不一样.5.3 迭代法的收敛条件与误差分析定义6[1]矩阵全部的特征值的模的最大值，叫做矩阵的谱半径，记作 ,即.定理7[1]对任意初始向量迭代过程收敛的充要条件是；当时，越小，那么其收敛的速度是越快的.由定理7可知，用雅可比迭代法求解时，其迭代的过程是收敛的，而用高斯-赛德尔迭代法来求解,其迭代的过程是发散的.在不同条件下，收敛的速度是不同的，对同一矩阵，一种方法是收敛的，一种方法发散.第 7 页。