熟悉雅可比法的帮忙了!!!!

bobnie

雅可比法求实对称矩阵所有特征值和对应特征向量
程序如下
#define EPSL 1.0e-6                                                            //迭代的精度水平
#define PI 3.1415926


#include <stdio.h>
#include <math.h>


void main()                                                                      //主程序开始；
{
        int size,//协方差矩阵的阶数
                row,col,k,
                p,q;//记录maxa的行与列
        double a[99][99],//存放协方差矩阵
                R[99][99]={0},//平面旋转变换矩阵
                RT[99][99],//矩阵R的转置
                S[99][99],//存放特征向量，赋初值为单位矩阵
                tan,//旋转角度的正切
                alfa,//旋转角度
                maxa;//非主对角线上绝对值最大元素
        FILE *fp;

        fp=fopen("sj.txt","r");
        fscanf(fp,"%d",&size);//读入协方差矩阵的阶数

        for(row=0;row<size;row++)
                for(col=0;col<size;col++)
                        fscanf(fp,"%lf",&a[row][col]);

        fclose(fp);

        for(row=0;row<size;row++)//给矩阵S赋初值为单位矩阵(经验证无误)
                for(col=0;col<size;col++)
                        if(row==col)
                                S[row][col]=1;
                        else
                                S[row][col]=0;


        while(fabs(maxa)>EPSL)
        {
                maxa=a[0][1];
                for(row=0;row<size;row++)
                        for(col=row+1;col<size;col++)
                                if(fabs(a[row][col])>maxa)
                                {
                                        maxa=a[row][col];//找出非主对角线上绝对值最大的元素
                                        p=row;q=col;//记录maxa的行与列数
                                }

                if(a[p][p]=a[q][q])
                        tan=PI/2.0;
                else
                        tan=2*a[p][q]/(a[p][p]-a[q][q]);

                alfa=atan(tan);//求出旋转角度

                R[p][p]=cos(alfa);
                R[q][q]=cos(alfa);
                R[p][q]=-sin(alfa);
                R[q][p]=sin(alfa);//形成矩阵R

                for(row=0;row<size;row++)
                        for(col=0;col<size;col++)
                                RT[col][row]=R[row][col];//形成R矩阵的转置矩阵RT

                for(row=0;row<size;row++)
                        for(col=0;col<size;col++)
                                for(k=0;k<size;k++)
                                        a[row][col]+=RT[row][k]*a[k][col];      //计算RT与a的乘积,并赋给a

                for(row=0;row<size;row++)
                        for(col=0;col<size;col++)
                                for(k=0;k<size;k++)
                                        a[row][col]+=a[row][k]*R[k][col];      //计算a与R的乘积,并赋给a

                for(row=0;row<size;row++)
                        for(col=0;col<size;col++)
                                for(k=0;k<size;k++)
                                        S[row][col]+=S[row][k]*R[k][col];      //计算S与R的乘积,并赋给S
        }

        fp=fopen("sc1.txt","w");

        for(row=0;row<size;row++)
        {
                fprintf(fp,"\n");
                for(col=0;col<size;col++)
                        fprintf(fp,"%10.4f",a[row][col]);
        }
        fclose(fp);

        fp=fopen("sc2.txt","w");

        for(row=0;row<size;row++)
        {
                fprintf(fp,"\n");
                for(col=0;col<size;col++)
                        fprintf(fp,"%10.4f",S[row][col]);
        }
        fclose(fp);
}


/*for(row=0;row<size;row++)//验证用语句
        {
                printf("\n");
                for(col=0;col<size;col++)
                        printf("%10.2f",S[row][col]);
        }
        printf("\n");*/

[ 本帖最后由 bobnie 于 2008-6-21 21:55 编辑 ]

dragon_

先把你用到的文件发上来再说。
哦对了，你还是新手，发不了附件……
文件应该是文本文件吧？贴上来。

[ 本帖最后由 dragon_ 于 2008-6-21 18:42 编辑 ]

bobnie

3
4 0 0
0 3 1
0 1 3

bobnie

关键时刻啊，明天就要上交了，这是最后一步了，虽然很多地方不够好，但毕竟是自己写的程序，老师问到什么地方能回答上来，不想却遇到了死循环，唉。。。。
上面数据第一个3表示数据矩阵是3阶的，下面是一个实对称矩阵。

dragon_

在Dev-C++上你的程序完全正常。只是最后的fscanf不对，应该是fprintf。

输出如下：
4.000000 0.000000 0.000000
0.000000 9.494414 -19.809553
0.000000 -19.809557 60.774910

1.000000 0.000000 0.000000
0.000000 0.681334 -3.946527
0.000000 1.670517 -1.484628

不过……貌似你只算了一步？

bobnie

程序已经调试成功，有需要的请留下E mail