數值分析實驗之矩陣特征值（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\n',i);
fprintf('主特征值為%g\n',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!\n");
        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;
}

 

運行結果：