雅克比迭代法求解线性方程组解

作者在 2012-02-16 21:28:50 发布以下内容
#include <iostream>
#include <math.h>
#include <stdlib.h>
#include <conio.h>  //包含


#define KMAX 1000  //给定最大迭代次数
#define EPS 0.000001  //给定精度要求
using namespace std;

static int n=3;//方程组阶数
static double aa[3][3]= {{8,-3,2},{4,11,-1},{6,3,12}}; //方程组系数矩阵
static double bb[3]= {20,33,36}; //非齐次线性方程组增广向量


double max_eps(double x[],double y[])
{
    int i;
    double eps=0.0;
    for(i=1; i<=n; i++)
    {
        if(fabs(y[i]-x[i])>eps)
            eps=fabs(y[i]-x[i]);
    }
    return eps;
}
int main()
{
    int k;//迭代次数
    int i,j;//控制迭代公式的实现
    double eps;//误差变量
    double a[n+1][n+1],b[n+1],x[n+1],y[n+1];

    for(i=1; i<=n; i++)
    {
        for(j=1; j<=n; j++)
        {
            a[i][j]=aa[i-1][j-1];
            if(j==i)
            {
                if(a[i][j]==0)
                {
                    cout<<"系数矩阵奇异,不可用雅克比迭代法求解"<<endl;
                    system("pause");
                    exit(-1);
                }
            }
        }
        b[i]=bb[i-1];
    }//把系数矩阵给a,把增广向量给b;且每组0元素空出,便于一下操作
    x[n+1]= {0}; //赋值0给迭代向量
    k=1;
lab_1:
    for(i=1; i<=n; i++)
    {
        y[i]=0;
        for(j=0; j<=n; j++)
        {
            if(i!=j)
                y[i]=y[i]-a[i][j]*x[j];
        }
        y[i]=(y[i]+b[i])/a[i][i];
    }
    k++;
    if(max_eps(x,y)<=EPS)
    {
        cout<<"经过"<<k<<"次雅克比迭代后得";
        cout<<"方程组解:"<<endl;
        for(int i=1; i<=n; i++)
        {
            cout<<"x"<<i<<"="<<y[i]<<"; ";
        }
        cout<<endl;
        return 0;
    }

    if(k>KMAX)
    {
        cout<<"在有效迭代次数里未找到满足要求解"<<endl;
        exit(-1);
    }
    else
    {
        for(int i=0; i<=n; i++)
        {
            x[i]=y[i];
        }
        goto lab_1;
    }

}
文章评论,共0条
游客请输入验证码
浏览66682次