eigenvalue / eigenvector

gistfile1.cpp

#include <iostream>
#include <vector>
#include <math.h>
#include <string>
#include <iomanip>
using namespace std;

void sub_into(vector< vector<double>> A, double r);
void gaussian_elimination(vector< vector<double>> *matrix, int m);

int main(){
	//輸入 2 * 2 矩陣
	int m = 2, n = 2;
	vector< vector<double>> A;
	A.resize(m);
	for (int r = 0; r < m; r++){
		cout << "row " << (r + 1) << ":";
		A[r].resize(n);
		for (int c = 0; c < n; c++){
			cin >> A[r][c];
		}
	}

	//顯示計算過程
	cout << "λI-A = { λ" << (-1 * A[0][0] > 0 ? "+" : "") << -1 * A[0][0] << "," << -1 * A[0][1] << ";" << endl;
	cout << "          " << -1 * A[1][0] << ",λ" << (-1 * A[1][1] > 0 ? "+" : "" ) <<  - 1 * A[1][1] << "}" << endl;
	cout << endl;

	double a = 1;
	double b = -1 * (A[0][0] + A[1][1]);
	double c = A[0][0] * A[1][1] - A[0][1] * A[1][0];
	cout << "det(λI-A) = λ^2 ";
	if (b != 0){
		if (b > 0){
			cout << "+";
		} else {
			cout << "-";
		}
		cout << " ";
		if (abs(b) != 1){
			cout << abs(b);
		}
		cout << "λ";
	}
	cout << " ";
	if (c != 1){
		if (c > 0){
			cout << "+";
		}
		else {
			cout << "-";
		}
		cout << " ";
		cout << abs(c);
	}
	cout << " = 0" << endl;
	cout << endl;

	//透過判斷式判斷是否有解
	if (b * b - 4 * a * c < 0){
		cout << "λ 無解" << endl;
		system("pause");
		return 0;
	}

	//透過公式解出兩根
	double r1 = (-b + sqrt(b * b - 4 * a * c)) / (2 * a);
	double r2 = (-b - sqrt(b * b - 4 * a * c)) / (2 * a);

	if (r1 == r2){
		cout << "λ = " << r1 << endl;
		sub_into(A, r1);
	} else {
		cout << "λ = " << r1 << " or " << r2 << endl;
		sub_into(A, r1);
		sub_into(A, r2);
	}
	system("pause");
	return 0;
}

//代數入字
void sub_into(vector< vector<double>> A, double r){
	cout << "eigenvector corresponding to λ = " << r << endl;
	vector< vector<double>> B = {
		{ r - A[0][0], 0 - A[0][1] },
		{ 0 - A[1][0], r - A[1][1] }
	};
	cout << " " << setw(4) << B[0][0] << ", " << setw(4) << B[0][1] << " *  x1 = 0" << endl;
	cout << " " << setw(4) << B[1][0] << ", " << setw(4) << B[1][1] << "    x2   0" << endl << endl;
	cout << "Gaussian => " << endl;
	gaussian_elimination(&B, 2);
	cout << endl;
	cout << " " << setw(4) << B[0][0] << ", " << setw(4) << B[0][1] << endl;
	cout << " " << setw(4) << B[1][0] << ", " << setw(4) << B[1][1] << endl;
	if (B[0][0] == 0){
		cout << "eigenvector = t" << endl;
		if (B[0][1] == 0){
			cout << "              s" << endl;
		} else {
			cout << "              0" << endl;
		}
	}
	else {
		if (B[1][1] == 1){
			cout << "eigenvector = 0" << endl;
			cout << "              0" << endl;
		} else if (B[0][1] == 0){
			cout << "eigenvector = 0" << endl;
			cout << "              t" << endl;
		} else {
			cout << "eigenvector = " << 0 - B[0][1] << "s" << endl;
			cout << "              " << B[0][0] << "s" << endl;
		}
	}
	cout << endl << endl;
}

//高斯消去法
void gaussian_elimination(vector< vector<double>> *matrix, int m)
{
	for (int i = 0; i < m; ++i)
	{
		// 如果上方row的首項係數為零，則考慮與下方row交換。
		if ((*matrix)[i][i] == 0){
			for (int j = i + 1; j < m; ++j){
				if ((*matrix)[j][i] != 0)
				{
					// 交換上方row與下方row。
					for (int k = i; k < m; ++k)
						swap((*matrix)[i][k], (*matrix)[j][k]);
					break;
				}
			}
		}

		if ((*matrix)[i][i] == 0){
			continue;
		}

		// 上方row的首項係數調整成一。目的是為了讓對角線皆為一。
		double t = (*matrix)[i][i];    // 首項係數
		for (int k = i; k < m; ++k){
			(*matrix)[i][k] /= t;
		}

		// 窮舉並消去下方row。
		for (int j = i + 1; j < m; ++j){
			if ((*matrix)[j][i] != 0)
			{
				double t = (*matrix)[j][i];
				for (int k = i; k < m; ++k){
					(*matrix)[j][k] -= (*matrix)[i][k] * t;
				}
			}
		}
	}
	
	//排序
	int k = 0;
	for (int j = 0; j < m; j++){
		for (int i = 0; i < m; ++i){
			if ((*matrix)[i][j] != 0){
				// 交換上方row與下方row。
				for (int n = i; n < m; ++n){
					swap((*matrix)[i][n], (*matrix)[k][n]);
				}
				break;
			}
		}
	}
}