#include <bits/stdc++.h>
using namespace std;
vector<double> mul(vector<vector<double>> A,vector<double> B) //A:M*N B:N*1
{
int m=A.size(),n=B.size();
vector<double>C(m);
for(int i=0;i<m;++i){
for(int j=0;j<n;++j){
C[i]+=A[i][j]*B[j];
}
}
return C;
}
vector<double> add(vector<double> A,vector<double> B)
{
int n=A.size();
for(int i=0;i<n;++i) A[i]+=B[i];
return A;
}
double fanshu(vector<double> A)
{
double big=0;
for(double x:A)
{
big=max(big,fabs(x));
}
return big;
}
vector<double> sub(vector<double> A,vector<double> B)
{
int n=A.size();
for(int i=0;i<n;++i) A[i]-=B[i];
return A;
}
vector<double>Jacobi(vector<vector<double>>A,vector<double>y,vector<double>X_0,int& times,double eps){ //A系数矩阵,X_0初始迭代解
int n=y.size();
vector<vector<double>>B;
vector<double>g(n);
for(int i=0;i<n;++i){
vector<double>B_row(n);
for(int j=0;j<n;++j){
if(i!=j)
B_row[j]=-A[i][j]/A[i][i];
}
B.emplace_back(B_row);
g[i]=y[i]/A[i][i];
}
int k=1;
vector<double>X_k= add(mul(B,X_0),g);
while(fanshu(sub(X_0,X_k))>eps && k<times)
{
k++;
X_0=X_k;
X_k = add(mul(B,X_0),g);
}
times=k;
return X_k;
}
vector<double> Gauss_Seidel(vector<vector<double>>A,vector<double>y,vector<double>X_0,int& times,double eps)
{
int n=y.size();
vector<vector<double>>B;
vector<double>g(n);
for(int i=0;i<n;++i){
vector<double>B_row(n);
for(int j=0;j<n;++j) //求B系数矩阵
{
if(i!=j)
B_row[j]=-A[i][j]/A[i][i];
}
B.emplace_back(B_row);
g[i]=y[i]/A[i][i];
}
vector<double>X_k(n);
int k;
for(k=1;k<times;++k){ //times 次迭代计算
fill(X_k.begin(),X_k.end(),0); //X_k重新赋值计算
for(int i=0;i<n;++i) //计算X_k[i]
{
for(int j=0;j<n;++j)
{
if(j<i) X_k[i] += X_k[j]*B[i][j]; //结合i之前刚算出的X_k[j]计算,加快收敛速度
else if(j>i) X_k[i] += X_0[j]*B[i][j]; //i之后的用上次算的X_0[j]
}
X_k[i] += g[i];
}
if(fanshu(sub(X_0,X_k))<eps) break;
X_0=X_k;
}
times=k;
return X_k;
}
vector<double> Relaxtion(vector<vector<double>>A,vector<double>y,vector<double>X_0,int& times,double eps,double w){
int n=y.size();
vector<vector<double>>B;
vector<double>g(n);
for(int i=0;i<n;++i){
vector<double>B_row(n);
for(int j=0;j<n;++j) //求B系数矩阵
{
if(i!=j)
B_row[j]=-A[i][j]/A[i][i];
}
B.emplace_back(B_row);
g[i]=y[i]/A[i][i];
}
vector<double>X_k(n);
int k;
for(k=1;k<times;++k){ //times 次迭代计算
fill(X_k.begin(),X_k.end(),0); //X_k重新赋值计算
for(int i=0;i<n;++i) //计算X_k[i]
{
for(int j=0;j<n;++j)
{
if(j<i) X_k[i] += X_k[j]*B[i][j]; //结合i之前刚算出的X_k[j]计算,加快收敛速度
else if(j>i) X_k[i] += X_0[j]*B[i][j]; //i之后的用上次算的X_0[j]
}
X_k[i] = (g[i]+X_k[i])*w + (1-w)*X_0[i]; //松弛迭代主要公式
}
if(fanshu(sub(X_0,X_k))<eps) break;
X_0=X_k;
}
times=k;
return X_k;
}
signed main()
{
vector<vector<double>> A={{2,-1,-1},{1,5,-1},{1,1,10}};
vector<double> g={-5,8,11};
int times=10000;
vector<double>ans=Relaxtion(A,g,{0,0,0},times,0.001,0.5);
for(int i=0;i<ans.size();++i){
cout<<ans[i]<<" ";
}
cout<<endl;
cout<<times<<endl;
system("pause");
return 0;
}
线性方程组迭代
于 2022-07-09 22:10:10 首次发布