MATLAB平方根法和改進(jìn)平方根法求解線性方程組例題與程序

1、（2）設(shè)對(duì)稱正定陣系數(shù)陣線方程組 2、 數(shù)學(xué)原理1、 平方根法解n階線性方程組Ax=b的choleskly方法也叫做平方根法，這里對(duì)系數(shù)矩陣A是有要求的，需要A是對(duì)稱正定矩陣，根據(jù)數(shù)值分析的相關(guān)理論，如果A對(duì)稱正定，那么系數(shù)矩陣就可以被分解為的形式，其中L是下三角矩陣，將其代入Ax=b中，可得：進(jìn)行如下分解：那么就可先計(jì)算y,再計(jì)算x，由于L是下三角矩陣，是上三角矩陣，這樣的計(jì)算比直接使用A計(jì)算簡(jiǎn)便，同時(shí)你應(yīng)該也發(fā)現(xiàn)了工作量就轉(zhuǎn)移到了矩陣的分解上面，那么對(duì)于對(duì)稱正定矩陣A進(jìn)行Cholesky分解，我再描述一下過程吧：如果你對(duì)原理很清楚那么這一段可以直接跳過的。設(shè)，即其中第1步，由矩陣乘法，故求

2、得一般的，設(shè)矩陣L的前k-1列元素已經(jīng)求出第k步，由矩陣乘法得于是2、 改進(jìn)平方根法在平方根的基礎(chǔ)上，為了避免開方運(yùn)算,所以用計(jì)算；其中，；得 按行計(jì)算的元素及對(duì)元素公式 對(duì)于. 計(jì)算出的第行元素后，存放在的第行相置，然后再計(jì)算的第行元素，存放在的第行.的對(duì)角元素存放在的相應(yīng)位置. 對(duì)稱正定矩陣按分解和按分解計(jì)算量差不多，但分解不需要開放計(jì)算。求解, 的計(jì)算公式分別如下公式。 3、 程序設(shè)計(jì)1、平方根法function x=pfpf(A,b)%楚列斯基分解 求解正定矩陣的線性代數(shù)方程 A=LL 先求LY=b 再用LX=Y 即可以求出解Xn,n=size(A);L(1,1)=sqrt(A(1,1

3、);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īng)的遞推公式如下，求出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,

4、k);end%解上三角方程組 LX=Y 遞推公式如下，可求出X矩陣x=zeros(n,1);U=L'%求上對(duì)角矩陣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

5、;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、改進(jìn)平方根法function x=improvecholesky(A,b,n) %用改進(jìn)平方根法求解Ax=bL=zeros(n,n); %L為n*n矩陣D=diag(n,0); %D為n*n的主對(duì)角矩陣S=L*D;for i=1:n %L的主對(duì)角元素均為1 L(i,i)=1;endfor i=1:n for j=1:n %驗(yàn)證A是否為對(duì)稱正定矩陣if (eig(A)<

6、;=0)|(A(i,j)=A(j,i) %A的特征值小于0或A非對(duì)稱時(shí)，輸出wrong disp('wrong');break;endendendD(1,1)=A(1,1); %將A分解使得A=LDLTfor i=2:n for j=1:i-1 S(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); end D(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=zer

7、os(n,1);for i=1:n y(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的值endfor i=n:-1:1 x(i)=y(i)-sum(L(i+1:n,i)'*x(i+1:n); %通過LTx=y解得x的值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,

8、-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、 結(jié)果分析和討論 平方根法和改進(jìn)平方根法求解線性方程組的解為x=（121.1481，-140.1127，29.7515，-60.1528，10.9120，-26.7963，5.4259，-2.0185）T。與精確解相比較也存在很大的誤差，雖然系數(shù)矩陣的對(duì)角元素都大于零，原則上可以不必選擇主元，但由于矩陣的數(shù)值問題較大，不選主元的結(jié)果就是產(chǎn)生很大的誤差，所以在求解的過程中還是應(yīng)該選擇主元以此消除誤差，提高精度。 5、 完成題目的體會(huì)與