MATLAB-平方根法和改进平方根法求解线性方程组例题与程序

（2）设对称正定阵系数阵线方程组12345678424024000221213206411418356200216143323218122410394334411142202531011421500633421945xxxxxxxx(1,1,0,2,1,1,0,2)Tx二、数学原理1、平方根法解n阶线性方程组Ax=b的choleskly方法也叫做平方根法，这里对系数矩阵A是有要求的，需要A是对称正定矩阵，根据数值分析的相关理论，如果A对称正定，那么系数矩阵就可以被分解为的TA=LL形式，其中L是下三角矩阵，将其代入Ax=b中，可得：TLLx=b进行如下分解：TLxLbyy那么就可先计算y,再计算x，由于L是下三角矩阵，是TL上三角矩阵，这样的计算比直接使用A计算简便，同时你应该也发现了工作量就转移到了矩阵的分解上面，那么对于对称正定矩阵A进行Cholesky分解，我再描述一下过程吧：如果你对原理很清楚那么这一段可以直接跳过的。设TA=LL，即1112111112112122221222221212....................................nnnnnnnnnnnnnnaaallllaaallllaaallll其中,,1,2,...,ijjiaaijn第1步，由矩阵乘法，211111111,iialall故求得11111111,,2,3,...iialalina一般的，设矩阵L的前k-1列元素已经求出第k步，由矩阵乘法得112211kkkkkmkkikimkmikkkmmallallll,　于是12111(2,3,...,n)1(),1,2,...kkkkkkmmkikikimkmmkklalklallikknl　　　2、改进平方根法在平方根的基础上，为了避免开方运算,所以用TLDLA计算；其中，1111122.........nnndddDDDddd；得1121212212111111nnnnndllldlAlld按行计算的L元素及D对元素公式对于ni,,2,111(1,21)jijijikjkktatlji…，./(1,2,)ijijjltdj…，i-1.11iiiiikikkdatl计算出LDT的第i行元素(1,2,i-1)ijtj…，后，存放在A的第i行相置，然后再计算L的第i行元素，存放在A的第i行.D的对角元素存放在A的相应位置.对称正定矩阵A按TLDL分解和按TLL分解计算量差不多，但TLDL分解不需要开放计算。求解bLy,yxDLT的计算公式分别如下公式。1111,iiiikhkybybly2,....,in1/,/nnnniiikikkixydxydlx1,....,2,1in三、程序设计1、平方根法function[x]=pfpf(A,b)%楚列斯基分解求解正定矩阵的线性代数方程A=LL’先求LY=b再用L’X=Y即可以求出解X[n,n]=size(A);L(1,1)=sqrt(A(1,1));fork=2:nL(k,1)=A(k,1)/L(1,1);endfork=2:n-1L(k,k)=sqrt(A(k,k)-sum(L(k,1:k-1).^2));fori=k+1:nL(i,k)=(A(i,k)-sum(L(i,1:k-1).*L(k,1:k-1)))/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的初始值fork=1:nj=1: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';%求上对角矩阵fork=n:-1:1j=k+1:n;x(k)=(y(k)-U(k,j)*x(j))/U(k,k);endA=[4,2,-4,0,2,4,0,02,2,-1,-2,1,3,2,0-4,-1,14,1,-8,-3,5,60,-2,1,6,-1,-4,-3,32,1,-8,-1,22,4,-10,-34,3,-3,-4,4,11,1,-40,2,5,-3,-10,1,14,20,0,6,3,-3,-4,2,19];b=[0;-6;20;23;9;-22;-15;45];x=pfpf(A,b)x=121.1481-140.112729.7515-60.152810.9120-26.79635.4259-2.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;fori=1:n%L的主对角元素均为1L(i,i)=1;endfori=1:nforj=1: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=LDLTfori=2:nforj=1:i-1S(i,j)=A(i,j)-sum(S(i,1:j-1)*L(j,1:j-1)');L(i,1:i-1)=S(i,1:i-1)/D(1:i-1,1:i-1);endD(i,i)=A(i,i)-sum(S(i,1:i-1)*L(i,1:i-1)');endy=zeros(n,1);%x，y为n*1阶矩阵x=zeros(n,1);fori=1:ny(i)=(b(i)-sum(L(i,1:i-1)*D(1:i-1,1:i-1)*y(1:i-1)))/D(i,i);%通过LDy=b解得y的值endfori=n:-1:1x(i)=y(i)-sum(L(i+1:n,i)'*x(i+1:n));%通过LTx=y解得x的值endA=[4,2,-4,0,2,4,0,02,2,-1,-2,1,3,2,0-4,-1,14,1,-8,-3,5,60,-2,1,6,-1,-4,-3,32,1,-8,-1,22,4,-10,-34,3,-3,-4,4,11,1,-40,2,5,-3,-10,1,14,20,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.1481-140.112729.7515-60.152810.9120-26.79635.4259-2.0185四、结果分析和讨论平方根法和改进平方根法求解线性方程组的解为x=（121.1481，-140.1127，29.7515，-60.1528，10.9120，-26.7963，5.4259，-2.0185）T。与精确解相比较也存在很大的误差，虽然系数矩阵的对角元素都大于零，原则上可以不必选择主元，但由于矩阵的数值问题较大，不选主元的结果就是产生很大的误差，所以在求解的过程中还是应该选择主元以此消除误差，提高精度。五、完成题目的体会与收获对称正定矩阵的平方根法及改进平方根法是目前解决这类问题的最有效的方法之一，合理利用的话，能够产生很好的求解效果。改进平方根法较平方根法，因为不用进行开方运算，所以具有一定的求解优势。通过求解此题，学会了平方根法和改进平方根法matlab编程，使我受益匪浅。