cholesky分解法

Cholesky分解對稱正定矩陣三角分解的一個基本方法。該方法是LU的特殊形式，其中L和U轉置。利用這個性質我們可以方便的解得L的值。然後再利用LU法的回代過程的到方程組的解，其代碼如下：

#include<iostream.h><br />#include<math.h><br />#include<process.h><br />class cholesky<br />{<br />private:<br />int i,j,k,n;<br />double sum,*b,*d,*x,**a,eps;<br />public:<br />void cholesky_input();<br />void cholesky_decomposition();<br />void cholesky_output();<br />~cholesky()<br />{<br />delete []b;<br />delete []d;<br />delete []x;<br />for(i=0;i<n;i++)<br />{<br />delete [] a[i];<br />}<br />delete []a;<br />}</p><p>};</p><p>void main()<br />{<br />cholesky solution;<br />solution.cholesky_input();<br />solution.cholesky_decomposition();<br />solution.cholesky_output();<br />}</p><p>void cholesky::cholesky_input()<br />{<br />cout<<"輸入方程的個數";<br />cin>>n;<br />b=new double[n];<br />d=new double[n];<br />x=new double[n];<br />a=new double*[n];</p><p>for(i=0;i<n;i++)<br />{<br />a[i] = new double[n];<br />}</p><p>for(i=0;i<n;i++)<br />for(j=0;j<n;j++)<br />{<br />cout<<"/n輸入a["<<i<<"]["<<j<<"]=";<br />cin>>a[i][j];<br />}<br />for(i=0;i<n;i++)<br />for(j=0;j<n;j++)<br />{<br />if(a[i][j] != a[j][i])<br />{<br />cout<<"/n係數矩陣不對稱.失敗..."<<endl;<br />exit(0);<br />}<br />}<br />for(i=0;i<n;i++)<br />{</p><p>cout<<"/n輸入 b["<<i<<"]=";<br />cin>>b[i];<br />}<br />cout<<"/n輸入最小主要元素";<br />cin>>eps;//輸入段結束<br />}<br />void cholesky::cholesky_decomposition()<br />{<br />for(i=0;i<n;i++)<br />for(j=0;j<n;j++)<br />{<br />sum=a[i][j];<br />for(k=0;k<i;k++)<br />{<br />sum -= a[i][k]*a[j][k];<br />}<br />if( i==j)<br />{<br />if(sum <= 0)<br />{<br />cout<<"/n矩陣非正定.失敗..."<<endl;<br />exit(0);<br />}<br /> d[i]=sqrt(sum);<br />}<br />else<br />{<br />a[j][i] = sum/d[i];<br />}</p><p>}</p><p>for(i=0;i<n;i++)<br />{<br />sum=b[i];<br />for(k=0;k<i;k++)<br />{<br />sum -= a[i][k]*x[k];<br />}<br />x[i] = sum/d[i];<br />}</p><p>for(i=(n-1);i>=0;i--)<br />{<br />sum = x[i];<br />for(k=(i+1);k<n;k++)<br />{<br />sum -= a[k][i]*x[k];<br />}<br />x[i] = sum/d[i];<br />}<br />}<br />void cholesky::cholesky_output()<br />{<br />cout<<"/n:結果是："<<endl;<br />for(i=0;i<n;i++)<br />{<br />cout<<"x["<<i<<"]="<<x[i]<<endl;<br />}</p><p>}