选主元的三角分解法
- 格式:doc
- 大小:469.00 KB
- 文档页数:8
第二章 线性方程组的数值解法在科技、工程技术、社会经济等各个领域中很多问题常常归结到求解线性方程组。
例如电学中的网络问题,样条函数问题,构造求解微分方程的差分格式和工程力学中用有限元方法解连续介质力学问题,以及经济学中求解投入产出模型等都导致求解线性方程组。
n 阶线性方程组的一般形式为⎪⎪⎩⎪⎪⎨⎧=+++=+++=+++nn nn n n n n n n b x a x a x a b x a x a x a b x a x a x a L K K K K L L 22112222212********* (1.1) 其矩阵形式为b Ax = (1.2) 其中⎥⎥⎥⎥⎦⎤⎢⎢⎢⎢⎣⎡=⎥⎥⎥⎥⎦⎤⎢⎢⎢⎢⎣⎡=⎥⎥⎥⎥⎦⎤⎢⎢⎢⎢⎣⎡=n n nn n n n n b b b b x x x x a a a a a a a a a A M M L K K K K L L 2121212222111211),,2,1,(n j i a ij L =,),,2,1(n i b i L =均为实数,i b 不全为0,且A 为非奇异。
关于线性方程组的数值解法一般分为两类:1.直接法 就是不考虑计算机过程中的舍入误差时,经有限次的四则运算得到方程组准确解的方法。
而实际中由于计算机字长的限制,舍入误差的存在和影响,这种算法也只能求得线性方程组的近似解。
本章将阐述这类算法中最基本的消去法及其某些变形。
这些方法主要用于求解低阶稠密系数矩阵方程组。
2.迭代法 从某个解的近似值出发,通过构造一个无穷序列,用某种极限过程去逐步逼近线性方程组的精确解的方法。
本章主要介绍迭代法与迭代法。
迭代法是解大型稀疏矩阵(矩阵阶数高而且零元素较多)的线性方程组的重要方法。
§1 高斯)(Gauss 消去法1.1 Gauss 消去法Gauss 消去法是将线性方程组化成等价的三角形方程组求解。
首先举例说明Gauss消去法的基本思想和过程。
初中数学因式分解的几种经典技巧初中数学因式分解的几种经典方法因式分解是初中数学的一个重点,涉及到分式方程和一元二次方程,因此学会一些基本的因式分解方法非常必要。
下面列举了九种方法,希望对大家的研究有所帮助。
1.提取公因式这种方法比较常规、简单,必须掌握。
常用的公式有完全平方公式、平方差公式等。
例如,对于方程2x-3x=0,可以进行如下因式分解:x(2x-3)=0,得到x=0或x=3/2.一个规律是:当一个方程有一个解x=a时,该式分解后必有一个(x-a)因式,这对我们后面的研究有帮助。
2.公式法将式子利用公式来分解,也是比较简单的方法。
常用的公式有完全平方公式、平方差公式等。
建议在使用公式法前先提取公因式。
例如,对于x^2-4,可以使用平方差公式得到(x+2)(x-2)。
3.十字相乘法是做竞赛题的基本方法,但掌握了这个方法后,做平时的题目也会很轻松。
关键是将二次项系数a分解成两个因数a1和a2的积a1.a2,将常数项c分解成两个因数c1和c2的积c1.c2,并使ac正好是一次项b,那么可以直接写成结果。
例如,对于2x^2-7x+3,可以使用十字相乘法得到(x-3)(2x-1)。
总结:对于二次三项式ax^2+bx+c(a≠0),如果二次项系数a可以分解成两个因数之积,即a=a1.a2,常数项c可以分解成两个因数之积,即c=c1.c2,那么可以使用十字相乘法进行因式分解。
文章中有一些格式错误,需要修正。
另外,第四段中的一些内容似乎有问题,建议删除。
改写后的文章如下:分解因式是数学中的一个重要概念,也是许多数学问题的基础。
在中学数学中,我们通常研究到七种分解因式的方法。
1.公因数法这种方法是最基础的方法之一,它的核心思想是找到表达式中的公因数。
例如,对于表达式6x+9y,我们可以找到它们的公因数3,然后将表达式简化为3(2x+3y)。
2.公式法公式法是通过运用数学公式来分解因式。
例如,对于二次三项式ax2+bx+c,我们可以使用求根公式来求出它的两个根,然后将表达式分解为(a(x-根1)(x-根2))的形式。
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;}}} 运⾏截图: ⾄此完毕。
计算方法实验报告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的语句还不是非常熟悉,因此在编程过程中,出现了许多错误提示;并且此次编程的两种方法对矩阵的运算也比较复杂;问题主要集中在循环控制中,循环次数多了一次或者缺少了一次,导致数据错误,一些基本的编程语句在语法上也会由于生疏而产生许多问题,但是语句的错误由于系统会提示,比较容易进行修改,数据计算过程中的一些逻辑错误,比如循环变量的控制,这些系统不会提示错误,需要我们细心去发现错误,不断修正,调试;。
列主元三角分解法步骤嘿,朋友们!今天咱来唠唠列主元三角分解法那些事儿。
你说这列主元三角分解法啊,就好比是搭积木,咱得一层一层稳稳地搭起来。
咱先说说啥是列主元。
这就好像是在一群数字小伙伴里,挑出那个最厉害、最能带头的数字来!它就是那个主元啦。
然后呢,咱就开始分解啦。
这就像是把一个大难题,一点点地拆解成小步骤,让它变得不再那么可怕。
你看啊,我们要通过巧妙的变换,让矩阵变得有秩序,就像整理房间一样,把东西都归置得整整齐齐。
比如说,咱有个矩阵,里面的数字就像一群调皮的小孩子,到处乱跑。
那我们就得想办法让它们听话呀!找到那个主元,把它放在合适的位置,然后其他数字就会跟着它的节奏来。
这过程中,可不能马虎哦!就像走路一样,一步一步都得走稳了。
要是有一步走错了,那可能就前功尽弃啦。
你想想,如果我们不仔细地找主元,随便乱来,那结果能对吗?那肯定不行呀!而且哦,这列主元三角分解法用处可大啦!在解决很多数学问题的时候,它就像是一把神奇的钥匙,能打开那扇困住我们的门。
它能让复杂的计算变得简单起来,就像给我们加了一双翅膀,能在数学的天空中自由翱翔。
咱再打个比方,它就像是一个好向导,带着我们在数学的迷宫里找到正确的路。
这可不是随便说说的呀!很多厉害的数学家都在用它呢!他们用它解决了一个又一个难题,创造了一个又一个奇迹。
咱虽然不是数学家,但咱也可以学会呀!只要我们有耐心,肯钻研,就一定能掌握这神奇的列主元三角分解法。
所以啊,朋友们,别害怕数学,别害怕这列主元三角分解法。
它其实没那么难,只要我们用心去学,去理解,就一定能学好它。
让我们一起加油,在数学的海洋里畅游吧!就这么定啦!。
线性方程组求解的直接法5.2线性方程组直接解法概述直接解法就是利用一系列公式进行有限步计算,直接得到方程组的精确解的方法.当然,实际计算结果仍有误差,譬如舍入误差,而且舍入误差的积累有时甚至会严重影响解的精度.这是一个众所周知的古老方法,但用在计算机上仍然十分有效.求解线性方程组最基本的一种直接法是消去法.消去法的基本思想是,通过将一个方程乘以或除以某个常数,以及将两个方程相加减这两种手段,逐步减少方程中的变元的数目,最终使每个方程仅含一个变元,从而得出所求的解.高斯(Gauss )消去法是其中广泛应用的方法,其求解过程分为消元过程和回代过程两个环节.消元过程将所给的方程组加工成上三角方程组,所归结的方程组再通过回代过程得出它的解.Gauss 消去法由于添加了回代的过程,算法结构稍复杂,但这种改进的算法明显减少了计算量.直接法比较适用于中小型方程组.对高阶方程组,即使系数矩阵是稀疏的,但在运算中很难保持稀疏性,因而有存储量大,程序复杂等不足.5.3直接解法5.3.1Gauss 消去法Gauss 消去法是一个古老的求解线性方程组的方法,由它改进而来的选主元法是目前计算机上常用的有效的求解低阶稠密矩阵线性方程组的方法.例5.1用Gauss 消去法解方程组1231231232221(5.3.1)1324 (5.3.2)2539(5.3.3)2x x x x x x x x x ⎧++=⎪⎪++=⎨⎪++=⎪⎩解〖JP4〗第1步,式35.3.12⨯-()()加到式(5.3.2)上,式()15.3.1()2⨯-加到式(5.3.3)上,得到等价方程组123232322211(5.4.4)282(5.4.5)x x x x x x x ⎧++=⎪⎪-+=-⎨⎪⎪+=⎩第2步,式()2⨯5.3.4加到式(5.3.5)上得等价的方程组12323322211100(5.3.6)x x x x x x ++=⎧⎪-+=-⎨⎪=⎩第3步,回代法求解方程组(5.3.6),即可求得该方程组的解为32110,1,.2x x x ===-.用矩阵描述其约化过程即为233(2)22221011100100r r r ⨯+⇒⎡⎤⎢⎥--⎢⎥⎢⎥⎣⎦→[]122133(1)3()21()222212221,3241/201111395/20282r r r r r r A b ⨯-+⇒⨯-+⇒⎡⎤⎡⎤⎢⎥⎢⎥=--⎢⎥⎢⎥⎢⎥⎢⎥⎣⎦⎣⎦→.这种求解过程称为具有回代的Gauss 消去法.由此例可见,Gauss 消去法的基本思想是:用矩阵的初等行变换将系数矩阵A 化为具有简单形式的矩阵(如上三角阵、单位矩阵等),而三角形方程组是很容易回代求解的.一般地,设有n 个未知数的线性方程组为11112211211222221122n n n n n n nn n na x a x a xb a x a x a x b a x a x a x b +++=⎧⎪+++=⎪⎨⎪⎪++=⎩L L MM M L (5.3.7)1212)(,,)(,,)T T ij n n n n A a X x x x b b b b ⨯===L L (,,,则方程组(5.3.7)化为AX b =.方便起见,记()(1)det 0A AA ==≠,(1)b b =,且()1A的元素记为()()11,ij a b ,的元素记为()1i b ,则消去法的步骤如下:第1步:1110a≠(),,计算(1)11(1)11(2,3,4),i i a m i n a ==L 用()1i m -乘方程组(5.3.7)中的第1个方程加到第i个方程中()2,3,i n =L ,即进行行初等变换()112,3,i i i R m R R i n -⋅→=L ,消去第2个到第n个方程中的未知数1,x ,得等价方程组111121(2)(2)(2)22222(2)(2)(2)2inn n n nn n x a a b x a a b ⎢⎥⎢⎥⎢⎥⎢⎥⎢⎥⎢⎥=⎢⎥⎢⎥⎢⎥⎢⎥⎢⎥⎢⎥⎢⎥⎢⎥⎣⎦⎣⎦⎣⎦LMM LM M L (5.3.8)记为(2)(2)A X b =,其中(2)(1)(1)(2)(1)(1)1111(,2,3),2,3,ij ij i j i i i a a m a i j n b b m b i n =-==-=L L ,,第k 步()1,2,1k n =-L:继续上述消元过程.第1步到第1k -步计算已完成,且得到与原方程组等价的方程组(1)(1)(1)(1)1112111(2)(2)(2)222223()()()()()()nn k k k kkkn k n k k k nk nn n a a a b x a a b xx aa b x a a b ⎡⎤⎡⎤⎡⎤⎢⎥⎢⎥⎢⎥⎢⎥⎢⎥⎢⎥⎢⎥⎢⎥⎢⎥=⎢⎥⎢⎥⎢⎥⎢⎥⎢⎥⎢⎥⎢⎥⎢⎥⎢⎥⎢⎥⎢⎥⎣⎦⎢⎥⎢⎥⎣⎦⎣⎦L L LLOM L M MMM L(5.3.9)记为()(()K k A X b =,进行第k 步消元:设()0k kka≠,计算乘数()()(1,)k ikk ik kka m k k n a ==+L ,用ik m -乘方程组(5.3.9)中第k 个方程加到第i 1)i k n =+L (,,,个方程上消去方程组(5.3.9)中第i 1)i k n =+L (,,个方程的未知数k x ,得到与原方程组等价的方程组:(1)()()(1)()()(1)(1)()(,1,)( 1.)k k k ij ij ik kj k k k i i ik k k k k k a a m a i j k n b b m b i k n A A k b b k ++++⎧=-=+⎪=-=+⎨⎪⎩L L ()与前行元素相同,与前个元素相同 (5.3.10) 记为(1)(1)k k A X b ++=其中(1)(1,k k A b ++)中元素计算公式为(1)()()(1)()()(1)(1)()(,1,)( 1.)k k k ij ij ik kj k k k i i ik k k k k k a a m a i j k n b b m b i k n A A k b b k ++++⎧=-=+⎪=-=+⎨⎪⎩L L ()与前行元素相同,与前个元素相同 (5.3.11)重复上述过程,且设()0(1,2,1)k kk a k n ≠=-L ,共完成1n -步消元计算,得到与方程组(5.3.7)等价的三角形方程组1111211(2)(2)(2)22222()()n n n n n nn n x a a b x a b ⎢⎥⎢⎥⎢⎥⎢⎥⎢⎥⎢⎥=⎢⎥⎢⎥⎢⎥⎢⎥⎢⎥⎢⎥⎢⎥⎢⎥⎣⎦⎣⎦⎣⎦LMOM M (5.3.12)再用回代法求方程组(5.3.12)的解,计算公式为()()()()1()(),(1,2,1)n n n nn n i i i ij j j i i i ii b x a b a x x i n n a =+⎧=⎪⎪⎨-⎪==--⎪⎩∑L (5.3.13)元素()k kka 称为约化的主元素.将方程组(5.3.7)化为方程组(5.3.12)的过程称为消元过程.方程组(5.3.12)的求解过程(5.3.13)称为回代过程.由消元过程和回代过程求解线性方程组的方法称为Gauss 消去法.定理5.1(Gauss 消去法)设AX b =。
§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列。
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 = B
Y[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 = Y
X[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];
}
//打印X
cout << "线性方程组的解(X1,X2,X3......Xn)为:"<<endl; for( i = 0; i < N ;i++)
{
cout << X[i] <<" ";
}
}
四、算法实现
例1. 用选主元(列)的三角分解法解
解:利用软件的解答方法为:
(1)先将程序中的维数N改为3。
(2)输入1 2 3 14回车,输入2 5 2 18 回车,输入3 1 5 20回车。
(3)显示结果 1 2 3
例2. 用选主元(列)的三角分解法解
解:利用软件的解答方法为:
(1)先将程序中的维数N改为3。
(2)输入0 3 4 1回车,输入1 -1 1 2回车,输入2 1 2 3回车。
(3)显示结果 1.27273 -1.18182 0.818182。