进击的巨人番外篇3:雅可比(Jacobi)迭代算法的C++实现 - 楼竞网站

来源:百度文库 编辑:九乡新闻网 时间:2024/05/02 19:10:29
楼竞网站君子以自强不息,君子以厚德载物。
  • 博客首页
  • 亲亲宝贝
  • 影像乐园
  • 数码世界
  • 程序设计
  • 学术科研
  • 朝花夕拾
  • 给我留言
上一篇 | 下一篇 程序设计

雅可比(Jacobi)迭代算法的C++实现

字体大小: 小 中 大 /*---------------------------------------------
假设有如下方程组:
Ax=b
用Jacobi迭代法求解方程组的解

方法:将A分裂为A=D-L-U,等价的迭代方程组为x=Bx+f。

有关算法的详细说明,请点击此地址查看。
---------------------------------------------*/

复制内容到剪贴板 程序代码
#include
#include
#include

/*楼竞网站www.LouJing.com
拥有该程序的版权,转载请保留该版权.
谢谢合作!*/
double* allocMem(int ); //分配内存空间函数
void GaussLineMain(double*,double*,double*,int );//采用高斯列主元素消去法求解x的初始向量值
void Jacobi(double*,double*,double*,double*,int,int);//利用雅可比迭代公式求解x的值

void main()
{
    short matrixNum; //矩阵的行数(列数)
    double *matrixA; //矩阵A,初始系数矩阵
    double *matrixD; //矩阵D为A中的主对角阵
    double *matrixL; //矩阵L为A中的下三角阵
    double *matrixU; //矩阵U为A中的上三角阵
    double *B;       //矩阵B为雅可比方法迭代矩阵
    double *f;       //矩阵f为中间的过渡的矩阵
    double *x;       //x为一维数组,存放结果
    double *xk;      //xk为一维数组,用来在迭代中使用
    double *b;       //b为一维数组,存放方程组右边系数

    int i,j,k;

    cout<<"<<请输入矩阵的行数(列数与行数一致)>>:";
    cin>>matrixNum;

    //分别为A、D、L、U、B、f、x、b分配内存空间
    matrixA=allocMem(matrixNum*matrixNum);
    matrixD=allocMem(matrixNum*matrixNum);
    matrixL=allocMem(matrixNum*matrixNum);
    matrixU=allocMem(matrixNum*matrixNum);
    B=allocMem(matrixNum*matrixNum);
    f=allocMem(matrixNum);
    x=allocMem(matrixNum);
    xk=allocMem(matrixNum);
    b=allocMem(matrixNum);

    //输入系数矩阵各元素值
    cout<>:"<    for(i=0;i    {
        cout<<"请输入矩阵中第 "<        for(j=0;j            cin>>*(matrixA+i*matrixNum+j);
    }

    //输入方程组右边系数b的各元素值
    cout<>:"<    for(i=0;i        cin>>*(b+i);

    /*  下面将A分裂为A=D-L-U */
    //首先将D、L、U做初始化工作
    for(i=0;i        for(j=0;j            *(matrixD+i*matrixNum+j)=*(matrixL+i*matrixNum+j)=*(matrixU+i*matrixNum+j)=0;
    //D、L、U分别得到A的主对角线、下三角和上三角;其中D取逆矩阵、L和U各元素取相反数
    for(i=0;i        for(j=0;j            if(i==j&&*(matrixA+i*matrixNum+j)) *(matrixD+i*matrixNum+j)=1/(*(matrixA+i*matrixNum+j));
            else if(i>j) *(matrixL+i*matrixNum+j)=-*(matrixA+i*matrixNum+j);
            else *(matrixU+i*matrixNum+j)=-*(matrixA+i*matrixNum+j);
    //求B矩阵中的元素
    for(i=0;i        for(j=0;j        {
            double temp=0;
            for(k=0;k                temp+=*(matrixD+i*matrixNum+k)*(*(matrixL+k*matrixNum+j)+*(matrixU+k*matrixNum+j));
            *(B+i*matrixNum+j)=temp;
        }
    //求f中的元素
    for(i=0;i    {
        double temp=0;
        for(j=0;j            temp+=*(matrixD+i*matrixNum+j)*(*(b+j));
        *(f+i)=temp;
    }

    /*  计算x的初始向量值 */
    GaussLineMain(matrixA,x,b,matrixNum);

    /* 利用雅可比迭代公式求解xk的值 */
    int JacobiTime;
    cout<>:";
    cin>>JacobiTime;
    while(JacobiTime<=0)
    {
        cout<<"迭代次数必须大于0,请重新输入:";
        cin>>JacobiTime;
    }
    Jacobi(x,xk,B,f,matrixNum,JacobiTime);

    //输出线性方程组的解 */
    cout<>"<    cout.precision(20); //设置输出精度,以此比较不同迭代次数的结果
    for(i=0;i        cout<<"x"<
    cout<
    //释放掉所有动态分配的内存
    delete [] matrixA;
    delete [] matrixD;
    delete [] matrixL;
    delete [] matrixU;
    delete [] B;
    delete [] f;
    delete [] x;
    delete [] xk;
    delete [] b;
}


/*--------------------
  分配内存空间函数
  www.LouJing.com
--------------------*/
double* allocMem(int num)
{
    double *head;
    if((head=new double[num])==NULL)
    {
        cout<<"内存空间分配失败,程序终止运行!"<        exit(0);
    }
    return head;
}


/*---------------------------------------------
  计算Ax=b中x的初始向量值,采用高斯列主元素消去法
  www.LouJing.com
---------------------------------------------*/
void GaussLineMain(double* A,double* x,double* b,int num)
{
    int i,j,k;

    //共处理num-1行
    for(i=0;i    {
        //首先每列选主元,即最大的一个
        double lineMax=fabs(*(A+i*num+i));
        int lineI=i;
        for(j=i;j            if(fabs(*(A+j*num+i))>fabs(lineMax)) lineI=j;

        //主元所在行和当前处理行做行交换,右系数b也随之交换
        for(j=i;j        {
            //A做交换
            lineMax=*(A+i*num+j);
            *(A+i*num+j)=*(A+lineI*num+j);
            *(A+lineI*num+j)=lineMax;
            //b中对应元素做交换
            lineMax=*(b+i);
            *(b+i)=*(b+lineI);
            *(b+lineI)=lineMax;
        }

        if(*(A+i*num+i)==0) continue; //如果当前主元为0,本次循环结束

        //将A变为上三角矩阵,同样b也随之变换
        for(j=i+1;j        {
            double temp=-*(A+j*num+i)/(*(A+i*num+i));
            for(k=i;k            {
                *(A+j*num+k)+=temp*(*(A+i*num+k));
            }
            *(b+j)+=temp*(*(b+i));
        }
    }

    /* 验证Ax=b是否有唯一解,就是验证A的行列式是否为0;
    如果|A|!=0,说明有唯一解*/
    double determinantA=1;
    for(i=0;i        determinantA*=*(A+i*num+i);
    if(determinantA==0)
    {
        cout<        exit(0);
    }

    /* 从最后一行开始,回代求解x的初始向量值 */
    for(i=num-1;i>=0;i--)
    {
        for(j=num-1;j>i;j--)
            *(b+i)-=*(A+i*num+j)*(*(x+j));
        *(x+i)=*(b+i)/(*(A+i*num+i));
    }
}


/*------------------------------------
利用雅可比迭代公式求解x的精确值
www.LouJing.com
迭代公式为:xk=Bx+f
------------------------------------*/
void Jacobi(double* x,double* xk,double* B,double* f,int num,int time)
{
    int t=1,i,j;
    while(t<=time)
    {
        for(i=0;i        {
            double temp=0;
            for(j=0;j                temp+=*(B+i*num+j)*(*(x+j));
            *(xk+i)=temp+*(f+i);
        }

        //将xk赋值给x,准备下一次迭代
        for(i=0;i            *(x+i)=*(xk+i);
        t++;
    }
}



[本日志由 楼竞 于 2010-02-15 12:19 AM 编辑]文章来自: 本站原创
引用通告: 查看所有引用 | 我要引用此文章
Tags: C++C++
相关日志:
  • 欧几里得算法[320]
  • Google Code Jam-中国编程挑战赛-资格赛-DiskClusters[225]
  • Google Code Jam-中国编程挑战赛-资格赛-BusStops[227]
  • 求一个大数的阶乘[306]
  • 单精度浮点数在VC++6.0中内存格式研究[237] 分页:[1]
  • 评论: 0 | 引用: 0 | 查看次数: 307 发表评论
    昵 称: 密 码: 游客发言不需要密码. 邮 箱: 支持Gravatar头像. 网 址: 输入网址便于回访. 内 容:
    正在加载编辑器... 验证码: 选 项:
    虽然发表评论不用注册,但是为了保护您的发言权,建议您注册帐号.
    字数限制 1000 字 | UBB代码 开启 | [img]标签 关闭

    topnav


    • 网站首页

    • 个人档案

    • 热门标签

    • 内容订阅

    Category

    博客首页
    亲亲宝贝 [68]
    影像乐园 [4]
    数码世界 [6]
    程序设计 [31]
    学术科研 [3]
    朝花夕拾 [12]
    给我留言
    Tags Cloud

    Calendar

    2011年1月
    • 26
    • 27
    • 28
    • 29
    • 30
    • 31
    • 1
    • 2
    • 3
    • 4
    • 5
    • 6
    • 7
    • 8
    • 9
    • 10
    • 11
    • 12
    • 13
    • 14
    • 15
    • 16
    • 17
    • 18
    • 19
    • 20
    • 21
    • 22
    • 23
    • 24
    • 25
    • 26
    • 27
    • 28
    • 29
    • 30
    • 31
    • 1
    • 2
    • 3
    • 4
    • 5

    Recent Comments

    折腾这个好纠结..... 这是在闸口? 哈哈~桂林好玩的~一次两次哪够啊! 楼老师什么时...楼总忙归忙,对宝宝的关心是一点也没少嘛。师母这个身材赶上陈慧琳了。。 长得不太像你哈 屋后的不知道是麦田还是稻田,您太幽默了 楼总想多了 岳岳,你儿子实在很会享受啊!亲爱的宝宝  睡吧    你会梦见花园里有朵红玫瑰...

    User Panel

    • 登录
    • 用户注册

    Geust Book

    好崇拜老师的学术研究。牛!!楼竞兄,您好! 我想学习一下你的 C++Buil... 网站不错... 加个朋友吧... jp...小孩脸肥嘟嘟的蛮好玩的,什么时候老师你小时候的照片...  楼老师! 支持一下~~~[b]要是宝宝照片能和楼老师小时候照片放在一起对比...楼宝宝挺象楼爸的.多放些宝妈的照片来欣赏一下吧.还是没更新啊既然来了,就要留言咯~~~ 网站办的还不错,就是...没更新啊

    Archive

    201002月200903月04月05月06月07月
    12月200803月200711月200601月02月03月04月200502月03月09月10月11月
    12月200401月07月09月12月200305月06月07月08月10月

    Support

     

    Powered By LouJing.Com CopyRight 2003 - 2011 , 楼竞网站 xhtml | css

    Processed in 0.250000 second(s) , 4 queries