废话不说,直接上代码。计算AX=b线性或超定方程组。原理就是将A分解成上三角和下三角阵,网上资料很多,这里就不贴了,直接上代码。代码都是经过多次测试的,请放心使用。
void LUD(double * a, int n, double * L, double * U, double * D)
{
for (int i = 0; i < n; ++i)
{
//先算U的第0~n-1列
for (int j = i; j < n; ++j)
{
double sum = 0;
//mulmatrix(L, n, n, U, n, n, i, j, &sum);
for (int k = 0; k < i; ++k)
sum += L[i * n + k] * U[k * n + j];
U[i * n + j] = a[i * n + j] - sum;
}
D[i] = U[i * n + i];
//再算j的第i列,
for (int j = i + 1; j < n; ++j) //这里j表示行,i表示列
{
double sum = 0;
for (int k = 0; k < i; ++k)
sum += L[j * n + k] * U[k * n + i];
L[j * n + i] = (a[j * n + i] - sum) / U[i * n + i];
}
L[i * n + i] = 1.;
}
}
//n是阶数
bool Improve_Square(double * a, int n, double * b, double x[])
{
double * L = (double *)malloc(sizeof(double) * n * n);
double * U = (double *)malloc(sizeof(double) * n * n);
memset(L, 0, sizeof(double) * n * n);
memset(U, 0, sizeof(double) * n * n);
double * D = (double *)malloc(sizeof(double) * n);
double * y = (double *)malloc(sizeof(double) * n);
if (L == NULL || D == NULL || y == NULL)
{
free(L); L = NULL;
free(D); D = NULL;
free(y); y = NULL;
free(U); U = NULL;
return false;
}
//if (!LDLT(a, n, L, D))
// return false;
LUD(a, n, L, U, D);
y[0] = b[0];
for (int k = 1; k < n; ++k)
{
double sum = 0;
for (int j = 0; j < k; ++j)
sum += L[k * n + j] * y[j];
y[k] = b[k] - sum;
}
//计算D*transpose(L)*x = y
x[n - 1] = y[n - 1] / D[n - 1];
for (int k = n - 2; k >= 0; --k)
{
double sum = 0;
for (int j = k + 1; j < n; ++j)
sum += L[j * n + k] * x[j];
x[k] = y[k] / D[k] - sum;
}
free(L); L = NULL;
free(D); D = NULL;
free(y); y = NULL;
free(U); U = NULL;
return true;
}
LUD是修正平方根法的核心,Improve_Square就是解算线性方程组的主体。Improve_Square形参的a和b并非原始的A和b,
a = AT*A; AT表示A的转置
b = AT*b;
x提前申请好空间,以便放入最终计算结果。
注意,AT*A一定是对称阵,但不一定是正定阵,这就是需要判断各阶逐子行列式是不是都>0。满足正定,解才是正确的,否则就是奇异解,不能使用。
如果不会c/c++下怎么做A的转置和矩阵乘法,请参考下面代码。
//矩阵乘法
bool mulmatrix(double * a, int ar, int ac, double * b, int br, int bc, double * c)
{
//c为输出矩阵
if (ac != br)
return false;
else
{
for (int i = 0; i < ar; ++i)
{
for (int j = 0; j < bc; ++j)
{
double sum = 0;
for (int k = 0; k < ac; ++k)
sum += a[i * ac + k] * b[k * bc + j];
c[i*bc + j] = sum;
}
}
return true;
}
}
//矩阵转置
void transpose(double * coef, int row, int col, double * trans_coef)
{
for (int r = 0; r < row; ++r)
for (int c = 0; c < col; ++c)
trans_coef[c * row + r] = coef[r * col + c];
}
有任何问题请下发留言或发360759935@qq.com邮箱,明天更新QR分解计算线性或超定方程组的代码文档,请大家多多支持。