* LU Decomposition
源自於
http://www.stat.nctu.edu.tw/misg/SUmmer_Course/C_language/Ch06/LUdecomposition.htm
將矩陣 A 分解成一上三腳矩陣 U 與一下三角矩陣 L所以方程式為:
令
則 x 由以下兩個方程式求解:
先算出b',再代回算出x.
(1)矩陣的分解
4X4的矩陣為例 子:
將 L 矩陣每一列與與 U 矩陣 的1,2,3,4行相乘則得到:
L11 = a11 , L21 = a21, L31 = a31, L41 = a41
將 L 矩陣的第一列和 U 矩陣的2,3,4行相乘則得到:
同理 U 的第二行乘以 L 的2,3,4列可得:
同理可得:
其一般化的式子為:
當 j=1 時, Li1 = ai1
當 i=1 時,
LU 分解法受歡迎的原因是不需要儲存 LU 中的 0 與 1,原矩陣A能表示LU兩個矩陣:
(2)演算法
(a) for i = 1 to n :do
L[i,1]=a[i,1]
(b) for j = 1 to n: do
U[1,j]=a[1,j]/L[1,1]
End i , j.
(c) for j = 2 to n : do
for i = j to n : do
sum=0.0
for k = 1 to j-1 : do
sum = sum + L[i,k]*U[k,j]
End k.
L[i,j] = a[i,j] - sum
End i.
U[j,j] = 1.
for i = j+1 to n : do
sum=0.0
for k = 1 to j-1 : do
sum = sum + L[j,k]*U[k,i]
End k.
U[j,i] = (A[j,i] - sum)/L[j,j]
End i.
End j.
(d)反代
(3)程式
#include <stdio.h>
#include <stdlib.h>
#include <math.h>
#define dim 10
double a[dim][dim], L[dim][dim], U[dim][dim], x[dim];
double b[dim], bp[dim];
mprint( int , int , double [][dim] , double [] );
main()
{
int n,i,j,k;
double temp;
printf("input the dimension of matrix\n");
scanf("%d",&n);
if(n>=dim){
printf("dimension over preset value %d\n",dim);
exit(0);
}
printf("Input the a matrix\n\n");
for(i=1; i<=n; i++){
for(j=1;j<=n;j++){
printf("input the element of a[%d,%d]\n",i,j);
scanf("%lf",&a[i][j]);
}
}
printf("\n Input the b vector\n\n");
for(i=1; i<=n ; i++){
printf("input the value of b[%d]\n",i);
scanf("%lf",&b[i]);
}
mprint(n,n,a,b);
/* start to the procedure of LU decomposition */
for(i=0; i<dim; i++){
for(j=0; j<dim;j++){
L[i][j]=0.0; U[i][j]=0.0;
}
}
for(i=1; i<=n; i++){L[i][1]=a[i][1]; }
for(j=1; j<=n;j++){U[1][j]=a[1][j]/L[1][1];}
for(j=2; j<=n; j++){
for(i=j; i<=n;i++){
temp=0.0;
for(k=1;k<=j-1;k++){
temp=temp+L[i][k]*U[k][j];
}
L[i][j]=a[i][j]-temp;
}
U[j][j]=1.0;
for(i=j+1;i<=n;i++){
temp=0.0;
for(k=1; k<=j-1;k++){
temp=temp+L[j][k]*U[k][i];
}
U[j][i]=(a[j][i]-temp)/L[j][j];
}
}
/* Backward substitution */
bp[1]=b[1]/L[1][1];
for(i=2; i<=n; i++){
temp=0.0;
for(k=1; k<=i-1;k++){
temp=temp+L[i][k]*bp[k];
}
bp[i]=(b[i]-temp)/L[i][i];
}
x[n]=bp[n];
for(j=n-1; j>=1; j--) {
temp=0.0;
for(k=j+1; k<=n; k++){
temp=temp+U[j][k]*x[k];
}
x[j]=bp[j]-temp;
}
mprint(n,n,l,b);
mprint(n,n,u,bp);
for(j=1; j<=n; j++){
printf(" x[%d]=%10.5lf\n",j,x[j]);
}
}/*End of main() */
mprint(int m,int n,double aa[][dim],double bb[])
{
int i,j;
for(i=1; i<=m; i++){
for(j=1; j<=n;j++){
printf("%10.5e ",aa[i][j]);
}
printf(" %10.5e \n",bb[i]);
}
}
沒有留言:
張貼留言