数值实验共轭梯度法matlab代码,数值线性代数实验-共轭梯度法

一开始用c++的运算符重载程序总是莫名其妙的崩掉，然后以为是运算符重载的问题就写了个class对矩阵重新封装，结果还是崩，然后好久才发现是我把空间开的太大导致程序崩掉，无语，这样就浪费了我一个上午。。。。

课本上的例题：

$$

\begin{equation}

{\left[ \begin{array}{ccc}

10&1&2&3&4\\

1& 9& -1& 2& -3\\

2& -1& 7& 3& -5\\

3& 2& 3& 12& -1\\

4& -3& -5 &-1 &15

\end{array}

\right ]}\times {

\left[ \begin{array}{ccc}

x_1\\

x_2\\

x_3\\

x_4\\

x_5

\end{array}

\right ]}

={

\left[ \begin{array}{ccc}

12\\ -27\\ 14\\ -17\\ 12

\end{array}

\right ]}

\end{equation}

$$

输入数据

5 5

10 1 2 3 4

1 9 -1 2 -3

2 -1 7 3 -5

3 2 3 12 -1

4 -3 -5 -1 15

12 -27 14 -17 12

#include "cstdio"#include"cstring"#include"cstdlib"#include"cmath"#include"iostream"

using namespacestd;const double eps = 1e-6;const int maxr = 10;const int maxc = 10;classMatrix {private:doublem[maxr][maxc];introw, col;public:

Matrixoperator -(Matrix b) {

Matrix c;

c.col= this->col; c.row = this->row;for (int i = 0; i < this->row; i++) {for (int j = 0; j < this->col; j++) {

c.m[i][j]= this->m[i][j]-b.m[i][j];

}

}returnc;

}

Matrixoperator +(Matrix b) {

Matrix c;

c.col= this->col; c.row = this->row;for (int i = 0; i < this->row; i++) {for (int j = 0; j < this->col; j++) {

c.m[i][j]= this->m[i][j] +b.m[i][j];

}

}returnc;

}

Matrixoperator *(Matrix b) {

Matrix c;

memset(c.m,0, sizeof(c.m));

c.row= this->row; c.col =b.col;for (int i = 0; i < this->row; i++) {for (int k = 0; k < this->col; k++) {for (int j = 0; j < b.col; j++) {

c.m[i][j]+= this->m[i][k]*b.m[k][j];

}

}

}returnc;

}

Matrixoperator & (doublea) {

Matrix c= *this;for (int i = 0; i < this->row; i++) {for (int j = 0; j < this->col; j++) {

c.m[i][j]= a*this->m[i][j];

}

}returnc;

}

Matrixoperator !() {

Matrix c= *this;

c.row= this->col, c.col = this->row;for (int i = 0; i < this->row; i++) {for (int j = 0; j < this->col; j++) {

c.m[j][i]= this->m[i][j];

}

}returnc;

}intdcmp() {for (int i = 0; i < this->row; i++) {for (int j = 0; j < this->col; j++) {if (fabs(this->m[i][j]) > eps) return 1;

}

}return 0;

}voidshow_Matrix() {for (int i = 0; i < this->row; i++) {

printf("%.6f", this->m[i][0]);for (int j = 1; j < this->col; j++) {

printf("%.8f", this->m[i][j]);

}

printf("\n");

}

}void set_matrix(int n, intm) {this->row = n, this->col =m;for (int i = 0; i < n; i++) {for (int j = 0; j < m; j++) {

scanf("%lf", &this->m[i][j]);

}

}

}void init(int n, intm) {this->row = n, this->col =m;for (int i = 0; i < n; i++) {for (int j = 0; j < m; j++) {this->m[i][j] = 1.0;

}

}

}int Col(){return this->col;}int Row(){return this->row;}double get(int i, int j) {return this->m[i-1][j-1];}

};

Matrix Conjugate_gradient(Matrix A, Matrix B) {

Matrix x0=B, p0, r1, p1;doublealf, bet;

x0.init(A.Col(),1);

Matrix r0= B-(A*x0);int k = 0;

r1=r0;while(r1.dcmp()) {

k++;if (k == 1) p1 =r0;else{

bet= (!r1*r1).get(1, 1)/(!r0*r0).get(1, 1);

p1= r1 + (p0&bet);

}

alf= (!r1*r1).get(1, 1)/(!p1*A*p1).get(1, 1);

x0= x0 + (p1&alf);

r0=r1;

r1= r1-((A&alf)*p1);

}returnx0;

}int main(int argc, char const *argv[])

{//freopen("in.txt", "r", stdin);

Matrix A, B, C;intn, m;

printf("输入系数矩阵的维数和矩阵A\n");

scanf("%d %d", &n, &m);

A.set_matrix(n, m);

printf("输入矩阵B\n");

B.set_matrix(n,1);

C=Conjugate_gradient(A, B);

printf("利用共轭梯度法得到的解为\n");

C.show_Matrix();return 0;

}