1、三角分解法基本思想
当线性方程组AX=b时:
- 先将非奇异矩阵A进行LU(这里指Doolittle分解)分解,使A=LU。则方程化为LUX=b。
- 从而将问题转化为两个简单地三角方程组:
LY=b 求解Y
UX=Y 求解X
2、Doolittle分解步骤
3、代码
#include<iostream>
#include<math.h>
#include<iomanip>
using namespace std;
//Doolittle分解法
void Doolittle(double** a, double* b, int n)
{
double** u = new double* [n];
for (int i = 0; i < n; i++)
u[i] = new double[n];
double** l = new double* [n];
for (int i = 0; i < n; i++)
l[i] = new double[n];
double* jie_x = new double[n];
double* jie_y = new double[n];
for (int k = 1; k < n + 1; k++)
{
for (int j = k, i = k + 1; j < n + 1, i < n + 1; j++, i++)
{
double sum_1 = 0;
for (int t = 1; t < k; t++)
{
sum_1 += l[k - 1][t - 1] * u[t - 1][j - 1];
}
u[k - 1][j - 1] = a[k - 1][j - 1] - sum_1;
double sum_2 = 0;
for (int t = 1; t < k; t++)
{
sum_2 += l[i - 1][t - 1] * u[t - 1][k - 1];
}
l[i - 1][k - 1] = (a[i - 1][k - 1] - sum_2) / u[k - 1][k - 1];
}
double sum_1 = 0;
for (int t = 1; t < k; t++)
{
sum_1 += l[k - 1][t - 1] * u[t - 1][n - 1];
}
u[k - 1][n - 1] = a[k - 1][n - 1] - sum_1;
}
jie_y[0] = b[0];
for (int i = 2; i < n + 1; i++)
{
double sum = 0;
for (int t = 1; t < i; t++)
sum += l[i - 1][t - 1] * jie_y[t - 1];
jie_y[i - 1] = b[i - 1] - sum;
}
jie_x[n - 1] = jie_y[n - 1] / u[n - 1][n - 1];
for (int i = n - 1; i > 0; i--)
{
double sum = 0;
for (int t = i + 1; t < n + 1; t++)
sum += u[i - 1][t - 1] * jie_x[t - 1];
jie_x[i - 1] = (jie_y[i - 1] - sum) / u[i - 1][i - 1];
}
cout << "使用Doolittle分解法的解为\n";
for (int i = 0; i < n - 1; i++)
cout << "x" << i + 1 << " = " << jie_x[i] << endl;
cout << "x" << n << " = " << jie_x[n - 1] << endl;
cout << endl;
delete[]jie_x;
delete[]jie_y;
for (int i = 0; i < n; i++)
delete[]u[i];
delete[]u;
for (int i = 0; i < n; i++)
delete[]l[i];
delete[]l;
}
void main()
{
//输入矩阵A和b
int n;
cout << "请输入系数矩阵的行数(或列数,列数等于行数):\n";
cin >> n;
cin.clear();//这3行代码重置输入
while (cin.get() != '\n')
continue;
cout << "请输入系数矩阵A:\n";
//要用到二维动态数组
double** a;
a = new double* [n]; //申请行的空间
for (int i = 0; i < n; i++)
a[i] = new double[n];
for (int i = 0; i < n; i++)
{
cout << "第 " << i + 1 << " 行:";
for (int j = 0; j < n; j++)
while (!(cin >> a[i][j]))
{
cin.clear();
while (cin.get() != '\n')
continue;
cout << "请输入数字";
}
cin.clear();
while (cin.get() != '\n')
continue;
}
cout << "请输入矩阵b:\n";
double* b = new double[n];
for (int j = 0; j < n; j++)
cin >> b[j];
cin.clear();
while (cin.get() != '\n')
continue;
//将系数矩阵A与矩阵b存储在临时矩阵中,临时矩阵中的元素参与行列变换
double** temp_a = new double* [n];
double* temp_b = new double[n];
for (int i = 0; i < n; i++)
temp_a[i] = new double[n];
for (int i = 0; i < n; i++)
for (int j = 0; j < n; j++)
temp_a[i][j] = a[i][j];
for (int i = 0; i < n; i++)
temp_b[i] = b[i];
Doolittle((double**)temp_a, temp_b, n);
for (int i = 0; i < n; i++)
delete[]temp_a[i];
delete[]temp_a;
delete[]temp_b;
for (int i = 0; i < n; i++)
delete[]a[i];
delete[]a;
delete[]b;
}
待博主有空后,再将上述基础知识补齐。