平方根法求解线性方程组
数值分析15-平方根法
误差分析
定 义 :如 果 矩 阵 A 或 常 数 项 b 的 微 小 变 化 , 引 起 方 程 组 A x b解 的 巨 大 变 化 , 则 称 此 方 程 组 为 “ 病 态 ” 方 程 组 , 矩 阵 A称 为 “ 病 态 ” 矩 阵 (相对于方程组而言)。 否 则 称 方 程 组 为 “ 良 态 ” 方 程 组 , A称 为 “良态”矩阵。 矩阵的“病态”性质是矩阵本身的特性。 为了定量刻划方程组的“病态”程度,下面 对 方 程 组 A x b就 系 数 矩 阵 或 右 端 项 分 别 有 扰 动 的两种情形进行讨论。
LU 分解唯一
A = LLT
U = LT
用平方根法解线性代数方程组的算法
(1)对矩阵A进行三角分解,即A=LLT,由矩阵乘法:
对于 i = 1, 2,…, n 计算
2 lii aii lik k 1
i 1
1 2
1 i l a l l l ij ij ik jk jj 1 k
1 2 v
A v ( v 1, 2 或
A
A
2
m ax ( A T A ) . T m in ( A A ) 1 , n
当 A为 对 称 矩 阵 时 , c o n d ( A ) 2
其 中 1, n 为 A的 绝 对 值 最 大 和 绝 对 值 最 小 的 特 征 值 。
x A 1 b A 1 b
x x
b
b
此 式 表 明 ,当 右 端 项 有 扰 动 时 ,解 的 相 对 误 差 不 超 过 右 端 项 的 相 对 误 差 的 A A 1 倍 。
误差分析
平方根法
实验名称 平方根法 小组成员 一、实验目的与内容 计算结果的分析 一、实验目的与内容 1.了解平方根法的原理和意义; 2.编程实现用平方根法求解线性方程组。 二、相关背景知识介绍 平方根法又叫 Cholesky 分解法,是求解对称正定线性方程组最常用的方法之一。 我们知道,对于一般方阵,为了消除 LU 分解的局限性和误差的过分积累,而采用了选主元的 方法。但对于对称正定矩阵而言,选主元却是完全不必要的。 若线性方程组 Ax=b 的系数矩阵是对称正定的,我们可按如下的步骤求其解: 1.求 A 的 Cholesky 分解: A LLT ;[1] 2.求解 Ly=b 得到 y, 3.将 y 回代求解 LT x y 得到 x。 三、代码 用平方根法求解下列方程组. n=5,10,100,…( 到你们小组计算能力的极限)求解,对计算解 和准确解比较,观察准确程度
10
15
20
五、计算结果分析 从数据结果可以看出,利用平方根法较为准确,且计算高阶矩阵时较快。在编程过程中 有的组员对计算出 L 阵后的计算不太清楚,造成了计算结果的错误,经过讨论对平方根法的 理解更加深刻了一层。
4教 师 评 语指导教 Nhomakorabea: 年 月 日
5
2.0000,-2.3094,-1.6330,-1.2649,-1.0328 -0.8729,-0.7559,-0.6667,-0.5963,-0.5394 -0.4924,-0.4529,-0.4193,-0.3904,3.7428 2.0000,-2.3094,-1.6330,-1.2649,-1.0328 -0.8729,-0.7559,-0.6667,-0.5963,-0.5394 -0.4924,-0.4529,-0.4193,-0.3904,-0.3651 -0.3430,-0.3234,-0.3059,-0.2902,3.8644
改进的平方根法解线性方程组
# include <stdio.h>
# define N 3 /*给定方程组的阶数*/
void main ( )
{
void input_2 (float a[N][N]); /*函数声明*/
void input_1 (float a[N]);
void output_1 (float a[N]);
float a[N][N]; /*定义A,b */
float b[N];
printf ("矩阵A please input data:\n");
input_2 (a); /*调用二维数组输入函数*/
printf ("向量b please input data:\n");
input_1 (b); /*调用一维数组输入函数*/
{
t=l[i][k]*d[k]*l[j][k];
s=s+t;
}
l[i][j]=(a[i][j]-s)/d[j];
s=0;
}
for (k=0;k<i;k++)
{
t=l[i][k]*l[i][k]*d[k];
s=s+t;
}
d[i]=a[i][i]-s;
s=0;
}
float y[N]; /*求解Ly=b */
{
for (k=i+1;k<N;k++)
{
t=l[k][i]*x[k];
s=s+t;
}
x[i]=y[i]/d[i]-s;
s=0;
}
printf ("方程组的解为:\n");
平方根法计求解线性方程组
解线性n 阶方程组直接法—Cholesky 方法解n 阶线性方程组Ax=b 的choleskly 方法也叫做平方根法,这里对系数矩阵A 是有要求的,需要A 是对称正定矩阵,根据数值分析的相关理论,如果A 对称正定,那么系数矩阵就可以被分解为的T A=L L ∙形式,其中L 是下三角矩阵,将其代入Ax=b 中,可得:T LL x=b进行如下分解:T L x L b y y ⎧=⎨=⎩那么就可先计算y,再计算x ,由于L 是下三角矩阵,是T L 上三角矩阵,这样的计算比直接使用A 计算简便,同时你应该也发现了工作量就转移到了矩阵的分解上面,那么对于对称正定矩阵A 进行Cholesky 分解,我再描述一下过程吧: 如果你对原理很清楚那么这一段可以直接跳过的。
设T A=L L ∙,即1112111112112122221222221212....................................n n n n n n nn n n nn nn a a a l l l l a a a l l l l a a a l l l l ⎡⎤⎡⎤⎡⎤⎢⎥⎢⎥⎢⎥⎢⎥⎢⎥⎢⎥=⎢⎥⎢⎥⎢⎥⎢⎥⎢⎥⎢⎥⎣⎦⎣⎦⎣⎦其中,,1,2,...,ij ji a a i j n ==第1步,由矩阵乘法,211111111,i i a l a l l == 故求得111111,2,3,...i i a l l i n a === 一般的,设矩阵L 的前k-1列元素已经求出第k 步,由矩阵乘法得112211k k kk km kkik im km ik kk m m a l l a l l l l --===+=+∑∑, 于是11(2,3,...,n)1(),1,2,...kk k ik ik im km m kk l k l a l l i k k n l -=⎧=⎪⎪=⎨⎪=-=++⎪⎩∑ 注意到21kkk km m a l ==∑,于是有21max km kk ii i nl a a ≤≤≤≤ 这充分说明分解过程中元素2km l 的平方不会超过系数矩阵A 的最大对角元,因而分解过程中舍入误差的放大收到了控制,用平方根法解对称正定方程组时可以不考虑选主元问题。
平方根法解线性方程组
a (k 1) ij
a / a (k1) ik
(k 1) kk
*
a(k kj
1)
,
i,
j
k
1,
,n
bi(k)
b(k 1) i
a / a (k1) ik
(k 1) kk
*bk(k1) , i k 1,
, n.
(3.9)
由(3.8)逐个回代,可得(3.1)的解
b(2) i
b (1) i
b2(1)mi2 ,
i,
j
3, 4,
, n.
重复上述过程,可以得到与(3.1)等价的上三角方程组:
a11x1
a12 x2 a12 x3
a (1) 22
x2
a (1) 22
x3
a (2) 33
x3
a1n xn b1
k 2
k 1
3
(n 1) (n 2)
2
1
1 2
n(n
1); (for
right
item)
除:(n 1) (n 2)
• 回代过程:
2 1
1 2
n(n
1)
1 2
(n
1)(multiply)
n(divide)
1 2
n(n
1).
• 运算量:1 n(n2 1) 1 n(n 1) 1 n(n 1) n3 3 n2 5 n
2. 即使高斯消元法可行,如果 ak(kk很1)小,运算中用它作分 母会导致其它元素数量计的严重增长和舍入误差的扩散。
平方根法和改进平方根法求解线性方程组例题与程序
平方根法和改进平方根法求解线性方程组例题与程序2、数学原理1、平方根法解n阶线性方程组Ax=b的choleskly方法也叫做平方根法,这里对系数矩阵A是有要求的,需要A是对称正定矩阵,根据数值分析的相关理论,如果A对称正定,那么系数矩阵就可以被分解为的形式,其中L是下三角矩阵,将其代入Ax二b 中,可得:进行如下分解:那么就可先计算y,再计算x,由于L 是下三角矩阵,是上三角矩阵,这样的计算比直接使用A计算简便,同时你应该也发现了工作量就转移到了矩阵的分解上面,那么对于对称正定矩阵A进行Cholesky分解,我再描述一下过程吧:如果你对原理很清楚那么这一段可以直接跳过的。
设,即其中第1步,由矩阵乘法,故求得一般的,设矩阵L的前k-l列元素已经求出第k步,由矩阵乘法得于是2、改进平方根法在平方根的基础上,为了避免开方运算,所以用计算;其中,;得按行计算的元素及对元素公式对于、、计算出的第行元素后,存放在的第行相置,然后再计算的第行元素,存放在的第行、的对角元素存放在的相应位置、对称正定矩阵按分解和按分解计算量差不多,但分解不需要开放计算。
求解, 的计算公式分别如下公式。
3、程序设计1、平方根法function[x]二pfpf(A,b)%楚列斯基分解求解正定矩阵的线性代数方程A二LL'先求LY二b再用L' X二Y即可以求出解X[n,n]=size(A) ;L(1, l)=sqrt(A(l, 1)) ;for k=2:nL(k, l)=A(k, 1)/L(1,1) ; endfor k=2: n-1 L(k, k) =sqrt (A(k, k)_sum(L(k, 1:k-1) > 2)) ; for i=k+l:n L(i,k) = (A(i,k)-stun(L(i, l:kT)、*L(k, l:kT)))/L(k,k); endendL (n, n)=sqrt (A(n, n) -sum(L(n, 1: n-1)、2)) ; %解下三角方程组Ly二b相应的递推公式如下,求出y矩阵y二zeros (n, 1) ;%先生成方程组的因变量的位置,给定y的初始值for k=l: n j=l: k-1; y (k) = (b(k)-L(k, j) *y (j))/L(k, k) ; end%解上三角方程组L' X=Y 递推公式如下,可求出X 矩阵x二zeros (n, 1) ;U=L;%求上对角矩阵for k=n: -1:1 j=k+l:n; x(k) = (y (k)-U(k,j)*x(j))/U(k,k);end » A二[4,2,-4,0,2,4,0,02,2,-1,- 2,1,3,2,01,14,1,-8,-3,5,6 0,-2,1,6,-1,-4,-3,32,1,-8,- 1,22,4,-10,-34,3,-3,-4,4,11,1,-4 0,2,5,-3,-10,1,14,2 0,0,6,3,-3,-4,2,19];>〉b=[0;-6;20;23;9;-22;-15;45];» x=pfpf(A,b)x =121、148160、152810、91202、01852、改进平方根法function[x]=improvecholesky(A,b,n)%用改进平方根法求解Ax=bL=zeros(n,n); %L为n*n矩阵D=diag(n, 0) ; %D 为n*n 的主对角矩阵S=L*D; for i=l: n %L 的主对角元素均为 1 L(i, i) = l; endfor i=l: n for j=l:n %验证 A 是否为对称正定矩阵if (eig(A)<=0) (A(i, j)~=A(j,i))%A的特征值小于0或A非对称时,输出wrongdisp(wrong) ; break; endendendD(1,1)=A(1, 1) ;%将A 分解使得A=LDLTfor i二2:n for j=l:i~l S(i,j)=A(i,j)-sum(S(i,l:j~ l)*L(j,l: j-1)); L(i,l:i-l)=S(i,l:i-l)/D(l:i-l,l:i-l); end D(i, i)=A(i, i)-sum(S(i, 1: i-l)*L(i, 1: i-1)) ; endy=zeros (n, 1) ; % x, y 为阶矩阵x=zeros (n, 1); for i=l: n y(i) = (b(i)-sum(L(i, 1: i-l)*D(l: i-1,1: i-l)*y (1: i-1)))/D(i, i) ;%通过LDy=b 解得y 的值endfor i=n:-1:1 x(i)=y(i)-sum(L(i+l:n, i)*x(i+l:n)) ;%通过LTx=y 解得x 的值end» A二[4,2,-4,0,2,4,0,02,2,-1,-2,1,3,2,01,14,1,-8,- 3,5,6 0,-2,1,6,-1,-4,-3,32,1,-8,-1,22,4,-10,-34,3,-3,- 4,4,11,1,-4 0,2,5,-3,-10,1,14,2 0,0,6,3,-3,-4,2,19];>> b=[0;-6;20;23;9;-22;-15;45];» n=8;»x=improvecholesky(A,b,n)x =121、148160、152810、91202、01854、结果分析和讨论平方根法和改进平方根法求解线性方程组的解为x二(121、1481, -140、1127,29、7515,-60、1528,10、9120,-26、7963,5、4259, -2、0185) To与精确解相比较也存在很大的误差,虽然系数矩阵的对角元素都大于零,原则上可以不必选择主元,但由于矩阵的数值问题较大,不选主元的结果就是产生很大的误差,所以在求解的过程中还是应该选择主元以此消除误差,提高精度。
第2章解线性方程组的直接方法5_6
25 29
9 10 9
y1
b1 l11
9 6
y2
b2
l21 l22
y1
2
b3 l3k yk
y3
k 1
l33
10 29
y1
b1 l11
i1
bi lik yk
yi
k 1
lii
10 7 * 9
6
3
29
174
6
13
即
y ( y1 , y2 , y3 )T
( 9 , 6
u~n1,n u~n1,n1
1
3
D diag(u~11 ,u~22 , ,u~nn ) [diag( u~11 , u~22 , , u~nn )]2
U~
DU
1
D
1 2
D
1
2U 1
11
A L~U~ L%(D 2 D 2U1)
11
而 A 为对称正定阵 , AT A U1T (D 2 D 2 L%T )
l11 li1 ln1
LT
lii
llnnni
------(15) i n 1, ,2,1
对称正定方程 组的平方根法
9
1. 输入系数矩阵A, 右端项b; 2. 作A的Cholesky分解:
for j=1,2,…n 2.1 如果j>1 则 for i= j,j+1,j+2,…,
j 1
§ 2.5 平方根法
一、对称正定矩阵的三角分解(Cholesky分解)
若n阶矩阵A为对称正定矩阵 则det( A) 0, AT A
且A的顺序主子式 det Ak 0, k 1,2, , n
因此A可以进行 LU分解(或Doolittle分解)
改进平方根法1
进平方根法解Ax=b一、题目:改进平方根法解线性方程组Ax=b二、算法:1.输入n阶系数矩阵A,元素数为n的向量b2.判断系数矩阵是否是对称正定矩阵:①a ij≠a ji,或特征值小于0跳出,则不是对称正定矩阵,程序终止,输出“该方程组不能用改进平方根法求解”;②如果是对称正定矩阵,向下进行第3步3.①d11=a11②对于i=2,3,4…nS ij=a ij-sum(S ik*l kj) (k=1 to j-1)l ij=S ij/d jjd ii=a ii-sum(S ik*l ik) (k=1 to i-1)③解下三角形方程组LDY=by i=(b i-sum(d kk*l ik*y k))/d ii(k=1 to i-1)④解上三角形方程组x i=(y i-sum(l ik*x k)) (k=i+1 to n ; i=n,n-1…1)4.输出解向量x1,x2…x n三、源程序(见附件gaijin.m)四、数值例子①输入A:[1 2 1 -3;2 5 0 -5;1 0 14 1;-3 -5 1 15]b:[1 2 16 8]输出:x=(1,1,1,1)T②输入A:[6 2 1 -1;2 4 1 0;1 1 4 -1;-1 0 -1 3] b:[6 -1 5 -5]输出:x=(1.0000,-1.0000,1.0000,-1.0000)T③输入A:[1 2 3;4 5 6;7 8 9]b:[4 4 4]输出:该方程组不可以用改进平方根法求解五、结果及其分析第一个例子来源于讲义,第二个例子来源于习题书,可以看出运行结果是准确的,但第一个例子的结果是精确的1,第二个是1.0000,我认为这是在运算过程中的存储有舍入误差导致最终结果是近似得来的1.000而不是精确的1。
第三个例子是自己构造的非对称正定矩阵,输出结果也是令人满意地。
在经历过列主元消去法程序编写之后,改进平方根法的编写就容易了许多,之前调试程序是出现的错误基本分为以下几类:1.没有给y或者D定义一个存储空间导致运行时寻找不到这些中间存储向量2.逗号或分号使用不当由于事先写出了算法,在循环中倒出问题比较少了。
(完整word版)线性方程组的平方根解法
浅析线性方程组的平方根解法在求解线性方程组时,直接解法有顺序高斯消元法、列主元高斯消元法、全主元高斯消元法、高斯约当消元法、消元形式的追赶法、LU 分解法、矩阵形式的追赶法,当我们遇到对称正定线性方程组时,我们就要用到平方根法(对称LLT 分解法)来求解,为了熟悉和熟练运用平方根法求解线性方程组,下面对运用平方根法求解线性方程组进行解析。
一、运用平方根法求解线性方程组涉及到的定理及定义我们在运用平方根法求解线性方程组时,要判定线性方程组Ax=b 的系数矩阵A 是否是对称正定矩阵,那么我们就要了解正定矩阵的性质和如下定理及定义:1、由线性代数知,正定矩阵具有如下性质:1) 正定矩阵A 是非奇异的2) 正定矩阵A 的任一主子矩阵也必为正定矩阵 3) 正定矩阵A 的主对角元素均为正数 4) 正定矩阵 A 的特征值均大于零 5) 正定矩阵A 的行列式必为正数定义一 线性方程组Ax=b 的系数矩阵A 是对称正定矩阵,那么Ax=b 是对称正定线性方程组。
定义二 如果方阵A 满足A=AT ,那么A 是对称阵。
2.1.4 平方根法和改进的平方根法如果A 是n 阶对称矩阵,由定理2还可得如下分解定理:定理2 若A 为n 阶对称矩阵,且A 的各阶顺序主子式都不为零,则A 可惟一分解为:A =LDLT ,其中L 为单位下三角阵,D 为对角阵。
证明 因为A 的各阶顺序主子式都不为零,所以A 可惟一分解为:A =LU 因为 ,所以可将 U 分解为:⎪⎪⎪⎪⎪⎭⎫ ⎝⎛=nn u u u U 2211⎪⎪⎪⎪⎪⎪⎪⎭⎫ ⎝⎛11122211112 u u u u u u n nn n 1DU =其中 D 为对角矩阵,U1为单位上三角阵.于是:A =LDU1=L(DU1)因为A 为对称矩阵,所以,A =AT =U1TDTLT =U1T(DLT),由 A 的 LU 分解的惟一性即得:L =U1T ,即U1=LT ,故A =LDLT 。
平 方 根 法
L是单位下三角阵, U是上三角阵, 将U再分解
u11
u 22
1
u
nn
u12
u11 1
u1n
u11
u n1,n
DU 0
u n1,n1
1
其中D为对角阵, U0为单位上三角阵,于是
A = L U = L D U0
又
A = AT = U0TD LT
数值计算方法
点l11是 需由a11要此进1例, 行可开l以21方看a运l1出211 算,11。平 1为,方避根免l法31 开解al方13正11 运定12算方 ,2程我组们的改缺
A LDL 用 解ll3232单成位a3三3a2l角2321 阵ll32222作1为11分24解T1阵4,13即的l把32 形对a式称32 l,正22l3其1定l21中矩 0阵11A分2 2
l11 l21 ln1
l22
l
n
2
l
nn
按矩阵乘法展开,可逐行求出分解矩阵L的元素,计
算公式是对于i=1,2,…,n
i1
1
lii (aii li2k ) 2
k 1
i 1
a ji l jk lik
l ji
k 1
lii
j=i+1, i+2,…,n
这一方法称为平方根法,又称乔累斯基(Cholesky)分
数值计算方法
平方根法 工程实际计算中,线性方程组的系数矩阵
常常具有对称正定性,其各阶顺序主子式及 全部特征值均大于0。矩阵的这一特性使它的 三角分解也有更简单的形式,从而导出一些 特殊的解法,如平方根法与改进的平方根法。
求解线性方程组的直接方法算法实验报告
求解线性方程组的直接方法(2学时)一 实验目的1.掌握求解线性方程组的高斯消元法及列主元素法;2. 掌握求解线性方程组的克劳特法;3. 掌握求解线性方程组的平方根法。
二 实验内容1.用高斯消元法求解方程组(精度要求为610-=ε):1231231233272212240x x x x x x x x x -+=⎧⎪-+-=-⎨⎪-+=⎩2.用克劳特法求解上述方程组(精度要求为610-=ε)。
3. 用平方根法求解上述方程组(精度要求为610-=ε)。
4. 用列主元素法求解方程组(精度要求为610-=ε):1231231233432222325x x x x x x x x x -+=⎧⎪-+-=⎨⎪--=-⎩三 实验步骤(算法)与结果1. #include<stdio.h>main(){ float a11,a12,a13,b1,a21,a22,a23,b2,a31,a32,a33,b3,l21,l31,l32,u11,u12,u13,u22,u23,u33,z1,z2,z3,x1,x2,x3;printf("enter a11,a12,a13,b1,a21,a22,a23,b2,a31,a32,a33,b3 :"); scanf("%f%f%f%f%f%f%f%f%f%f%f%f/n",&a11,&a12,&a13,&b1,&a21,&a22, &a23,&b2,&a31,&a32,&a33,&b3 );l21=a12/a11;l31=a31/a11;u22=a22-l21*a12;l32=(a32-l31*a12)/u22;u23=a23-l21*a13;u33=a33-l31*a13-l32*u23;z2=b2-l21*b1;z3=b3-l31*b1-l32*z2;printf("u22=%fu23=%fz2=%fu33=%fz3=%f`````",u22,u23,z2,u33,z3); x3=z3/u33;x2=(z2-u23*x3)/u22;x1=(b1-a13*x3-a12*x2)/a11;printf("x1=%f x2=%f x3=%f ",x1,x2,x3);return 0;}2.#include<stdio.h>main(){float a11,a12,a13,b1,a21,a22,a23,b2,a31,a32,a33,b3,l22,l32,l33,u11,u12,u13,u22,u23,u33,z1,z2,z3,x1,x2,x3;printf("enter a11,a12,a13,b1,a21,a22,a23,b2,a31,a32,a33,b3 :");scanf("%f%f%f%f%f%f%f%f%f%f%f%f/n",&a11,&a12,&a13,&b1,&a21,&a22,&a23,&b2,&a31,&a32,&a33,&b3 );u11=1;u22=1; u33=1;u12=a12/a11;u13=a13/a11;z1=b1/a11;l22=a22-a21*u12;u23=(a23-a21*u13)/l22;z2=(b2-a21*z1)/l22;l32=a32-a31*u12;l33=a33-a31*u13-l32*u23;z3=(b3-a31*z1-l32*z2)/l33;printf("u11=%f u12=%f u13=%f z1=%f u22=%f u23=%f z2=%f u33=%f z3=%f------",u11,u12,u13,z1,u22,u23,z2,u33,z3);x3=z3;x2=z2-u23*x3;x1=z1-u13*x3-u12*x2;printf("x1=%f x2=%f x3=%f ",x1,x2,x3);getch();return 0;}3. #include<stdio.h>#include<math.h>{float a11,a12,a13,b1,a21,a22,a23,b2,a31,a32,a33,b3,l11,l12,l13,l23,l21,l22,l31,l32,l33,z1,z2,z3,x1,x2,x3;printf("enter a11,a12,a13,b1,a21,a22,a23,b2,a31,a32,a33,b3 :");scanf("%f%f%f%f%f%f%f%f%f%f%f%f/n",&a11,&a12,&a13,&b1,&a21,&a22, &a23,&b2,&a31,&a32,&a33,&b3 );l11=sqrt(a11);l21=a21/l11; l31=a31/l11;l22=sqrt(a22-l21*l21);l32=(a32-l21*l31)/l22;l33=sqrt(a33-l31*l31-l32*l32);z1=b1/l11;z2=(b2-l21*z1)/l22;z3=(b3-l31*z1-l32*z2)/l33;printf("l11=%f z1=%f l22=%f z2=%f l33=%fz3=%f---",l11,z1,l22,z2,l33,z3);x3=z3/l33;x2=(z2-l32*x3)/l22;x1=(z1-l31*x3-l21*x2)/l11;printf("x1=%f x2=%f x3=%f ",x1,x2,x3);getch();return 0;}4. #include "stdio.h"#include "math.h"main(){ float a11,a12,a13,b1,a21,a22,a23,b2,a31,a32,a33,b3,l21,l31,A22,A23,d1,A32,A33,d2,l32,a,d3,x1,x2,x3,A,B,C,D;printf("enter a11,a12,a13,b1,a21,a22,a23,b2,a31,a32,a33,b3:"); scanf("%f%f%f%f%f%f%f%f%f%f%f%f", &a11,&a12,&a13,&b1,&a21,&a22,&a23,&b2,&a31,&a32,&a33,&b3);if(fabs(a11)<fabs(a21)){ if(fabs(a11)>fabs(a31))A=a11;a11=a31;a31=A;B=a12;a12=a32;a32=B;C=a13;a13=a33;a33=C;D=b1;b1=b3;b3=D ;}if (fabs(a11)<fabs(a21)){if(fabs(a21)>fabs(a31)){A=a11;a11=a21;a21=A;B=a12;a12=a22;a22=B;C=a13;a13=a23;a23=C;D=b1;b1=b2;b2=D ;}elseA=a11;a11=a31;a31=A;B=a12;a12=a32;a32=B;C=a13;a13=a33;a33=C;D=b1;b1=b3;b1=D ;}printf("now a11=%f a12=%f a13=%f b1=%f\n",a11,a12,a13,b1); printf("now a21=%f a22=%f a23=%f b2=%f\n",a21,a22,a23,b2); printf("now a31=%f a32=%f a33=%f b3=%f\n",a31,a32,a33,b3);l21=a21/a11; l31=a31/a11;A22=a22-l21*a12;A23=a23-l21*a13;d1=b2-l21*b1;A32=a32-l31*a12; A33=a33-l31*a13;d2=b3-l31*b1;if(fabs(A22)>fabs(A32)){ l32=A32/A22;a=A33-l32*A23;d3=d2-l32*d1;x3=d3/a;x2=(d1-A23*x3)/A22;x1=(b1-a13*x3-a12*x2)/a11;}else l32=A22/A32;a=A23-l32*A33;d3=d1-l32*d2;x3=d3/a;x2=(d2-A33*x3)/A32;x1=(b1-a13*x3-a12*x2)/a11;printf("x1=%f x2=%f x3=%f\n",x1,x2,x3);getch(); return 0; }。
MATLAB-平方根法和改进平方根法求解线性方程组例题与程序演示教学
M A T L A B-平方根法和改进平方根法求解线性方程组例题与程序(2)设对称正定阵系数阵线方程组12345678424024000221213206411418356200216143323218122410394334411142202531011421500633421945x x x x x x x x -⎡⎤⎡⎤⎡⎤⎢⎥⎢⎥⎢⎥---⎢⎥⎢⎥⎢⎥⎢⎥⎢⎥⎢⎥----⎢⎥⎢⎥⎢⎥----⎢⎥⎢⎥⎢⎥=⎢⎥⎢⎥⎢⎥----⎢⎥⎢⎥⎢⎥----⎢⎥⎢⎥⎢⎢⎥⎢⎥⎢---⎢⎥⎢⎥⎢--⎢⎥⎢⎢⎥⎣⎦⎣⎦⎣⎦⎥⎥⎥⎥ (1,1,0,2,1,1,0,2)T x *=--二、数学原理 1、平方根法解n 阶线性方程组Ax=b 的choleskly 方法也叫做平方根法,这里对系数矩阵A 是有要求的,需要A 是对称正定矩阵,根据数值分析的相关理论,如果A 对称正定,那么系数矩阵就可以被分解为的T A=L L •形式,其中L 是下三角矩阵,将其代入Ax=b 中,可得:T LL x=b 进行如下分解:T L xL by y ⎧=⎨=⎩ 那么就可先计算y,再计算x ,由于L 是下三角矩阵,是T L 上三角矩阵,这样的计算比直接使用A 计算简便,同时你应该也发现了工作量就转移到了矩阵的分解上面,那么对于对称正定矩阵A 进行Cholesky 分解,我再描述一下过程吧: 如果你对原理很清楚那么这一段可以直接跳过的。
设T A=L L •,即1112111112112122221222221212....................................n n n n n n nn n n nn nn a a a l l l l aa a l l l l a a a l l l l ⎡⎤⎡⎤⎡⎤⎢⎥⎢⎥⎢⎥⎢⎥⎢⎥⎢⎥=⎢⎥⎢⎥⎢⎥⎢⎥⎢⎥⎢⎥⎣⎦⎣⎦⎣⎦其中,,1,2,...,ij ji a a i j n ==第1步,由矩阵乘法,211111111,i i a l a l l ==g 故求得111111,2,3,...i i a l l i n a === 一般的,设矩阵L 的前k-1列元素已经求出 第k 步,由矩阵乘法得112211k k kk kmkkik im km ik kkm m a l l a l l l l --===+=+∑∑, 于是11(2,3,...,n)1(),1,2,...kk k ik ik im km m kk l k l a l l i k k n l -=⎧=⎪⎪=⎨⎪=-=++⎪⎩∑ 2、改进平方根法在平方根的基础上,为了避免开方运算,所以用TLDL A =计算;其中,11122.........n d D D D d ⎤⎤⎡⎤⎥⎥⎢⎥⎢⎥⎢⎥⎢⎥⎢⎥⎢⎥⎢⎥===⎢⎥⎢⎥⎢⎥⎢⎥⎢⎥⎢⎥⎢⎢⎢⎥⎣⎦⎣⎣;得1121212212111111n n n n n d l l l d l A l l d ⎡⎤⎡⎤⎡⎤⎢⎥⎢⎥⎢⎥⎢⎥⎢⎥⎢⎥=⎢⎥⎢⎥⎢⎥⎢⎥⎢⎥⎢⎥⎣⎦⎣⎦⎣⎦L L M MO O O M L按行计算的L 元素及D 对元素公式 对于n i ,,2,1Λ=11(1,21)j ij ij ik jk k t a t l j i -==-=-∑…,./(1,2,)ij ij j l t d j ==…,i-1.11i i ii ik ikk d a t l -==-∑计算出LD T =的第i 行元素(1,2,i-1)ij t j =…,后,存放在A 的第i 行相置,然后再计算L 的第i 行元素,存放在A 的第i 行.D 的对角元素存放在A 的相应位置.对称正定矩阵A 按T LDL 分解和按T LL 分解计算量差不多,但T LDL 分解不需要开放计算。
用改进的平方根法解方程组例题
改进的平方根法是一种用于解方程组的数值方法,它可以帮助我们更快速地找到方程组的解。
在本文中,我们将首先介绍改进的平方根法的原理和基本概念,然后结合例题来具体说明如何应用这一方法。
通过本文的阐述,相信你对改进的平方根法会有更加深入的理解。
1. 改进的平方根法原理与基本概念改进的平方根法是一种用于解线性方程组的数值方法,它通过对系数矩阵进行分解,将原方程组转化为两个容易求解的三角方程组,再通过回代的方式得到原方程组的解。
改进的平方根法相较于传统的平方根法,能够在计算过程中减小误差,提高计算精度。
2. 应用改进的平方根法解方程组例题考虑以下线性方程组:```3x + 2y - z = 12x - 2y + 4z = -2- x + 0.5y - z = 0```我们需要对系数矩阵进行特征值分解,将其转化为一个上三角矩阵和一个下三角矩阵的乘积。
通过回代的方式求解两个三角方程组,最终得到原方程组的解。
3. 改进的平方根法的个人观点和理解改进的平方根法是一种高效而有效的数值方法,它在解线性方程组时能够降低计算误差,提高解的精度。
相较于其他数值方法,改进的平方根法在实际应用中具有一定的优势,特别是对于大型稀疏方程组的求解。
总结回顾通过本文对改进的平方根法的介绍和例题分析,相信你对这一数值方法有了更深入的理解。
在实际应用中,我们可以根据方程组的特点和求解需求选择合适的数值方法,以便更快速、准确地得到方程组的解。
通过这篇文章,我希望读者能够对改进的平方根法有一个全面、深刻和灵活的理解。
希望这篇文章对你有所帮助,如果有任何问题或疑惑,欢迎随时交流讨论。
改进的平方根法是一种用于解方程组的数值方法,它能够帮助我们更快速地找到方程组的解,并且在计算过程中减小误差,提高计算精度。
在本文中,我们将继续深入探讨改进的平方根法的原理和基本概念,并结合更多例题来具体说明如何应用这一方法。
我们还将探讨改进的平方根法在实际应用中的优势和局限性,以及一些其他相关的知识点。
解对称正定矩阵线性方程组的平方根法
i− i−1
k =1 n
例4.4.3(P91) 4.4.3(P91)
优点: 数值稳定。 优点:1、数值稳定。 3 计算量小, 次乘除法,是一般矩阵A的 分解 2、计算量小,大约为 n 6次乘除法,是一般矩阵 的LU分解 计算量的一半。 计算量的一半。 缺点:计算l 时要开平方。 缺点:计算 ii时要开平方。 6.3 改进的平方根法 T 阶对称正定矩阵A有分解 为单位下三角阵, 设n阶对称正定矩阵 有分解A = LDL。其中L为单位下三角阵, 阶对称正定矩阵 D为对角阵。即 为对角阵。 为对角阵
11
21
l 22
L l n2 O M l nn
n1
则aij = ∑ lik l jk = ∑ lik l jk + l ij l jj
L 的元素 l ij。
经 (当 k > j时, l jk = 0 ), n 步可直接求得
求解对称正定方程组Ax=b的平方根法(计算公式) : 的平方根法(计算公式) 求解对称正定方程组 n 1.分解计算 A = LLT:
j −1 k =1
1) t ij = a ij − ∑ t ik l jk , ( j = 1,2,L , i − 1)
2) l ij = t ij d j , ( j = 1,2, L , i − 1)
3) d i = a ii − ∑ t ik l ik
k =1 i −1
a11 a21 A= L a n1
a11 a21 A 令 = L a n1
并有d 1 = a11, 再由 a ij = ∑ ( LD ) ik ( LT ) kj = ∑ ( l ik d k )l jk + l ij d j (当 k > j 时, l jk = 0 )
(完整word版)线性方程组的平方根解法
浅析线性方程组的平方根解法在求解线性方程组时,直接解法有顺序高斯消元法、列主元高斯消元法、全主元高斯消元法、高斯约当消元法、消元形式的追赶法、LU 分解法、矩阵形式的追赶法,当我们遇到对称正定线性方程组时,我们就要用到平方根法(对称LLT 分解法)来求解,为了熟悉和熟练运用平方根法求解线性方程组,下面对运用平方根法求解线性方程组进行解析。
一、运用平方根法求解线性方程组涉及到的定理及定义我们在运用平方根法求解线性方程组时,要判定线性方程组Ax=b 的系数矩阵A 是否是对称正定矩阵,那么我们就要了解正定矩阵的性质和如下定理及定义:1、由线性代数知,正定矩阵具有如下性质:1) 正定矩阵A 是非奇异的2) 正定矩阵A 的任一主子矩阵也必为正定矩阵 3) 正定矩阵A 的主对角元素均为正数 4) 正定矩阵 A 的特征值均大于零 5) 正定矩阵A 的行列式必为正数定义一 线性方程组Ax=b 的系数矩阵A 是对称正定矩阵,那么Ax=b 是对称正定线性方程组。
定义二 如果方阵A 满足A=AT ,那么A 是对称阵。
2.1.4 平方根法和改进的平方根法如果A 是n 阶对称矩阵,由定理2还可得如下分解定理:定理2 若A 为n 阶对称矩阵,且A 的各阶顺序主子式都不为零,则A 可惟一分解为:A =LDLT ,其中L 为单位下三角阵,D 为对角阵。
证明 因为A 的各阶顺序主子式都不为零,所以A 可惟一分解为:A =LU 因为 ,所以可将 U 分解为:⎪⎪⎪⎪⎪⎭⎫ ⎝⎛=nn u u u U 2211⎪⎪⎪⎪⎪⎪⎪⎭⎫ ⎝⎛11122211112 u u u u u u n nn n 1DU =其中 D 为对角矩阵,U1为单位上三角阵.于是:A =LDU1=L(DU1)因为A 为对称矩阵,所以,A =AT =U1TDTLT =U1T(DLT),由 A 的 LU 分解的惟一性即得:L =U1T ,即U1=LT ,故A =LDLT 。
线性方程组的特殊解法——平方根法
, =口/l 23 一n j nl , n l )
从 而确定 了第 的算式 为 : 列
r J一 _ 2 ) 厂 - , 口 … 广 — —
l= ∑[ ) (.l+,, ≠ l (一 ・ / f,' 2 , 口 正 2 =+_ …l j 『 √
因此 , 方根 法 求解对 称正 定矩 阵线 性方 程组 AX= : 用平 b ① 利用 递 推 公 式 ()() 矩 阵 A分 解 为 A 工 , 1 、2将 ②对称正定矩阵线性方程组中图分类号o1222文献标识码a文章编号1674098x200812a023901得到l用公式3解lyb得所以原方程组的解为4结语此算法的优点是在计算过程中不需选主元当n较大时约需n次乘除法运算相当于高斯消元法计算量的一半并且数值稳定储存量小
峰 朱 论 坛
Sic d e nlyn v cnea c og Io e n T h o n
一
I
0
有 方 一 平 根 。 效 法 一 方 法
可以证明对于对称正定矩阵A, 可以唯一地分解成 =工 r 其中L £,
; = 云
兰
,
得到L =
l、 一 3
, =/ 32 : ,
2平方根法递推公式
是非奇 异 下三 角形矩 阵 。因 篇幅有 限 , 明从 略) ( 证 下 面给 出平方 根 法 的递 推算 法 。 设
…
用公式(解 】西得,= , o 3 ) , , I 2 2 . _
再用公式( 解方程 4 )
LX T得 = l , ; , ; , r= , 一/ 3 一
…
=
;
:
…
,
O …
根 据矩 阵乘法 法
线性方程组的平方根解法
浅析线性方程组的平方根解法在求解线性方程组时,直接解法有顺序高斯消元法、列主元高斯消元法、全主元高斯消元法、高斯约当消元法、消元形式的追赶法、LU 分解法、矩阵形式的追赶法,当我们遇到对称正定线性方程组时,我们就要用到平方根法(对称LLT 分解法)来求解,为了熟悉和熟练运用平方根法求解线性方程组,下面对运用平方根法求解线性方程组进行解析。
一、运用平方根法求解线性方程组涉及到的定理及定义我们在运用平方根法求解线性方程组时,要判定线性方程组Ax=b 的系数矩阵A 是否是对称正定矩阵,那么我们就要了解正定矩阵的性质和如下定理及定义:1、由线性代数知,正定矩阵具有如下性质:1) 正定矩阵A 是非奇异的2) 正定矩阵A 的任一主子矩阵也必为正定矩阵 3) 正定矩阵A 的主对角元素均为正数 4) 正定矩阵 A 的特征值均大于零 5) 正定矩阵A 的行列式必为正数定义一 线性方程组Ax=b 的系数矩阵A 是对称正定矩阵,那么Ax=b 是对称正定线性方程组。
定义二 如果方阵A 满足A=AT ,那么A 是对称阵。
2.1.4 平方根法和改进的平方根法如果A 是n 阶对称矩阵,由定理2还可得如下分解定理:定理2 若A 为n 阶对称矩阵,且A 的各阶顺序主子式都不为零,则A 可惟一分解为:A =LDLT ,其中L 为单位下三角阵,D 为对角阵。
证明 因为A 的各阶顺序主子式都不为零,所以A 可惟一分解为:A =LU 因为 ,所以可将 U 分解为:⎪⎪⎪⎪⎪⎭⎫ ⎝⎛=nn u u u U O 2211⎪⎪⎪⎪⎪⎪⎪⎭⎫⎝⎛11122211112M O ΛΛu u u u u u n nn n 1DU =其中 D 为对角矩阵,U1为单位上三角阵.于是:A =LDU1=L(DU1)因为A 为对称矩阵,所以,A =AT =U1TDTLT =U1T(DLT),由 A 的 LU 分解的惟一性即得:L =U1T ,即U1=LT ,故A =LDLT 。