main.cpp

研究生期间上《数值计算方法》课的作业

#include <iostream>
#include <stdlib.h>
#include <math.h>
#define MaxD 10

using namespace std;

float Matrix[MaxD][MaxD];
float Vk[2][MaxD];
float Uk[2][MaxD];
float Lk;
int N;
float e;

//使用幂法计算矩阵最大特征值及对应的特征向量
void CalcuMax()
{
    int i=0,j=0,k=0;           //计算次数 
    int kl = 0,kr = 0;         //控制Vk和Uk的行数 
    float temp = 0,TempLk = 0;
    do
    {
        //求Vk 
        kl = k%2;
        k = k+1;
        kr = k%2;
        for(i=0; i<N; i++)
        {
            temp = 0;
            for(j=0; j<N; j++)
            {
                temp = temp + Matrix[i][j]*Uk[kl][j];
            }
            Vk[kr][i] = temp;
        }
        
        //求本次运算中的最大Vk 
        TempLk = Lk;
        Lk = Vk[kr][0];
        for(i=0; i<N; i++)
        {
            if(fabs(Vk[kr][i]) > fabs(Lk))
            {
                Lk = Vk[kr][i];
            } 
        }
        //求Uk 
        for(i=0; i<N; i++) 
        {
            Uk[kr][i] = Vk[kr][i]/Lk;
        }
        
        //打印本次运算结果
        printf("计算次数:%d ,本次计算中:\n",k);
        printf("Lk = %5f\n",Lk);
        printf("Vk = {");
        for(i=0; i<N; i++)
        {
            if(i==N-1)
                printf("%5f",Vk[kr][i]);
            else
                printf("%5f, ",Vk[kr][i]);
        } 
        printf("}\n");
        
        printf("Uk = {");
        for(i=0; i<N; i++)
        {
            if(i==N-1)
                printf("%5f",Uk[kr][i]);
            else
                printf("%5f, ",Uk[kr][i]);
        } 
        printf("}\n\n");
                
    }while(fabs(Lk-TempLk) > e);
    printf("最大特征值为: %5f ;\n",Lk);
    printf("对应的特征向量为: ");
    printf("{");
    for(i=0; i<N; i++)
    {
        if(i==N-1)
            printf("%5f",Uk[kr][i]);
        else
            printf("%5f, ",Uk[kr][i]);
    } 
    printf("}\n");
};

//-------------------------------------------------------------------------
//以下为对矩阵进行DOOLITTLE分解,以求解方程的解 
float Sumof(int cyc, int i, int j)
{
    float fReturnVal = 0;
    for(int k=0; k<cyc; k++)
    {
        fReturnVal += Matrix[i][k]*Matrix[k][j];
    }
    return fReturnVal;
}
void Doolittle()
{
    int i = 0,j=0;
    if(Matrix[0][0] == 0.0)
    {
        printf("不能进行LU分解");
        return; 
    }
    
    for(i=1; i<N; i++)
    {
        Matrix[i][0] = Matrix[i][0]/Matrix[0][0];
    }
    
    for(i=1; i<N; i++)
    {
        Matrix[i][i] = Matrix[i][i] - Sumof(i,i,i);
        if(Matrix[i][i] == 0.0)
        {
            printf("不能进行LU分解");
            return;
        }
        for(j=i+1; j<N; j++)
        {
            Matrix[i][j] = Matrix[i][j] - Sumof(i,i,j);
            Matrix[j][i] = (Matrix[j][i] - Sumof(i,j,i))/Matrix[i][i];
        }
    }
    /*
    for(i=0; i<N; i++)
    {
        for(j=0; j<N; j++)
        {
            printf("Matrix[%d][%d]=%5f    ",i,j,Matrix[i][j]);
        }
        printf("\n");
    }
    */
}
//DOOLITTLE分解 
//--------------------------------------------------------------------------
//计算Ax=b   A=LU    LUx = b    Ly = b  Ux = y 
void CalculateVk(int now,int next)
{
    int i=0,j=0;
    float temp[MaxD];
    float ft = 0;
    //Ly = b
    temp[0] = Uk[now][0];
    for(i=1; i<N; i++)
    {
        ft = 0;
        for(j=0; j<i; j++)
        {
            ft += Matrix[i][j]*temp[j]; 
        }
        temp[i] = (Uk[now][i]-ft)/Matrix[i][i];
    }
    //Ux = y
    Vk[next][N-1] = temp[N-1]/Matrix[N-1][N-1];
    for(i=N-2; i>=0; i--)
    {
        ft = 0;
        for(j=i+1; j<N; j++)
        {
            ft += Matrix[i][j]*Vk[next][j]; 
        }
        Vk[next][i] = (temp[i]-ft)/Matrix[i][i];
    }
}

//使用反幂法计算矩阵的最小特征值和相应的特征向量 
void CalcuMin()
{
    int i=0,j=0,k=0;           //计算次数 
    int kl = 0,kr = 0;         //控制Vl和Uk的行数 
    int index = 0;
    float temp = 0,TempLk = 0;
    
    Doolittle();
    //计算最小的特征值及其对应的特征向量 
    k = 0;  
    do
    {
        //求Vk 
        kl = k%2;
        k = k+1;
        kr = k%2;
        
        CalculateVk(kl,kr);
        //求本次运算中的最大Vk 
        TempLk = Lk;
        Lk = Vk[kr][0];
        for(i=0; i<N; i++)
        {
            if(fabs(Vk[kr][i]) > fabs(Lk))
            {
                Lk = Vk[kr][i];
            } 
        }
        //求Uk 
        for(i=0; i<N; i++) 
        {
            Uk[kr][i] = Vk[kr][i]/Lk;
        }
        
        //打印本次运算结果
        printf("计算次数:%d ,本次计算中:\n",k);
        printf("Lk = %5f\n",Lk);
        printf("Vk = {");
        for(i=0; i<N; i++)
        {
            if(i==N-1)
                printf("%5f",Vk[kr][i]);
            else
                printf("%5f, ",Vk[kr][i]);
        } 
        printf("}\n");
        
        printf("Uk = {");
        for(i=0; i<N; i++)
        {
            if(i==N-1)
                printf("%5f",Uk[kr][i]);
            else
                printf("%5f, ",Uk[kr][i]);
        } 
        printf("}\n\n");
                
    }while(fabs(Lk-TempLk) > e);
    printf("最小特征值为: %5f ;\n",(1/Lk));
    printf("对应的特征向量为: ");
    printf("{");
    for(i=0; i<N; i++)
    {
        if(i==N-1)
            printf("%5f",Uk[kr][i]);
        else
            printf("%5f, ",Uk[kr][i]);
    } 
    printf("}\n");
};

int main(int argc, char *argv[])
{
    //初始化相关参数 
    int i=0,j=0;         //循环变量 
    N = 3;               //矩阵维数 
    e = 0.001;           //误差系数 
    Lk = 0;              //最终特征值 
    
    do
    {
        printf("请输入矩阵维数 不大于10\n");
        scanf("%d",&N);
    }while(N > 10);
  
    printf("请输入允许误差 e=");
    scanf("%f",&e);  
  
    printf("请输入矩阵元素:\n");
    for(i=0; i<N; i++)
    {
        for(j=0; j<N; j++)
        {
            printf("Matrix[%d][%d]=",i,j);
            scanf("%f",&Matrix[i][j]);
        }
    }
  
    for(i=0; i<N; i++)
    {
        for(j=0; j<N; j++)
        {
            printf("Matrix[%d][%d]=%5f    ",i,j,Matrix[i][j]);
        }
        printf("\n");
    }
  
    for(i=0; i<MaxD; i++)
    {
        Vk[0][i] = 1.0;
        Uk[0][i] = 1.0;
    }  
    printf("-------------------------------------------\n");
    CalcuMax(); 
    printf("-------------------------------------------\n");
    for(i=0; i<MaxD; i++)
    {
        Vk[0][i] = 1.0;
        Uk[0][i] = 1.0;
    }  
    CalcuMin();
    system("PAUSE");	
    return 0;
}