c_cpp 将方程式解析成矩阵
Posted
tags:
篇首语：本文由小常识网(cha138.com)小编为大家整理，主要介绍了c_cpp 将方程式解析成矩阵相关的知识，希望对你有一定的参考价值。
#include <iostream>
#include <vector>
#include <string>
#include <cmath>
#include <iomanip>
#include <cstdio>
#include <regex>
using namespace std;
class Matrix {
vector< vector<double> > data;
public:
Matrix() {};
Matrix(vector< vector<double> > &_data) {
data = _data;
}
Matrix(size_t rowSize, size_t colSize) : data(rowSize, vector<double>(colSize, 0)) {}
void setElement(size_t rowNum, size_t colNum, double value) {
data[rowNum][colNum] = value;
}
void setRow(size_t rowNum, vector<double> v) {
data[rowNum] = v;
}
double getElement(size_t rowNum, size_t colNum) {
return data[rowNum][colNum];
}
vector<double> getRow(size_t rowNum) {
return data[rowNum];
}
size_t colSize() {
return data[0].size();
}
size_t rowSize() {
return data.size();
}
friend ostream& operator<<(ostream& os, Matrix& m) {
size_t rowSize = m.rowSize(), colSize = m.colSize();
for (int i = 0; i<rowSize; ++i) {
for (int j = 0; j<colSize; ++j) {
double value = m.getElement(i, j);
os << setprecision(log10(abs(value)) + 3) << value;
if (j != colSize - 1) {
os << " ";
}
}
os << endl;
}
return os;
}
friend istream& operator>>(istream& is, Matrix& m) {
size_t rowSize = m.rowSize(), colSize = m.colSize();
for (int i = 0; i<rowSize; ++i) {
for (int j = 0; j<colSize; ++j) {
double value;
is >> value;
m.setElement(i, j, value);
}
}
return is;
}
Matrix operator+(Matrix &m2) {
if (rowSize() != m2.rowSize() || colSize() != m2.colSize()) {
throw "matrix size is not match";
}
size_t rowSize = this->rowSize(), colSize = this->colSize();
Matrix result(rowSize, colSize);
for (size_t i = 0; i<rowSize; ++i) {
for (size_t j = 0; j<colSize; ++j) {
double eleValue = getElement(i, j) + m2.getElement(i, j);
result.setElement(i, j, eleValue);
}
}
return result;
}
Matrix operator-(Matrix &m2) {
if (rowSize()!=m2.rowSize() || colSize() != m2.colSize()) {
throw "matrix size is not match";
}
size_t rowSize=this->rowSize(), colSize=this->colSize();
Matrix result(rowSize, colSize);
for (size_t i=0; i<rowSize; ++i) {
for (size_t j=0; j<colSize; ++j) {
double eleValue = getElement(i, j) - m2.getElement(i, j);
result.setElement(i, j, eleValue);
}
}
return result;
}
Matrix operator*(Matrix &m2) {
if (colSize() != m2.rowSize()) {
throw "matrix size is not match";
}
size_t rowSize = this->rowSize(), colSize = m2.colSize();
Matrix result(rowSize, colSize);
for (size_t i = 0; i<rowSize; ++i) {
for (size_t j = 0; j<colSize; ++j) {
double eleValue = 0;
for (size_t k = 0; k<colSize; ++k) {
eleValue += getElement(i, k) * m2.getElement(k, j);
}
result.setElement(i, j, eleValue);
}
}
return result;
}
Matrix operator-() {
for (int i=0; i<rowSize(); ++i) {
for (int j=0; j<colSize(); ++j) {
double value = getElement(i, j);
setElement(i, j, -value);
}
}
return *this;
}
};
Matrix operator*(double c, Matrix m) {
for (int i=0; i<m.rowSize(); ++i)
for (int j=0; j<m.colSize(); ++j)
m.setElement(i, j, m.getElement(i, j)*c);
return m;
}
Matrix operator/(Matrix m, double c) {
for (int i=0; i<m.rowSize(); ++i)
for (int j=0; j<m.colSize(); ++j)
m.setElement(i, j, m.getElement(i, j)/c);
return m;
}
Matrix transpose(Matrix & m)
{
size_t rowSize = m.rowSize(), colSize = m.colSize();
Matrix result(colSize, rowSize);
for (int i = 0; i<rowSize; ++i) {
for (int j = 0; j<colSize; ++j) {
double value = m.getElement(i, j);
result.setElement(j, i, value);
}
}
return result;
}
bool Matrix_Inv(double **naturalmat, int num, double **InvMat)
{
int i, j, k = 0;
double **MatEnhanced;
MatEnhanced = (double**)malloc(num*sizeof(double*));
for (i = 0; i<num; i++)
MatEnhanced[i] = (double*)malloc(2 * num*sizeof(double));
double *temp;
temp = (double*)malloc(2 * num*sizeof(double));
double xishu = 1;
for (i = 0; i<num; i++)
{
for (j = 0; j<num; j++)
MatEnhanced[i][j] = naturalmat[i][j];
}
for (i = 0; i<num; i++)
{
for (j = num; j<2 * num; j++)
MatEnhanced[i][j] = 0;
MatEnhanced[i][i + num] = 1;
}
for (i = 0; i<num; i++)
{
if (MatEnhanced[i][i] == 0)
{
if (i == num - 1)
return false;
for (j = i; j<num; j++)
{
if (MatEnhanced[j][i] != 0)
{
k = j;
break;
}
}
temp = MatEnhanced[k];
MatEnhanced[k] = MatEnhanced[i];
MatEnhanced[i] = temp;
}
for (j = 0; j<num; j++)
{
if (j != i)
{
if (MatEnhanced[j][i] != 0)
{
xishu = MatEnhanced[j][i] / MatEnhanced[i][i];
for (k = i; k<2 * num; k++)
{
MatEnhanced[j][k] -= xishu*MatEnhanced[i][k];
}
}
}
}
xishu = MatEnhanced[i][i];
for (j = i; j<2 * num; j++)
if (xishu != 0)
MatEnhanced[i][j] /= xishu;
}
for (i = 0; i<num; i++)
{
for (j = 0; j<num; j++)
InvMat[i][j] = MatEnhanced[i][j + num];
}
free(MatEnhanced);
free(temp);
return true;
}
double regexParse(string s, string es) {
std::regex e (es);
std::smatch sm;
std::regex_search(s, sm, e);
if (sm[0] == "") {
string s1 = sm[1];
string s2 = sm[2];
string s3 = sm[3];
return 0;
}
if (sm[1] == ""){
string s1 = sm[1];
string s2 = sm[2];
string s3 = sm[3];
return 1;
}
if (sm[1] == "-") {
return -1;
}
return stod(sm[1]);
}
/* x^2 x x^0.5 1
y^2 0 0 0 0
y 0 0 0 0
y^0.5 0 0 0 0
1 0 0 0 0
*/
Matrix parse(string poly) {
string token11 = "(\\-?\\d*\\.?\\d*)x\\^2\\*?y\\^2(?=$|\\+|\\-)";
string token12 = "(\\-?\\d*\\.?\\d*)x\\*?y\\^2(?=$|\\+|\\-)";
string token13 = "(\\-?\\d*\\.?\\d*)x\\^\\(0+.5\\)\\*?y\\^2(?=$|\\+|\\-)";
string token14 = "(\\-?\\d*\\.?\\d*)y\\^2(?=$|\\+|\\-)";
string token21 = "(\\-?\\d*\\.?\\d*)x\\^2\\*?y(?=$|\\+|\\-)";
string token22 = "(\\-?\\d*\\.?\\d*)x\\*?y(?=$|\\+|\\-)";
string token23 = "(\\-?\\d*\\.?\\d*)x\\^\\(0+.5\\)\\*?y(?=$|\\+|\\-)";
string token24 = "(\\-?\\d*\\.?\\d*)x\\^2\\*?y(?=$|\\+|\\-)";
string token31 = "(\\-?\\d*\\.?\\d*)x\\^2\\*?y\\^\\(0+.5\\)(?=$|\\+|\\-)";
string token32 = "(\\-?\\d*\\.?\\d*)x\\*?y\\^\\(0+.5\\)(?=$|\\+|\\-)";
string token33 = "(\\-?\\d*\\.?\\d*)x\\^\\(0+.5\\)\\*?y\\^\\(0+.5\\)(?=$|\\+|\\-)";
string token34 = "(\\-?\\d*\\.?\\d*)x\\^2\\*?y\\^\\(0+.5\\)(?=$|\\+|\\-)";
string token41 = "(\\-?\\d*\\.?\\d*)x\\^2(?=$|\\+|\\-)";
string token42 = "(\\-?\\d*\\.?\\d*)x(?=$|\\+|\\-)";
string token43 = "(\\-?\\d*\\.?\\d*)x\\^\\(0+.5\\)(?=$|\\+|\\-)";
string token44 = "[^\\^](\\-?\\d+\\.?\\d*(?=$|\\+|\\-))";
// 搜尋出來可能僅有一個 - 號代表 -1, 空為 1
Matrix result(4, 4);
result.setElement(0, 0, regexParse(poly, token11));
result.setElement(0, 1, regexParse(poly, token12));
result.setElement(0, 2, regexParse(poly, token13));
result.setElement(0, 3, regexParse(poly, token14));
result.setElement(1, 0, regexParse(poly, token21));
result.setElement(1, 1, regexParse(poly, token22));
result.setElement(1, 2, regexParse(poly, token23));
result.setElement(1, 3, regexParse(poly, token24));
result.setElement(2, 0, regexParse(poly, token31));
result.setElement(2, 1, regexParse(poly, token32));
result.setElement(2, 2, regexParse(poly, token33));
result.setElement(2, 3, regexParse(poly, token34));
result.setElement(3, 0, regexParse(poly, token41));
result.setElement(3, 1, regexParse(poly, token42));
result.setElement(3, 2, regexParse(poly, token43));
result.setElement(3, 3, regexParse(poly, token44));
return result;
}
enum DPS{
/// DIFFERENTIAL POWER SIGN
NONE = 0,
FIRST = 1,
SECOND = 2
};
void check_err_NS(Matrix FM, Matrix ini_pnt){
/// NEGATIVE VALUE IN SQRT
bool err_sign = 0;
if(ini_pnt.getElement(0, 0) < 0)
for(int j=0; j<4; j++)
if(FM.getElement(j, 2) != 0)
err_sign = 1;
if(ini_pnt.getElement(1, 0) < 0)
for(int i=0; i<4; i++)
if(FM.getElement(2, i) != 0)
err_sign = 1;
if(err_sign)
throw "error: ini value with negative and have sqrt";
return;
}
Matrix set_DC_X(DPS sign, Matrix ini_pnt){
/// SET DIFFERENTIAL COEFFICIENT OF X
/// initial the DC_X
Matrix DC_X(4, 1);
double x = ini_pnt.getElement(0, 0);
/// set DC_X
if(sign == NONE){
DC_X.setElement(0, 0, x*x );
DC_X.setElement(1, 0, x );
if (x>=0) DC_X.setElement(2, 0, sqrt(x) );
DC_X.setElement(3, 0, 1 );
}else if(sign == FIRST){
DC_X.setElement(0, 0, 2*x );
DC_X.setElement(1, 0, 1 );
if (x>=0) DC_X.setElement(2, 0, 0.5*pow(x, -0.5));
DC_X.setElement(3, 0, 0 );
}else if(sign == SECOND){
DC_X.setElement(0, 0, 2 );
DC_X.setElement(1, 0, 0 );
if (x>=0) DC_X.setElement(2, 0, 0.25*pow(x, -1.5));
DC_X.setElement(3, 0, 0 );
}else{
throw "error: out of Second differential.";
}
return DC_X;
}
Matrix set_DC_Y(DPS sign, Matrix ini_pnt){
/// SET DIFFERENTIAL COEFFICIENT OF X
/// initial the DC_Y
Matrix DC_Y(1, 4);
double y = ini_pnt.getElement(1, 0);
/// set DC_Y
if(sign == NONE){
DC_Y.setElement(0, 0, y*y );
DC_Y.setElement(0, 1, y );
if (y>=0) DC_Y.setElement(0, 2, sqrt(y) );
DC_Y.setElement(0, 3, 1 );
}else if(sign == FIRST){
DC_Y.setElement(0, 0, 2*y );
DC_Y.setElement(0, 1, 1 );
if (y>=0) DC_Y.setElement(0, 2, 0.5*pow(y, -0.5));
DC_Y.setElement(0, 3, 0 );
}else if(sign == SECOND){
DC_Y.setElement(0, 0, 2 );
DC_Y.setElement(0, 1, 0 );
if (y>=0) DC_Y.setElement(0, 2, 0.25*pow(y, -1.5));
DC_Y.setElement(0, 3, 0 );
}else{
throw "error";
}
return DC_Y;
}
Matrix set_FODM(Matrix FM, Matrix ini_pnt){
/// SET FIRST ORDER DIFFERENTIAL MATRIX
/// filter the impossible ini value in sqrt
check_err_NS(FM, ini_pnt);
/// initial FODM
Matrix FODM(2, 1);
/// calculate df/dx
Matrix DC_X = set_DC_X(FIRST, ini_pnt);
Matrix DC_Y = set_DC_Y(NONE, ini_pnt);
Matrix value_M = DC_Y*FM*DC_X;
double value = value_M.getElement(0, 0);
FODM.setElement(0, 0, value);
/// calculate df/dy
DC_X = set_DC_X(NONE, ini_pnt);
DC_Y = set_DC_Y(FIRST, ini_pnt);
value_M = DC_Y*FM*DC_X;
value = value_M.getElement(0, 0);
FODM.setElement(1, 0, value);
return FODM;
}
Matrix set_SODM(Matrix FM, Matrix ini_pnt){
/// SET SECOND ORDER DIFFERENTIAL MATRIX(Hessian matrix)
/// filter the impossible ini value in sqrt
check_err_NS(FM, ini_pnt);
/// initial SODM
Matrix SODM(2, 2);
/// calculate df/dx^2
Matrix DC_X = set_DC_X(SECOND, ini_pnt);
Matrix DC_Y = set_DC_Y(NONE, ini_pnt);
Matrix value_M = DC_Y*FM*DC_X;
double value = value_M.getElement(0, 0);
SODM.setElement(0, 0, value);
/// calculate df/dxy
DC_X = set_DC_X(FIRST, ini_pnt);
DC_Y = set_DC_Y(FIRST, ini_pnt);
value_M = DC_Y*FM*DC_X;
value = value_M.getElement(0, 0);
SODM.setElement(0, 1, value);
SODM.setElement(1, 0, value);
/// calculate df/dy^2
DC_X = set_DC_X(NONE, ini_pnt);
DC_Y = set_DC_Y(SECOND, ini_pnt);
value_M = DC_Y*FM*DC_X;
value = value_M.getElement(0, 0);
SODM.setElement(1, 1, value);
return SODM;
}
double get_fuction_value(Matrix FM, Matrix ini_pnt){
/// GET VALUE OF FUNCTION AT GIVEN POINT
/// filter the impossible ini value in sqrt
check_err_NS(FM, ini_pnt);
/// calculate f(x)
Matrix DC_X = set_DC_X(NONE, ini_pnt);
Matrix DC_Y = set_DC_Y(NONE, ini_pnt);
Matrix value_M = DC_Y*FM*DC_X;
double value = value_M.getElement(0, 0);
return value;
}
double get_fuction_value(Matrix FM, Matrix ini_pnt, Matrix S, double test_alpha){
Matrix tmp = test_alpha*S;
Matrix test_point = ini_pnt + tmp;
return get_fuction_value(FM, test_point);
}
double get_alpha(Matrix &FM, Matrix &ini_pnt, Matrix &S, double a, double b, double c) {
const double tau = 0.001;
const static double phi = (1 + sqrt(5) / 2);
const static double resphi = 2 - phi;
double x;
if (c-b > b-a)
x = b + resphi * (c-b);
else {
x = b - resphi * (b-a);
}
if (abs(c-a) < tau * abs(b) + abs(x))
return (c+a)/2;
double f_x = get_fuction_value(FM, ini_pnt, S, x);
double f_b = get_fuction_value(FM, ini_pnt, S, b);
if (f_x < f_b) {
if (c-b > b-a)
return get_alpha(FM, ini_pnt, S, b, x, c);
else {
return get_alpha(FM, ini_pnt, S, a, x, b);
}
} else {
if (c-b > b-a)
return get_alpha(FM, ini_pnt, S, a, b, x);
else
return get_alpha(FM, ini_pnt, S, x, b, c);
}
}
//double goldenSectionSearch(double a, double b, double c, double tau) {
// const static double phi = (1 + sqrt(5) / 2);
// const static double resphi = 2 - phi;
// double x;
// if (c-b > b-a)
// x = b + resphi * (c-b);
// else {
// x = b - resphi * (b-a);
// }
//
// if (abs(c-a) < tau * abs(b) + abs(x))
// return (c+a)/2;
// if (f(x) < f(b)) {
// if (c-b > b-a)
// return goldenSectionSearch(b, x, c, tau);
// else {
// return goldenSectionSearch(a, x, b, tau);
// }
// } else {
// if (c-b > b-a)
// return goldenSectionSearch(a, b, x, tau);
// else
// return goldenSectionSearch(x, b, c, tau);
// }
//}
//void Conjugate_Gradient(Matrix FM, Matrix x_now, size_t n, ){
//
// /// step 1, initialize the variable and value
// const size_t N = 1000; /// initial limit_of_literations_times
// double f_X = get_fuction_value(FM, x_now);/// initial value of f(X)
// double stop_1 = 0.1; /// initial the stop sign
// double stop_2 = 0.1;
// double stop_3 = 0.1;
// vector<Matrix> s;
//
//
// /// step 2,
// for(size_t i=0; i<N; i++){
//
// if(i==1){
// Matrix s_now = -set_FODM(FM, x_now); /// initial the FODE Matrix
// s.push_back(s_now);
// } else {
// Matrix Beta = get_Beta(X_pre, X_now);
// double S_now = -set_FODM(FM, x_now).getElement(0, 0) + Beta*S_pre;
// }
//
// Matrix alpha = get_alpha(double value);
// X_next = X_now + alpha*s[i];
// = X_next;
//
// }
//}
Matrix inverse(Matrix m) {
size_t colSize = m.colSize();
double **ma = new double*[colSize];
double **mb = new double*[colSize];
for (int i=0; i<colSize; ++i) {
ma[i] = new double[colSize];
mb[i] = new double[colSize];
for (int j=0; j<colSize; ++j)
ma[i][j] = m.getElement(i, j);
}
if (Matrix_Inv(ma, colSize, mb)) {
for (int i=0; i<colSize; ++i) {
for (int j=0; j<colSize; ++j) {
m.setElement(i, j, mb[i][j]);
}
}
return m;
}
throw "error";
}
Matrix upperTriangle(Matrix & m)
{
size_t colSize = m.colSize(), rowSize = m.rowSize();
// for column
for (int i = 0; i<rowSize; ++i) {
int notZeroRow = -1;
for (int k = i; k<rowSize; ++k) {
double value = m.getElement(k, i);
if (value != 0) {
notZeroRow = k;
break;
}
}
if (notZeroRow == -1) {
continue;
}
else {
// 交換
for (int j = 0; j<colSize; ++j) {
double v1 = m.getElement(notZeroRow, j);
double v2 = m.getElement(i, j);
m.setElement(i, j, v1);
m.setElement(notZeroRow, j, v2);
}
}
for (int j = i + 1; j<rowSize; ++j) {
double c;
if (m.getElement(notZeroRow, i) != 0) {
c = m.getElement(j, i) / m.getElement(notZeroRow, i);
}
else {
c = 0;
}
for (int k = 0; k<colSize; ++k) {
double value = m.getElement(j, k) - m.getElement(notZeroRow, k)*c;
m.setElement(j, k, value);
}
}
}
return m;
}
double determinant(Matrix & m)
{
if (m.rowSize() != m.colSize()) {
throw "m is not square matrix";
}
size_t colSize = m.rowSize();
if (colSize == 1) {
return m.getElement(0, 0);
}
if (colSize == 2) {
return m.getElement(0, 0) * m.getElement(1, 1) - m.getElement(0, 1) * m.getElement(1, 0);
}
double result = 1;
Matrix r = upperTriangle(m);
for (int i = 0; i<colSize; ++i) {
result *= r.getElement(i, i);
}
return result;
}
Matrix newton(Matrix FM, Matrix ini_pnt){
/// First part: initialize
Matrix POINT_NOW = ini_pnt;
double stop_range = 0.0001;
int counter = 0;
/// Second part: newton way to calculate the point
while(true){
double function_value = get_fuction_value(FM, POINT_NOW);
Matrix FODM = set_FODM(FM, POINT_NOW);
Matrix SODM = set_SODM(FM, POINT_NOW);
if(counter > 100)
throw "Too many times in newton.";
if(abs(FODM.getElement(0, 0)) < stop_range && abs(FODM.getElement(1, 0)) < stop_range && determinant(SODM) > 0)
break;
// if(determinant(SODM) <= 0)
// throw "Hessian is less than 0";
Matrix tmp = inverse(SODM)* FODM;
POINT_NOW = POINT_NOW - tmp;
counter++;
std::cout << "Hessian:\n" << SODM << std::endl;
tmp = inverse(SODM);
std::cout << "Hessian^(-1):\n" << tmp << std::endl << std::endl;
std::cout << "f\':\n" << FODM << std::endl;
std::cout << "Delta:\n" << tmp << std::endl;
std::cout << "Point:\n" << POINT_NOW << std::endl;
}
// std::cout << "Hessian:\n" << SODM << std::endl;
// std::cout << "Hessian^(-1):\n" << inverse(SODM) << std::endl << std::endl;
// std::cout << "f\':\n" << FODM << std::endl;
// std::cout << "Point:\n" << POINT_NOW << std::endl;
/// Third part: return the point
return POINT_NOW;
}
Matrix quasi_newton(Matrix FM, Matrix ini_pnt, double up_bound, double low_bound){
/// step 1: initialize
Matrix POINT_NOW = ini_pnt;
double stop_range = 0.0001;
Matrix F(2, 2);
F.setElement(0, 0, 1);
F.setElement(0, 1, 0);
F.setElement(1, 0, 0);
F.setElement(1, 1, 1);
Matrix FODM_NOW = set_FODM(FM, POINT_NOW);
// std::cout << "F:\n" << F << std::endl;
// std::cout << "f\':\n" << FODM_NOW << std::endl;
/// step 2:
while(abs(FODM_NOW.getElement(0, 0)) > stop_range && abs(FODM_NOW.getElement(1, 0)) > stop_range){
/// calculate FODM, POINT, their DELTA
Matrix d = -F*FODM_NOW;
double alpha = get_alpha(FM, POINT_NOW, d, up_bound, (up_bound+low_bound)/2, low_bound);
Matrix DELTA_X = alpha*d;
POINT_NOW = POINT_NOW + DELTA_X;
Matrix FODM_NEXT = set_FODM(FM, POINT_NOW);
Matrix DELTA_FODM = FODM_NEXT - FODM_NOW;
FODM_NOW = FODM_NEXT;
/// calculate new F
Matrix tmp1 = transpose(DELTA_X)*DELTA_FODM;
Matrix tmp2 = transpose(DELTA_FODM)*F*DELTA_FODM;
Matrix tmp3 = F*DELTA_FODM;
Matrix t = transpose(tmp3);
Matrix tmp4 = tmp3* t;
Matrix t2 =transpose(DELTA_X);
Matrix ty = DELTA_X*t2/tmp1.getElement(0, 0);
F = F + ty;
ty = tmp4/tmp2.getElement(0, 0);
F = F - ty;
/// output
std::cout << "F:\n" << F << std::endl;
std::cout << "f\':\n" << FODM_NOW << std::endl;
std::cout << "alpha: " << alpha << std::endl;
std::cout << "point\n" << POINT_NOW << std::endl;
}
return POINT_NOW;
}
Matrix pointToMatrix(double x, double y) {
Matrix p(2, 1);
p.setElement(0, 0, x);
p.setElement(1, 0, y);
return p;
}
int main() {
// string s = "-x+13xy+1.2";
// Matrix t = parse(s);
//
// cout << t;
newton(parse("x^2-2xy+4y^2"), pointToMatrix(-3, 1));
}以上是关于c_cpp 将方程式解析成矩阵的主要内容，如果未能解决你的问题，请参考以下文章