怎样用高斯-赛德尔迭代法求解矩阵方程组

发布网友 发布时间：2022-05-11 21:25

我来回答

共1个回答

热心网友 时间：2023-10-21 03:21

// Seidel.h: interface for the CSeidel class.
//
//////////////////////////////////////////////////////////////////////

#if !defined(AFX_SEIDEL_H__35754D65_C3B8_4814_B9D7_8DE3BA72EFF3__INCLUDED_)
#define AFX_SEIDEL_H__35754D65_C3B8_4814_B9D7_8DE3BA72EFF3__INCLUDED_

#if _MSC_VER > 1000
#pragma once
#endif // _MSC_VER > 1000

//Seidel算法
//方法取自《计算方法引论》(第二版)徐萃薇、孙绳武著 高等教育出版社 第187页
//Matrix--系数矩阵，Y--常数项，X0--初始值，dimension--方程的阶数，error--误差；
//count--计算次数，达到次数即使不满足精度也返回积分值
//计算结果在X0中。
class CSeidel
{
public:
static bool Seidel(double *Matrix, double *Y, double *X0, int dimension, double error, int count);
CSeidel();
virtual ~CSeidel();

};

#endif // !defined(AFX_SEIDEL_H__35754D65_C3B8_4814_B9D7_8DE3BA72EFF3__INCLUDED_)

// Seidel.cpp: implementation of the CSeidel class.
//
//////////////////////////////////////////////////////////////////////

#include "stdafx.h"
//#include "NumericalMethods.h"
#include "Seidel.h"

#ifdef _DEBUG
#undef THIS_FILE
static char THIS_FILE[]=__FILE__;
#define new DEBUG_NEW
#endif

//////////////////////////////////////////////////////////////////////
// Construction/Destruction
//////////////////////////////////////////////////////////////////////

CSeidel::CSeidel()
{

}

CSeidel::~CSeidel()
{

}
//Seidel算法
//方法取自《计算方法引论》(第二版)徐萃薇、孙绳武著 高等教育出版社 第187页
//Matrix--系数矩阵，Y--常数项，X0--初始值，dimension--方程的阶数，error--误差；
//count--计算次数，达到次数即使不满足精度也返回积分值
//计算结果在X0中。
#define Matrix(row,col) (*(Matrix+(row)*dimension+col))

bool CSeidel::Seidel(double *Matrix, double *Y, double *X0, int dimension, double error, int count)
{
int i,j,k=1;
double *X=new double[dimension];
double sum;

do
{
sum=0.0f;
for(i=1;i<dimension;i++)
sum+=Matrix(0,i)*X0[i];
X[0]=(Y[0]-sum)/Matrix(0,0);
for(i=1;i<dimension-1;i++)
{
sum=0.0f;
for(j=0;j<dimension;j++)
sum+=(j==i?0:(j<i?Matrix(i,j)*X[j]:Matrix(i,j)*X0[j]));
X[i]=(Y[i]-sum)/Matrix(i,i);
}
sum=0.0f;
for(i=0;i<dimension-1;i++)
sum+=Matrix(dimension-1,i)*X[i];
X[dimension-1]=(Y[dimension-1]-sum)/Matrix(dimension-1,dimension-1);
sum=0.0f;
for(i=0;i<dimension;i++)//计算误差的范数，这里使用欧氏范数
sum+=(X[i]-X0[i])*(X[i]-X0[i]);
for(i=0;i<dimension;i++)
X0[i]=X[i];
if(sum<error*error)
{
delete[] X;
return true;
}
if(k++>count)
{
delete[] X;
return false;
}
}while(true);
}