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

（2）设对称正定阵系数阵线方程组 2、 数学原理1、 平方根法解n阶线性方程组Ax=b的choleskly方法也叫做平方根法，这里对系数矩阵A是有要求的，需要A是对称正定矩阵，根据数值分析的相关理论，如果A对称正定，那么系数矩阵就可以被分解为的形式，其中L是下三角矩阵，将其代入Ax=b中，可得：进行如下分解：那么就可先计算y,再计算x，由于L是下三角矩阵，是上三角矩阵，这样的计算比直接使用A计算简便，同时你应该也发现了工作量就转移到了矩阵的分解上面，那么对于对称正定矩阵A进行Cholesky分解，我再描述一下过程吧：如果你对原理很清楚那么这一段可以直接跳过的。设，即其中第1步，由矩阵乘法，故求得一般的，设矩阵L的前k-1列元素已经求出第k步，由矩阵乘法得于是2、 改进平方根法在平方根的基础上，为了避免开方运算,所以用计算；其中，；得 按行计算的元素及对元素公式 对于. 计算出的第行元素后，存放在的第行相置，然后再计算的第行元素，存放在的第行.的对角元素存放在的相应位置. 对称正定矩阵按分解和按分解计算量差不多，但分解不需要开放计算。求解, 的计算公式分别如下公式。 3、 程序设计1、平方根法function x=pfpf(A,b)%楚列斯基分解 求解正定矩阵的线性代数方程 A=LL 先求LY=b 再用LX=Y 即可以求出解Xn,n=size(A);L(1,1)=sqrt(A(1,1);for k=2:n L(k,1)=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+1:n L(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的初始值for k=1:n j=1:k-1; y(k)=(b(k)-L(k,j)*y(j)/L(k,k);end%解上三角方程组 LX=Y 递推公式如下，可求出X矩阵x=zeros(n,1);U=L;%求上对角矩阵for k=n:-1:1 j=k+1:n; x(k)=(y(k)-U(k,j)*x(j)/U(k,k);end A=4,2,-4,0,2,4,0,0 2,2,-1,-2,1,3,2,0 -4,-1,14,1,-8,-3,5,6 0,-2,1,6,-1,-4,-3,3 2,1,-8,-1,22,4,-10,-3 4,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.1481 -140.1127 29.7515 -60.1528 10.9120 -26.7963 5.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;for i=1:n %L的主对角元素均为1 L(i,i)=1;endfor i=1:n for j=1:n %验证A是否为对称正定矩阵if (eig(A) A=4,2,-4,0,2,4,0,0 2,2,-1,-2,1,3,2,0 -4,-1,14,1,-8,-3,5,6 0,-2,1,6,-1,-4,-3,3 2,1,-8,-1,22,4,-10,-3 4,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.1481 -140.1127 29.7515 -60.1528 10.9120 -26.7963 5.4259 -2.01854、 结果分析和讨论 平方根法和改进平方根法求解线性方程组的解为x=（121.1481，-140.1127，29.7515，-60.1528，10.9120，-26.7963，5.4259，-2.0185）T。与精确解相比较也存在很大的误差，虽然系数矩阵的对角元素都大于零，原则上可以不必选择主元，但由于矩阵的数值问题较大，不选主元的结果就是产生很大的误差，所以在求解的过程中还是应该选择主元以此消除误差，提高精度。 5、 完成题目的体会与收获 对称正定矩阵的平方根法及改进平方根法是目前解决