数值分析实验之矩阵特征值（Python代码）

一、实验目的

 1.求矩阵的部分特征值问题具有十分重要的理论意义和应用价值；

 2.掌握幂法、反幂法求矩阵的特征值和特征向量以及相应的程序设计；

 3.掌握矩阵QR分解

二、实验原理

  幂法是一种计算矩阵主特征值(矩阵按模最大的特征值)及对应特征向量的迭代方法, 特别是用于大型稀疏矩阵。设实矩阵A=[aij]n×n有一个完全的特征向量组，其特征值为λ1 ，λ2 ,…,λn,相应的特征向量为x1 ,x2 ,…,xn.已知A的主特征值是实根，且满足条件

      |λ1 |>|λ2 |≥|λ3 |≥…≥|λn |

现讨论求λ1 的方法。

幂法的基本思想是任取一个非零的初始向量ν0,由矩阵A构造一向量序列，称为迭代向量。由假设，ν0 可表示为

       ν0 =α1 x1 +α2 x2 + … +αn xn (α≠0 ),

于是得到序列vk=Avk-1,序列νk /λ1k 越来越接近A的对应于λ1 的特征向量

三、实验内容

选取五级矩阵如下：

    

  

四、实验要求

  利用幂法、反幂法求某个5阶矩阵的主特征值和特征向量，利用QR分解求一个5阶矩阵的所有特征值和特征向量

五、实验代码

  幂法（Python）

#-*- coding:utf-8 -*-
import numpy as np


def Solve(mat, max_itrs, min_delta):
    """
    mat 表示矩阵
    max_itrs 表示最大迭代次数
    min_delta 表示停止迭代阈值
    """
    itrs_num = 0
    delta = float('inf')
    N = np.shape(mat)[0]
    # 所有分量都为1的列向量
    x = np.ones(shape=(N, 1))
    #x = np.array([[0],[0],[1]])
    while itrs_num < max_itrs and delta > min_delta:
        itrs_num += 1
        y = np.dot(mat, x)
        #print(y)
        m = y.max()
        #print("m={0}".format(m))
        x = y / m
        print("***********第{}次迭代*************".format(itrs_num))
        print("y = ",y)
        print("m={0}".format(m))
        print("x^T为：",x)
        print(
            "——————————————分割线——————————————")


IOS = np.array([[2, 3, 4, 5, 6],
                [4, 4, 5, 6, 7],
                [0, 3, 6, 7, 8],
                [0, 0, 2, 8, 9],
                [0, 0, 0, 1, 0]], dtype=float)
Solve(IOS, 100, 1e-3)

   运行结果：

       

  

  反幂法（MATLAB）

A =[2,3,4,5,6;4,4,5,6,7;0,3,6,7,8;0,0,2,8,9;0,0,0,1,0];
I = eye(5,5);
p=13;
u0 = [1;1;1;1;1];
v = inv(A - p * I) * u0;
u = v / norm(v, inf);
i = 0;
while norm(u - u0, inf) > 1e-5
    u0 = u;
    v = inv(A - p * I) * u0;
    u = v / norm(v, inf);
    i=i+1;
end
x = p + 1 / norm(v, inf);
fprintf('迭代次数为%g
',i);
fprintf('主特征值为%g
',x);
u

　　运行结果：

           

 QR分解（C）

void QR(Matrix* A, Matrix* Q, Matrix* R)
{
    unsigned  i, j, k, m;
    unsigned size;
    const unsigned N = A->row;
    matrixType temp;

    Matrix a, b;

    if (A->row != A->column || 1 != A->height)
    {
        printf("ERROE: QR() parameter A is not a square matrix!
");
        return;
    }

    size = MatrixCapacity(A);
    if (MatrixCapacity(Q) != size)
    {
        DestroyMatrix(Q);
        SetMatrixSize(Q, A->row, A->column, A->height);
        SetMatrixZero(Q);
    }
    else
    {
        Q->row = A->row;
        Q->column = A->column;
        Q->height = A->height;
    }

    if (MatrixCapacity(R) != size)
    {
        DestroyMatrix(R);
        SetMatrixSize(R, A->row, A->column, A->height);
        SetMatrixZero(R);
    }
    else
    {
        R->row = A->row;
        R->column = A->column;
        R->height = A->height;
    }

    SetMatrixSize(&a, N, 1, 1);
    SetMatrixSize(&b, N, 1, 1);

    for (j = 0; j < N; ++j)
    {
        for (i = 0; i < N; ++i)
        {
            a.array[i] = b.array[i] = A->array[i * A->column + j];
        }

        for (k = 0; k < j; ++k)
        {
            R->array[k * R->column + j] = 0;

            for (m = 0; m < N; ++m)
            {
                R->array[k * R->column + j] += a.array[m] * Q->array[m * Q->column + k];
            }

            for (m = 0; m < N; ++m)
            {
                b.array[m] -= R->array[k * R->column + j] * Q->array[m * Q->column + k];
            }
        }

        temp = MatrixNorm2(&b);
        R->array[j * R->column + j] = temp;

        for (i = 0; i < N; ++i)
        {
            Q->array[i * Q->column + j] = b.array[i] / temp;
        }
    }

    DestroyMatrix(&a);
    DestroyMatrix(&b);
}

void Eigenvectors(Matrix* eigenVector, Matrix* A, Matrix* eigenValue)
{
    unsigned i, j, q;
    unsigned count;
    int m;
    const unsigned NUM = A->column;
    matrixType eValue;
    matrixType sum, midSum, mid;
    Matrix temp;

    SetMatrixSize(&temp, A->row, A->column, A->height);

    for (count = 0; count < NUM; ++count)
    {
        //计算特征值为eValue，求解特征向量时的系数矩阵
        eValue = eigenValue->array[count];
        CopyMatrix(&temp, A, 0);
        for (i = 0; i < temp.column; ++i)
        {
            temp.array[i * temp.column + i] -= eValue;
        }

        //将temp化为阶梯型矩阵
        for (i = 0; i < temp.row - 1; ++i)
        {
            mid = temp.array[i * temp.column + i];
            for (j = i; j < temp.column; ++j)
            {
                temp.array[i * temp.column + j] /= mid;
            }

            for (j = i + 1; j < temp.row; ++j)
            {
                mid = temp.array[j * temp.column + i];
                for (q = i; q < temp.column; ++q)
                {
                    temp.array[j * temp.column + q] -= mid * temp.array[i * temp.column + q];
                }
            }
        }
        midSum = eigenVector->array[(eigenVector->row - 1) * eigenVector->column + count] = 1;
        for (m = temp.row - 2; m >= 0; --m)
        {
            sum = 0;
            for (j = m + 1; j < temp.column; ++j)
            {
                sum += temp.array[m * temp.column + j] * eigenVector->array[j * eigenVector->column + count];
            }
            sum = -sum / temp.array[m * temp.column + m];
            midSum += sum * sum;
            eigenVector->array[m * eigenVector->column + count] = sum;
        }

        midSum = sqrt(midSum);
        for (i = 0; i < eigenVector->row; ++i)
        {
            eigenVector->array[i * eigenVector->column + count] /= midSum;
        }
    }
    DestroyMatrix(&temp);
}
int main()
{
    const unsigned NUM = 50; //最大迭代次数

    unsigned N = 5;
    unsigned k;

    Matrix A, Q, R, temp;
    Matrix eValue;

    //分配内存
    SetMatrixSize(&A, N, N, 1);
    SetMatrixSize(&Q, A.row, A.column, A.height);
    SetMatrixSize(&R, A.row, A.column, A.height);
    SetMatrixSize(&temp, A.row, A.column, A.height);
    SetMatrixSize(&eValue, A.row, 1, 1);

    A.array[0] = 2;A.array[1] = 3;A.array[2] = 4;A.array[3] = 5;A.array[4] = 6;
    A.array[5] = 4;A.array[6] = 4;A.array[7] = 5;A.array[8] = 6;A.array[9] = 7;
    A.array[10] = 0;A.array[11] = 3;A.array[12] = 6;A.array[13] = 7;A.array[14] = 8;
    A.array[15] = 0;A.array[16] = 0;A.array[17] = 2;A.array[18] = 8;A.array[19] = 9;
    A.array[20] = 0;A.array[21] = 0;A.array[22] = 0;A.array[23] = 1;A.array[24] = 0;

    //拷贝A矩阵元素至temp
    CopyMatrix(&temp, &A, 0);

    //初始化Q、R所有元素为0
    SetMatrixZero(&Q);
    SetMatrixZero(&R);
    //使用QR分解求矩阵特征值
    for (k = 0; k < NUM; ++k)
    {
        QR(&temp, &Q, &R);
        MatrixMulMatrix(&temp, &R, &Q);
    }
    //获取特征值，将之存储于eValue
    for (k = 0; k < temp.column; ++k)
    {
        eValue.array[k] = temp.array[k * temp.column + k];
    }

    //对特征值按照降序排序
    SortVector(&eValue, 1);

    //根据特征值eValue，原始矩阵求解矩阵特征向量Q
    Eigenvectors(&Q, &A, &eValue);

    //打印特征值
    printf("特征值：");
    PrintMatrix(&eValue);

    //打印特征向量
    printf("特征向量");
    PrintMatrix(&Q);
    DestroyMatrix(&A);
    DestroyMatrix(&R);
    DestroyMatrix(&Q);
    DestroyMatrix(&eValue);
    DestroyMatrix(&temp);

    return 0;
}

运行结果：