选主元的三角分解法

2012-2013（1）专业课程实践论文选主元（列）的三角分解法

范俊，0818180124，R数学08-1班

毛勇，0818180117，R数学08-1班

一、算法理论

从直接三角分解公式可看出，当时计算将中断或者当绝对值很小时，按分解公式计算可能要求舍入误差的累积，但如果非奇异，就可通过交换的行实现矩阵的分解，因此可采用与列主元素消去法类似的方法

（可以证明下述方法与列主元素消去法等价），将直接三角分解法修改为（部分）选主元的三角分解法。

设第步分解已完成，这时有

第步分解需要用到式

及式

为了避免用小的数作除数，引进量

于是有

,,

用作为，交换的行与行元素(将位置的新元素仍记作及

)，于是有。由此在进行第步分解计算。

该程序即利用选主元的三角分解法解方程求方程的根。先选择列主元，再构造LU矩阵，通过求解LY = B和UX = Y得出需要的解。

注：方程维数在程序中需按题意自定义。

否

源代码源代码：

LU_Decomposition.cpp

#include

#include

#define N 4 //矩阵维数，可自定义

static double A[N][N]; //系数矩阵

static double B[N]; //右端项

static double Y[N]; //中间项

static double X[N]; //输出

static double S[N]; //选取列主元的比较器

int i,j,k; //计数器

void main()

{

cout << "请输入线性方程组（ai1，ai2，ai3......ain, yi）："<

for ( i = 0; i < N ;i++)

{

for (int j = 0; j< N ;j++ )

cin >> A[i][j];

cin >>B[i];

}

for (k = 0 ; k < N; k++)

{

//选列主元

int index =k;

for (i = k ; i < N; i++)

{

double temp = 0;

for (int m = 0 ; m < k ;m++)

{

temp = temp + A[i][m]*A[m][k];

}

S[i] = A[i][k] - temp;

if(S[index] < S[i])

{

index = i;

}

} //交换行

double temp;

for( i = k ; i < N ;i++ )

{

temp = A[index][i];

A[index][i] = A[k][i];

A[k][i] = temp;

}

temp = B[index];

B[index] = B[k];

B[k] = temp; // 构造L、U矩阵

for (j = k ; j < N; j++)

{

double temp = 0;

for (int m = 0 ;m < k; m++ )

{

temp = temp + A[k][m] * A[m][j];

}

A[k][j] = A[k][j] - temp; //先构造U一行的向量}

for( i = k+1; i < N; i++)

{

double temp = 0;

for (int m =0 ; m < k; m++ )

{

temp = temp + A[i][m] * A[m][k];

}

A[i][k] = (A[i][k] - temp)/A[k][k]; //再构造L一列的向量}

}

//求解L Y = B

Y[0] = B[0];

for (i = 1; i < N ; i++)

{

double temp = 0;

for (int j =0 ; j < i; j++ )

{

temp = temp + A[i][j] * Y[j];

}

Y[i] = B[i] - temp;

}

//求解UX = Y

X[N-1] = Y[N-1]/A[N-1][N-1];

for (i = N-2 ; i >= 0 ; i-- )

{

double temp = 0;

for (int j =i+1 ; j < N; j++ )

{

temp = temp + A[i][j] * X[j];

}

X[i] = (Y[i] - temp)/A[i][i];

}

//打印X

cout << "线性方程组的解（X1，X2，X3......Xn）为："<

{

cout << X[i] <<" ";

}

}

四、算法实现

例1. 用选主元（列）的三角分解法解

解：利用软件的解答方法为：

（1）先将程序中的维数N改为3。

（2）输入1 2 3 14回车，输入2 5 2 18 回车，输入3 1 5 20回车。

（3）显示结果 1 2 3

例2. 用选主元（列）的三角分解法解

解：利用软件的解答方法为：

（1）先将程序中的维数N改为3。

（2）输入0 3 4 1回车，输入1 -1 1 2回车，输入2 1 2 3回车。

（3）显示结果 1.27273 -1.18182 0.818182