c++ 反矩陣解矩陣方程式

#include <math.h>

#include <stdio.h>

#include <vector>

using namespace std;

using F_matrix = vector<vector<double>>;// matrix store raw vectors

void vectorShow(vector<double> v) {

for (int i = 0; i < v.size(); i ++) { printf("%10.4f", v[i]); }

printf("\n");

}

void matrixShow(F_matrix matrix) {

for (int j = 0; j < matrix.size(); j ++) { vectorShow(matrix[j]); }

printf("\n");

}

F_matrix arrayClone(double *array, int w, int h) {// copy array elements one by one

F_matrix temp(h);

for (int i = 0; i < h; i++) {

temp[i].resize(w);

for (int n = 0; n < w; n++) { temp[i][n]= *array ++; }

}

return temp;

}

F_matrix subMatrix(F_matrix &mat, int size, int xcol, int xrow = 0) {// remove xrow and xcol

int nn = mat.size();

if (size <= 0 || size > nn) return F_matrix(0);

F_matrix temp(size, vector<double>(size));

for (int j = 0, xj = 0; j < nn; j ++) {// row scan, skip 1st row

if (j == xrow) continue; // skip row xrow

for (int i = 0, xi = 0; i < nn; i++) { // column scan

if (i == xcol) continue;// skip column xcol

temp[xj][xi ++] = mat[j][i];// copy element one by one,

}

xj ++;

}

return temp;

}

double determinant(F_matrix &x, int s) {

if (s == 0) return 0.0;

int n = x.size();

if (n != x[0].size()) return 0;

if (n == 1) return x[0][0];

if (n == 2) return x[0][0] * x[1][1] - x[0][1] * x[1][0];

int det = 0;

bool positive = true;

for (int k = 0; k < n; k++) {

F_matrix subXk = subMatrix(x, n - 1, k);// remove column k, row 0

double temp = x[0][k] * determinant(subXk, n - 1);// recursive

det += positive ? temp : -temp;

positive = ! positive; // toggle the sign

}

return det;

}

double determinant(F_matrix &x) { return determinant(x, x.size()); }

F_matrix inverse(F_matrix &x) {// cramer's rule [x]⁻ = adj([x])ᵗ/det([x])

double det = determinant(x);

if (fabs(det) < 1e-16) return F_matrix(0);

int xdim = x.size();

F_matrix inv(xdim, vector<double>(xdim));

int k = xdim - 1; // to find cofactor

for (int j = 0; j < xdim; j ++) {

for (int i = 0; i < xdim; i ++) {

F_matrix xijSub = subMatrix(x, k, i, j);// remove col i and row j

double cji = determinant(xijSub); // cofactor of x[j][i]

if ((j + i) % 2) cji = - cji; // sign if odd

inv[i][j] = cji / det;// transpose and divdide determinant at the same time

}

}

return inv;

}

vector<double> xSolve(F_matrix &A, vector<double> b) { // x̂ᵗ = A⁻¹b̂ᵗ

int n = A.size();

vector<double> ans(A.size());

if (n == b.size()) {

F_matrix invA = inverse(A);

if (invA.size() > 0) {

for(int j = 0; j < n; j++) {

double temp = 0;

for(int i=0; i < n; i++) { // Aj⁻¹ • b̂ᵗ

temp += invA[j][i] * b[i];

}

ans[j] = temp;

}

}

}

return ans;

}

int main(int argc, char **argv) {

double a[] = { // Axᵗ = bᵗ, xᵗ = [w,y,z,u]ᵗ, bᵗ = [1,3,4,5]ᵗ, A:

1 , 2, 3, 4, // 1w + 2y + 3z + 4u = 1

4 , 5, 6, 2, // 4w + 5y + 6z + 2u = 3

17, 8, 9, 3, // 17w + 8y + 9z + 3u = 4

14, 7, 5, 6 // 14w + 7y + 5z + 6u = 5

};

vector<double> b{1, 3, 4, 5};

F_matrix A = arrayClone(a, 4, 4);

printf("\nSolve equation: Ax̂ᵗ = b̂ᵗ\n");

matrixShow(A);

vectorShow(b);

printf("\nAnswer is x̂ᵗ = A⁻¹b̂ᵗ:\n");

vector<double> x = xSolve(A, b);

vectorShow(x);

}