数值分析——三角矩阵排序向量及求解过程

三角矩阵排序向量及求解过程

问题

1. 对于一个三角矩阵，如何得到排序向量 $p$ ？给出算法并用程序实现它；
2. 请实现三角系统的向前带入算法，给出程序并通过算例验证；
3. 请给出三角系统的向后带入算法，给出程序并通过算例验证；

分析

1. 由于给出矩阵默认为三角矩阵，因此通过排序后只有上三角或下三角两种可能。因此我们先以下三角矩阵为例进行分析，然后再考虑如何判断三角矩阵的类型；
2. 对于一个可化为下三角矩阵的三角矩阵$A_{n\times n}$，其第 $i$ 列向量至少有 $（n-i+1）$ 个 $0$ 元。而当我们从后向前遍历时，第 $i$ 列出现的非零元所在行向量即为排序后第 $i$ 行（需排除从 $i+1$ 到 $n$ 已找出的行向量），上三角矩阵反之；
3. 通过2所述方法我们可以找到排序向量，但倘若采用暴力遍历的方式，其复杂度为 $O(n^2)$ 对于大矩阵来说是不太理想的，考虑采用队列的方式优化，即用 $(1\to n)$ 到初始化队列，然后按队列里数据 $i$ 从矩阵最后一列开始遍历，第 $i$ 行数值非零即弹出否则再入队尾；
4. 由于不清楚三角阵的类型，我们可以先假设其为下三角阵，如果该矩阵本为上三角阵或存在主元为0，则迭代过程中就会出现：在排除从 $i+1$ 到 $n$ 已找出的行向量后剩余行向量在第 $i$ 列的值均为0的情况，即找不出非零元。那么可以假设其为上三角阵再计算一遍。若仍存在找不到非零元情况，则可以判断存在0主元排序向量不唯一且向前或向后带入法不可解；
5. 最后考虑到在c++中输入矩阵可视化较差且出错不易修改，考虑在txt文件中输入矩阵，然后通过算法读入并存储。

实验代码

#include<cstdio>
#include<vector>
#include<fstream>
#include<iostream>
#include<string>
#include<queue>
#include <algorithm>
using namespace std;

bool under_tag(vector<vector<int>>& matrix, vector<int>& p, int n) {
															//假设为下三角阵
	queue<int> init;										//带检测向量行编号
	for (int i = 0; i < n; i++) {
		init.push(i);
	}
	int l = n - 1;											//从最后一列开始遍历
	int rept = 0;											//迭代步
	while (!init.empty() && rept < n) {						//迭代步大于n,则存在0主元或为上三角阵
		if (matrix[init.front()][l] != 0) {					//避免全遍历而采用的加速
			p.push_back(init.front());
			init.pop();
			l--;
		}
		else {
			rept++;
			init.push(init.front());
			init.pop();
		}
	}
	if (rept < n) { return true; }
	else {
		while (!p.empty())	p.pop_back();					//清空
		return false;
	}
}

bool up_tag(vector<vector<int>>& matrix, vector<int>& p, int n) {
															//假设为上三角阵
	queue<int> init;										//带检测向量行编号
	for (int i = 0; i < n; i++) {
		init.push(i);
	}
	int l = 0;												//从第一列开始遍历
	int rept = 0;											//迭代步
	while (!init.empty() && rept < n) {						//迭代步大于n,则存在0主元或为下三角阵
		if (matrix[init.front()][l] != 0) {					//避免全遍历而采用的加速
			p.push_back(init.front());
			init.pop();
			l++;
		}
		else {
			rept++;
			init.push(init.front());
			init.pop();
		}
	}
	if (rept < n) { return true; }
	else {
		while (!p.empty())	p.pop_back();					//清空
		return false;
	}
}

void function1(vector<double>& ans, vector<vector<int>>& matrix, vector<int>& p, int n) {	//向前带入算法
	for (int i = 0; i < n; i++) {
		if (i == 0) {
			ans[i] = (double)matrix[p[i]][n] / matrix[p[i]][i];
		}
		else {
			double temp = 0;
			for (int j = 0; j < i; j++) {
				temp += (double)matrix[p[i]][j] * ans[j];
			}
			ans[i] = ((double)matrix[p[i]][n] - temp) / matrix[p[i]][i];
		}
	}
}

void function2(vector<double>& ans, vector<vector<int>>& matrix, vector<int>& p, int n) {	//向后带入算法
	for (int i = n - 1; i >= 0; i--) {
		if (i == n - 1) {
			ans[i] = (double)matrix[p[i]][n] / matrix[p[i]][i];
		}
		else {
			double temp = 0;
			for (int j = i + 1; j < n; j++) {
				temp += (double)matrix[p[i]][j] * ans[j];
			}
			ans[i] = ((double)matrix[p[i]][n] - temp) / matrix[p[i]][i];
		}
	}
}

int main() {
	int n;
	ifstream read_file("data1.txt");				//读入矩阵数据，注：数据需用空格分隔，且每行最后一个数据后不加空格
	if (!read_file.is_open()) {
		cout << "can't open this file" << endl;		//若文件读取失败，请核对是否将文件放入工程文件所在目录且格式为可读
		return 0;
	}
	string	line;
	vector<string> arr;
	while (getline(read_file, line)) {
		line += " ";										//统一格式，以数据+空格形式分割
		arr.push_back(line);
	}

	n = arr.size();											//获得矩阵的行数，即维度
	vector<vector<int>> matrix(n, vector<int>(n + 1, 0));	//用二维vector存储增广矩阵
	for (int i = 0; i < n; i++) {
		int m = arr[i].size();
		int temp = 0;
		int line = 0;
		for (int j = 0; j < m; j++) {						//将初始使用string存储的每行数据，分割转化为独立的数据
			if (arr[i][j] != ' ') {
				temp = temp * 10 + arr[i][j] - '0';
			}
			else {
				matrix[i][line] = temp;
				line++;
				temp = 0;
			}
		}
	}

	cout << "输入矩阵为：" << endl;							//将读入矩阵输出
	for (int i = 0; i < n; i++) {
		for (int j = 0; j < n + 1; j++) {
			cout << matrix[i][j] << " ";
		}
		cout << endl;
	}
	cout << endl;



	vector<int> p;											//定义排序向量
	while (!p.empty())	p.pop_back();						//初始化
	vector<double> ans;										//存储解并初始化
	for (int i = 0; i < n; i++) ans.push_back(0);

	if (under_tag(matrix, p, n) == true) {					//矩阵为下三角阵
		reverse(p.begin(), p.end());						//原p序列的倒序即为排序向量
		cout << "排序向量P为：";
		for (int i = 0; i < p.size(); i++) {
			cout << p[i] + 1 << " ";
		}
		cout << endl << endl;
		cout << "采用向前带入算法" << endl << endl;				//采用向前带入算法
		function1(ans, matrix, p, n);
	}
	else if (up_tag(matrix, p, n) == true) {				  //矩阵为上三角阵
		cout << "排序向量P为：";
		for (int i = 0; i < p.size(); i++) {
			cout << p[i] + 1 << " ";
		}
		cout << endl << endl;
		cout << "采用向后带入算法" << endl << endl;				//采用向后带入算法
		function2(ans, matrix, p, n);
	}
	else {
		cout << "存在主元为0，排序向量不唯一,,且向前或向后带入法不可解" << endl;
		return 0;
	}

	cout << "方程的解x1到x" << n << "为：";					//输出解
	for (int i = 0; i < n; i++) {
		cout << ans[i] << " ";
	}
	cout << endl;
	return 0;
}

结语

最近好忙啊，没时间再作修改了，就按照上机报告的样子拷贝过来了。有什么错误或考虑不周全的忘各位看官们指出。另外要是有人能教教我Typora代码块导出错位如何解决就好了，哈哈哈，不过我这小博客应该没几个人看吧~