经典问题用高斯约当算法求解线性方程组。这里要求对任意形式的线性方程组都能够妥善处理，不能只适用于方程个数和未知量数目相等的特殊情形。
先用循环结构将增广矩阵转换为阶梯形矩阵，循环结束时得到阶梯型矩阵非零行行数，同时得到一个链表其中存放有各非零行主元的列标，列标在链表中按从左到右的顺序依次递减。然后根据线性代数中线性方程组的解的情况及判别准则判断方程是否有解，有多少个解。当线性方程组有解时，需要用convert函数将其转换为简化行阶梯型矩阵，然后输出唯一解或一般解
更新:定义了新的有理数类,使得求解算法可以处理方程中常数项为有理数情形:
C语言代码如下：
#include <stdio.h>
#include <malloc.h>
#include <cmath>
#include <vector>
#include "RationalNumber.h"
using std::vector;
struct maincol
{
int col; //存放各主元下标的结构体类型
struct maincol* next;
};
typedef struct maincol mc1;
int test(int s, int t, vector<vector<RationalNumber>>& a, int& _j, const RationalNumber& zero); //判断增广矩阵的s行至M行与t列至N列相交形成的子矩阵是否为零矩阵，若是返回0，若不是返回第一个不为零的列的列标
void add(mc1* head, int col, mc1** tail); //函数，用于新建一节点，其中保存有主元col列标，然后按递减顺序将其插入maincol类型的链表
void convert(vector<vector<RationalNumber>>& a, int row, mc1* head); //函数，用于将阶梯型矩阵转化为简化行阶梯形矩阵
#define N 5 //增广矩阵列数
#define M 3 //增广矩阵行数
int main()
{
int matrix[M][N] = { {2, 1, -3, 5, 6}, {-3, 2, 1, -4, 5}, {-1, 3, -2, 1, 11} }; //增广矩阵
int i, j;
int row, col; //子矩阵左上角元素行列下标,循环结束后row用于存放阶梯形矩阵非零行行数
RationalNumber swap;
mc1* head, * tail, * psnew;
vector<vector<RationalNumber>> a(M, vector<RationalNumber>(N));
for (i = 0; i < M; ++i)
{
for (j = 0; j < N; ++j)
{
a[i][j] = RationalNumber(matrix[i][j], 1);
}
}
head = (mc1*)malloc(sizeof(mc1));
head->next = NULL;
tail = head;
row = col = 1; //子矩阵即为增广矩阵本身，用增广矩阵左上角元素行列标初始化s,t
RationalNumber zero(0, 1);
while ((col = test(row, col, a, j, zero)) != 0) //子矩阵不为零矩阵
{
if (row == M) //增广矩阵已化为阶梯形矩阵
{
add(head, col, &tail); //最后一个非零行主元列标放入maincol类型链表
break; //结束循环
}
else
{
int temp_j = j;
for (i = j + 1; i < M; i++)
{
if (a[j][col - 1].abs() < a[i][col - 1].abs()) //列选主元
j = i;
}
if (row - 1 != j)
{
for (i = col - 1; i < N; i++)
{
swap = a[j][i];
a[j][i] = a[row - 1][i]; //列选主元
a[row - 1][i] = swap;
}
}
if (col == N) //增广矩阵已经化为阶梯形矩阵
{
add(head, col, &tail); //最后一个非零行主元列标放入maincol类型链表
break; //结束循环
}
if (temp_j == row - 1)
++temp_j;
for (j = temp_j; j < M; j++)
{
if (a[j][col - 1] != zero)
{
a[j][col - 1] = -(a[j][col - 1] / a[row - 1][col - 1]);
for (i = col; i < N; i++) //消元
{
if (a[row - 1][i] != zero)
{
a[j][i] = a[j][col - 1] * a[row - 1][i] + a[j][i];
}
}
}
}
add(head, col, &tail); //将消元后得到的主元列标col放入maincol类型链表
++row;
++col; //更新s,t,使s,t成为消元后得到的新的子矩阵的左上角元素行列标,为test函数操作新的子矩阵作准备
}
}
if (col == 0) //从循环控制条件退出循环
--row; //此时增广矩阵已成为阶梯形矩阵,非零行函数就是row-1
if (row == 0) //利用线性方程组解的判别准则判断是否有解，有多少个解
{
printf("线性方程组有无穷多组解\n"); //增广矩阵为零矩阵，无穷多组解
printf("一般解:\n");
for (i = 1; i < N; i++)
printf("x%d=t%d\n", i, i); //输出解
}
else
{
psnew = head->next;
if (psnew->col == N) //阶梯形矩阵最后一主元在最后一列，无解
printf("线性方程组无解\n");
else
{
convert(a, row, head); //用convert函数将阶梯形矩阵进一步化为简化行阶梯形矩阵
if (row == N - 1) //非零行行数等于未知量个数，有唯一解
{
printf("线性方程组有唯一解:\n");
for (i = 1; i <= row; i++)
std::cout << "x" << i << "=" << a[i - 1][N - 1] << std::endl;//输出唯一解
}
else //非零行行数小于未知量个数，有无穷多组解
{
int* m;
m = (int*)malloc((N - 1 - row) * sizeof(int));
i = N - 1 - row;
for (j = N - 1; j >= 1; j--)
{
if (j != psnew->col)
{
m[--i] = j; //从所有未知量标号中筛选出自由未知量标号
if (i == 0)
break;
}
else
{
if (psnew->next != NULL)
psnew = psnew->next;
}
}
printf("线性方程组有无穷多组解\n");
printf("一般解:\n");
i = row;
for (psnew = head->next; psnew != NULL; psnew = psnew->next)
{
std::cout << "x" << psnew->col << "=" << a[i - 1][N - 1]; //输出一般解
for (j = 0; j < N - 1 - row; j++)
{
if (m[j] < psnew->col)
{
printf("-%dx%d", 0, m[j]);
}
else
{
std::cout << "-" << a[i - 1][m[j] - 1] << "x" << m[j];
}
}
printf("\n");
i--;
}
}
}
}
return 0;
}
int test(int s, int t, vector<vector<RationalNumber>>& a, int& _j, const RationalNumber &zero) //判断增广矩阵的s行至M行与t列至N列相交形成的子矩阵是否为零矩阵，若是返回0，若不是返回第一个不为零的列的列标
{
int i, j;
for (j = t - 1; j < N; j++)
{
for (i = s - 1; i < M; i++)
{
if (a[i][j] != zero)
{
_j = i;
return (j + 1);
}
}
}
return (0);
}
void add(mc1* head, int col, mc1** tail) //函数，用于新建一节点，其中保存有主元col列标，然后按递减顺序将其插入maincol类型的链表
{
mc1* psnew;
psnew = (mc1*)malloc(sizeof(mc1));
psnew->col = col;
if (head->next == NULL)
{
psnew->next = NULL;
head->next = psnew;
*tail = psnew;
}
else
{
psnew->next = head->next;
head->next = psnew;
}
}
void convert(vector<vector<RationalNumber>> &a, int row, mc1* head) //函数，用于将阶梯型矩阵转化为简化行阶梯形矩阵
{
mc1* psnew, * pq;
int i, j, k, m;
RationalNumber one(1, 1);
RationalNumber zero(0, 1);
psnew = head->next;
for (i = row - 1; i >= 0; i--)
{
if (a[i][psnew->col - 1] != one) //各非零行主元化为1
{
for (j = psnew->col; j < N; j++)
a[i][j] = a[i][j] / a[i][psnew->col - 1];
}
psnew = psnew->next;
}
psnew = head->next; //向上消元把除第一个主元外其余主元所在列中在主元上方的部分变为零
for (i = row - 1; i >= 1; i--)
{
m = N - psnew->col - (row - i - 1); //获取未知量标号1,2,--,N-1中位于i+1非零行主元列标号右边的自由未知量标号个数
pq = head->next; //pq指向存放最后一个非零行主元列标号的节点
for (k = N; k > psnew->col; k--)
{
if (k != pq->col)
{
if (a[i][k - 1] != zero)
{
for (j = i - 1; j >= 0; j--)
{
if (a[j][psnew->col - 1] != zero)
{
a[j][k - 1] = -(a[i][k - 1] * a[j][psnew->col - 1]) + a[j][k - 1]; //从右向左作初等行变换直至i+1行主元所在列右边的列位置，期间跳过i+2----row行主元所在的列
}
}
}
m--;
if (m == 0)
break;
}
else
{
if (pq->next != psnew)
pq = pq->next;
}
}
psnew = psnew->next; //递进至上一行主元，为新一轮向上消元作准备
}
}
RationalNumber.h内容如下:
#pragma once
#include <iostream>
#include <cmath>
#include <algorithm>
class RationalNumber
{
public:
RationalNumber(long long m, long long d);
RationalNumber() = default;
long long computeGCD();
bool operator<(const RationalNumber& opNum);
bool operator<=(const RationalNumber& opNum);
bool operator>(const RationalNumber& opNum);
bool operator>=(const RationalNumber& opNum);
bool operator==(const RationalNumber& opNum);
bool operator!=(const RationalNumber& opNum);
RationalNumber operator+(const RationalNumber& opNum);
RationalNumber operator-(const RationalNumber& opNum);
RationalNumber operator*(const RationalNumber& opNum);
RationalNumber operator/(const RationalNumber& opNum);
RationalNumber operator-() { RationalNumber t = *this; t.positive = !t.positive; return t; }
RationalNumber abs() { RationalNumber t = *this; t.positive = true; return t; }
public:
bool positive;
long long molecule;
long long denominator;
};
std::ostream& operator<<(std::ostream& out, const RationalNumber& output)
{
if (output.positive == false && output.molecule != 0)
out << "-";
if (output.molecule == 0)
out << "0";
else
{
if (output.denominator == 1)
out << output.molecule;
else
out << output.molecule << "/" << output.denominator;
}
return out;
}
bool RationalNumber::operator<(const RationalNumber& opNum)
{
if (molecule == 0)
{
if (opNum.molecule == 0)
{
return false;
}
else
{
if (opNum.positive)
return true;
return false;
}
}
else
{
if (opNum.molecule == 0)
{
if (positive == false)
return true;
return false;
}
else
{
if (!positive && opNum.positive)
{
return true;
}
else if (positive && !opNum.positive)
{
return false;
}
else if (positive)
{
if (denominator == opNum.denominator)
{
return molecule < opNum.molecule;
}
else
{
if (molecule == opNum.molecule)
return denominator > opNum.denominator;
if (molecule < opNum.molecule && denominator > opNum.denominator)
return true;
if (molecule > opNum.molecule && denominator < opNum.denominator)
return false;
return molecule * opNum.denominator < opNum.molecule* denominator;
}
}
else
{
if (denominator == opNum.denominator)
{
return molecule > opNum.molecule;
}
else
{
if (molecule == opNum.molecule)
return denominator < opNum.denominator;
if (molecule < opNum.molecule && denominator > opNum.denominator)
return false;
if (molecule > opNum.molecule && denominator < opNum.denominator)
return true;
return molecule * opNum.denominator > opNum.molecule * denominator;
}
}
}
}
}
bool RationalNumber::operator<=(const RationalNumber& opNum)
{
return !operator>(opNum);
}
bool RationalNumber::operator>(const RationalNumber& opNum)
{
if (molecule == 0)
{
if (opNum.molecule == 0)
{
return false;
}
else
{
if (opNum.positive)
return false;
return true;
}
}
else
{
if (opNum.molecule == 0)
{
if (positive == false)
return false;
return true;
}
else
{
if (!positive && opNum.positive)
{
return false;
}
else if (positive && !opNum.positive)
{
return true;
}
else if (positive)
{
if (denominator == opNum.denominator)
{
return molecule > opNum.molecule;
}
else
{
if (molecule == opNum.molecule)
return denominator < opNum.denominator;
if (molecule < opNum.molecule && denominator > opNum.denominator)
return false;
if (molecule > opNum.molecule && denominator < opNum.denominator)
return true;
return molecule * opNum.denominator > opNum.molecule * denominator;
}
}
else
{
if (denominator == opNum.denominator)
{
return molecule < opNum.molecule;
}
else
{
if (molecule == opNum.molecule)
return denominator > opNum.denominator;
if (molecule < opNum.molecule && denominator > opNum.denominator)
return true;
if (molecule > opNum.molecule && denominator < opNum.denominator)
return false;
return molecule * opNum.denominator < opNum.molecule* denominator;
}
}
}
}
}
bool RationalNumber::operator>=(const RationalNumber& opNum)
{
return !operator<(opNum);
}
bool RationalNumber::operator==(const RationalNumber& opNum)
{
if (molecule == 0)
{
if (opNum.molecule == 0)
{
return true;
}
else
{
return false;
}
}
else
{
if (opNum.molecule == 0)
{
return false;
}
else
{
if (!positive && opNum.positive || positive && !opNum.positive)
{
return false;
}
else
{
if (denominator == opNum.denominator)
{
return molecule == opNum.molecule;
}
else if (molecule == opNum.molecule)
{
return false;
}
else
{
return molecule * opNum.denominator == opNum.molecule * denominator;
}
}
}
}
}
bool RationalNumber::operator!=(const RationalNumber& opNum)
{
return !operator==(opNum);
}
RationalNumber RationalNumber::operator/(const RationalNumber& opNum)
{
if (opNum.molecule == 0)
{
std::cout << "Error,除数不能为零!" << std::endl;
exit(-1);
}
if (molecule == 0)
return RationalNumber(0, 1);
RationalNumber temp = opNum;
std::swap(temp.molecule, temp.denominator);
return operator*(temp);
}
RationalNumber RationalNumber::operator*(const RationalNumber& opNum)
{
if (molecule == 0 || opNum.molecule == 0)
{
return RationalNumber(0, 1);
}
if (positive && !opNum.positive || !positive && opNum.positive)
{
return RationalNumber(-(molecule * opNum.molecule), denominator * opNum.denominator);
}
else
{
return RationalNumber(molecule * opNum.molecule, denominator * opNum.denominator);
}
}
RationalNumber RationalNumber::operator-(const RationalNumber& opNum)
{
RationalNumber temp = opNum;
temp.positive = !temp.positive;
return operator+(temp);
}
RationalNumber RationalNumber::operator+(const RationalNumber& opNum)
{
if (molecule == 0)
{
return opNum;
}
if (opNum.molecule == 0)
{
return *this;
}
if (denominator == opNum.denominator)
{
if (positive)
{
if (opNum.positive)
{
return RationalNumber(molecule + opNum.molecule, denominator);
}
else
{
return RationalNumber(molecule - opNum.molecule, denominator);
}
}
else
{
if (opNum.positive)
{
return RationalNumber(opNum.molecule - molecule, denominator);
}
else
{
return RationalNumber(-(opNum.molecule + molecule), denominator);
}
}
}
if (positive)
{
if (opNum.positive)
{
return RationalNumber(molecule * opNum.denominator + opNum.molecule * denominator, denominator * opNum.denominator);
}
else
{
return RationalNumber(molecule * opNum.denominator - opNum.molecule * denominator, denominator * opNum.denominator);
}
}
else
{
if (opNum.positive)
{
return RationalNumber(opNum.molecule * denominator - molecule * opNum.denominator, denominator * opNum.denominator);
}
else
{
return RationalNumber(-(opNum.molecule * denominator + molecule * opNum.denominator), denominator * opNum.denominator);
}
}
}
long long RationalNumber::computeGCD()
{
long long d = denominator;
long long m = molecule;
while (m != 0)
{
long long temp = m;
m = d % m;
d = temp;
}
return d;
}
RationalNumber::RationalNumber(long long m, long long d)
{
if (d == 0)
{
std::cout << "ERROR,分母不能为零!" << std::endl;
exit(-1);
}
if (m < 0 && d>0 || m > 0 && d < 0)
{
positive = false;
}
else
{
positive = true;
}
denominator = std::abs(d);
molecule = std::abs(m);
if (molecule != 0)
{
if (molecule > 1 && denominator > 1)
{
unsigned long long gcd = computeGCD();
denominator /= gcd;
molecule /= gcd;
}
}
}
扫一扫
1727
被折叠的 条评论
为什么被折叠?
到【灌水乐园】发言
