精华内容
下载资源
问答
  • 基于微分方程组理论和矩阵理论,采用按列比较方法和待定矩阵方法,给出了非齐次项为二次多项式与指数函数乘积的一类三维二阶常系数线性微分方程组的特解公式。对特殊情况进行了讨论,并通过算例验证了微分方程组特解...
  • 线性代数,多项式

    复制以下的代码,选择运行的函数,验证有关矩阵、多项式的操作。
    相关参考:Python语言程序设计(上海交通大学出版社 赵璐主编)<<----传送门
    原谅我没有找到电子版,不然一定爬下来

    准备工作:恶补一下矩阵相关内容,高中课本
    在这里插入图片描述
    解释和输出都在注释中

    # -*- coding: utf-8 -*-
    import numpy as np
    
    
    # 线性代数的应用-----------------------------------------------------------------------------------
    def linearAlgebra():
        """
        linalg.det()            计算行列式
    
        linalg.inv()            计算逆矩阵
    
        linalg.solve()          多元一次方程组求根
    
        linalg.eig()            返回特征值和特征向量构成的元组
    
        linalg.eigvals()        计算特征值
    
        linalg.svd()            矩阵的奇异值分解
    
        linalg.pinv()           广义逆矩阵
    
        """
    
        Array1 = np.arange(6).reshape((2, 3))
        Array2 = np.arange(6).reshape((3, 2))
        print(Array1.dot(Array2))  # 矩阵内积
        # [[0 1 2]
        #  [3 4 5]]
    
        Array3 = np.arange(4).reshape((2, 2))
        detArray3 = np.linalg.det(Array3)  # 计算行列式
        print(detArray3)  # -2.0
    
        invArray3 = np.linalg.inv(Array3)  # 计算逆矩阵
        print(invArray3)
        # [[-1.5  0.5]
        #  [1.    0.]]
    
        eyeArray = Array3.dot(invArray3)  # 互逆矩阵内积运算,验证结果
        print(eyeArray)
        # [[1. 0.]
        #  [0. 1.]]
    
        Array4 = np.array([[-1, 0], [0, 1]])
        eigenvalues, eigenvectors = np.linalg.eig(Array4)  # 计算特征值和特征向量
        print(eigenvalues)  # 特征值构成的数组,[-1.  1.]
    
        print(eigenvectors)  # 特征向量构成的数组
        # [[1. 0.]
        #  [0. 1.]]
    
        # 验证特征值和特征向量,满足:Aα=λα
        # 其中A是n阶方阵,λ是特征值,非零向量α是矩阵A属于特征值λ的特征向量
        print(Array4.dot(eigenvectors[0]) ==
              eigenvalues[0]*eigenvectors[0])  # [ True  True]
        print(Array4.dot(eigenvectors[1]) ==
              eigenvalues[1]*eigenvectors[1])  # [ True  True]
    
        # 在线性代数中,矩阵可以对向量进行线性变换,这对数学中的线性方程组。
        # numpy.linalg中的sovle()函数可以求解形如Ax=b的线性方程组,其中A为
        # 系数矩阵,b为一维或二维的数组,x是未知量
        # 如求解鸡兔同笼问题:上有35头,下有94足,求鸡兔个数
    
        heads, foots = 35, 94
        A = np.array([[1, 1], [2, 4]])  # 方程组的系数矩阵
        b = np.array([heads, foots])  # 方程组右侧的常数矩阵
        X = np.linalg.solve(A, b)  # sovle()函数返回方程组的解
        print('鸡:{},兔:{}'.format(X[0], X[1]))  # 鸡:23.0,兔:12.0
    
    
    # 多项式的应用-------------------------------------------------------------------------------------------
    def polynomial():
        """  
        poly1d(A)               利用系数数组A生成多项式
    
        polyval(p,k)            求多项式p在x=k时的值
    
        polyder(p,m=1)          求多项式p的m阶导数,m默认值为1
    
        polyint(p,m=1)          求多项式p的m重积分,m默认值为1
    
        polyadd(p1,p2)          多项式求和,等价于p1+p2
    
        polysub(p1,p2)          多项式求差p1-p2
    
        polymul(p1,p2)          多项式求积p1*p2
    
        polydiv(p1,p2)          多项式求商p1/p2,结果为商和余数构成的元组,商和余数都用多项式表示
    
        polyfit(x,y,k)          多项式拟合,x,y分别为要拟合的两组数据,k为拟合多项式中最高次幂
    
        """
    
        A = np.array([1, 0, -2, 1])  # 系数数组,没有出现的系数项用0补齐
        f = np.poly1d(A)
        print(f)  # 输出多项式f的数学表达式
        #     3
        # 1 x - 2 x + 1
        print(type(f))  # <class 'numpy.poly1d'>
        print(f(1))  # 当x=1时,输出多项式的值  0
        print(f(2))  # 5
        print(f([1, 2]))  # [0,5]
    
        print(np.polyval(f, 1))  # 当x=1时,输出多项式的值  0
        print(np.polyval(f, 2))  # 5
        print(np.polyval(f, [1, 2]))  # [0,5]
    
        fder1 = np.polyder(f)  # 求多项式f的一阶导数
        print(fder1)
        #     2
        # 3 x - 2
        fder2 = np.polyder(f, 2)  # 求多项式f的二阶导数
        print(fder2)  # 6 x
    
        fint1 = np.polyint(f)  # 求多项式f的一重积分
        print(fint1)
        #        4     2
        # 0.25 x - 1 x + 1 x
        print(f == np.polyder(fint1))  # 求一重积分的一阶导数来验证是否正确  True
        fint2 = np.polyint(f, 2)  # 求多项式f的二重积分
        print(fint2)
        #        5          3       2
        # 0.05 x - 0.3333 x + 0.5 x
    
        p1 = np.poly1d(np.array([1, 2, 3]))
        p2 = np.poly1d(np.array([1, 2, 3, 4]))
        print(p1)
        #     2
        # 1 x + 2 x + 3
        print(p2)
        #     3     2
        # 1 x + 2 x + 3 x + 4
        print(p1+p2)
        #     3     2
        # 1 x + 3 x + 5 x + 7
        print(np.polyadd(p1, p2))
        #     3     2
        # 1 x + 3 x + 5 x + 7
    
        print(p2-p1)
        #     3     2
        # 1 x + 1 x + 1 x + 1
        print(np.polysub(p2, p1))
        #     3     2
        # 1 x + 1 x + 1 x + 1
    
        print(np.polymul(p1, p2))  # 等价于p1*p2
        #     5     4      3      2
        # 1 x + 4 x + 10 x + 16 x + 17 x + 12
    
        print(np.polydiv(p2, p1))  # 等价于p2/p1
        #(poly1d([1., 0.]), poly1d([4.]))
    
        x = np.array([0.0, 1.0, 2.0, 3.0, 4.0, 5.0])
        y = np.array([0.0, 0.8, 0.9, 0.1, -0.8, -1.0])
        parray = np.polyfit(x, y, 3)  # 用 polyfit 返回一个拟合多项式的系数数组
        print(parray)
        # [0.08703704 - 0.81349206  1.69312169 - 0.03968254]
        p = np.poly1d(parray)
        print(p)
        #     3          2
        # 0.08704 x - 0.8135 x + 1.693 x - 0.03968
    
    
    # 选择要验证的函数-----------------------------------------------------------------------------------------
    def main():
        # linearAlgebra()
        # polynomial()
    
    
    if __name__ == "__main__":
        main()
    
    

    相关参考:Python语言程序设计(上海交通大学出版社 赵璐主编)<<----传送门

    如有错误,欢迎私信纠正
    技术永无止境,谢谢支持!

    展开全文
  • 该函数计算拉格朗日多项式相对于节点基“InterpolationPoints”的二阶导数矩阵。 在点“EvaluationPoints”处评估多项式多项式按列堆叠,行与评估点相对应。 假设输入是一维数组。
  • 在多项域里构造正交小波滤波器组等价于构造仿酉矩阵,而仿酉矩阵的构造涉及非线性方程组...通过对 Cayley 变换的研究,给出构造仿酉矩阵的条件,利用这些条件给出元素为二元一次多项式二阶仿酉矩阵的通解,并给出了算例.
  • 矩阵论笔记】最小多项式与Jordan型的关系

    千次阅读 多人点赞 2020-05-07 17:02:01
    Jordan块的最小多项式是他的特征多项式,阶数不能再降了。 例题 A1A_1A1​是一个Jordan块,他的最小多项式不将阶。 这时候A2A_2A2​的最小多项式只有两种情况,从1阶开始试。然后取最小共倍数,构成原来矩阵的.....

    最小多项式

    方阵A的次数最低、且首一的零化多项式称为A的最小多项式。
    在这里插入图片描述

    最小多项式的一般形式

    在这里插入图片描述
    算这个没什么办法,只能暴力计算,从m=1开始算,把A带进去是不是等0。

    Jordan块的最小多项式是他的特征多项式,阶数不能再降了。

    在这里插入图片描述

    例题

    在这里插入图片描述
    在这里插入图片描述
    A 1 A_1 A1是一个Jordan块,他的最小多项式不将阶。
    在这里插入图片描述
    这时候 A 2 A_2 A2的最小多项式只有两种情况,从1阶开始试。然后取最小共倍数,构成原来矩阵的最小多项式。
    在这里插入图片描述
    有了最小多项式就可以简化方阵多项式的计算。

    定理:相似矩阵有相似的最小多项式
    N阶方阵可对角化的充要条件是最小多项式无重根。

    最小多项式与Jordan的关系

    在这里插入图片描述

    性质

    在这里插入图片描述

    上面的 n i n_i ni是代数重数,是子Jordan的阶数。
    子Jordan矩阵是由Jordan块构成,Jordan块的个数是由 λ i \lambda_i λi的几何重数构成 k i k_i ki

    在这里插入图片描述
    这些Jordan块的阶和最小多项式有关,最小多项式的幂就是Jordan块的阶。
    在这里插入图片描述

    例题

    代数重数就是特征多项式的幂,阶数就是总的幂数之和。
    在这里插入图片描述
    在这里插入图片描述
    根据最小多项式得出子Jordan块的最高阶数
    在这里插入图片描述
    在这里插入图片描述

    简化方阵多项式的计算

    1、相似化简成对角阵
    2、通过零化多项式对多项式降阶

    展开全文
  • 第8章 矩阵特征值问题计算

    千次阅读 2012-06-05 21:59:08
    第8章 矩阵特征值问题计算 ...物理、力学和工程技术中的很多问题在数学上都归结为求矩阵特征值问题。例如,振动问题(大型桥梁或建筑物...为的特征多项式. 的特征方程 (8.1.1) 一般有个根(实的或复的
     
      
       
        
       

    8 矩阵特征值问题计算

    8.1

    物理、力学和工程技术中的很多问题在数学上都归结为求矩阵特征值问题。例如,振动问题(大型桥梁或建筑物的振动、机械的振动、电磁振荡等),物理学中某些临界值的确定,这些问题都归结为下述数学问题。

    定义1 (1)已知,则称

    特征多项式

    的特征方程

    (8.1.1)

    一般有个根(实的或复的,重根按重数计算.当时,为实系数次代数方程,其复根共轭成对出现),称为的特征值.用表示的所有特征值的集合,并称之为的谱集.

    (2)设的特征值,相应的齐次方程组

                 (8.1.2)

    的非零解称为的对应于特征向量

    例1 求的特征值及特征向量,其中

     矩阵的特征方程为

    求得的特征值为:.对应的各特征向量是对应齐次方程

    的解.计算得到:

    下面我们来看几个有关特征值的结论.

    定理8.1.1 设的特征值,且,其中0,则

    (1)特征值(为非零常数);

    (2)的特征值,即

    (3)的特征值;

    (4)设为非奇异矩阵,那么的特征值,即

    定理8.1.2 设阶矩阵的特征值,则

    (1) 是矩阵的迹);

    (2)

    定理8.1.3 设,则

    定理8.1.4 设为分块上三角阵,即

    其中每个对角块均为方阵,则

    定理8.1.5 设为相似矩阵(即存在非奇异矩阵使),则

    (1)有相同的特征值;

    (2)如果的特征向量,则的特征向量.

    定理8.1.5说明,一个矩阵经过相似变换,则的特征值不变.

    定义8.1.2 设,如果有一个重数为的特征值且对应于矩阵的线性无关的特征向量个数少于,称亏损矩阵

    一个亏损矩阵是一个没有足够特征向量的矩阵,亏损矩阵在理论上和计算上都存在困难.

    定理8.1.6(1)可对角化,即存在非奇异矩阵使

    的充要条件是具有个线性无关的特征向量.

    (2)如果个(不同的特征值,则对应的特征向量线性无关.

    定理8.1.7(对称矩阵的正交约化)设为对称矩阵,则:

    (1)的特征值均为实数;

    (2)个线性无关的特征向量;

    (3)存在一个正交矩阵使

    特征值,而的列向量的对应于的特征向量.

    下面讨论矩阵特征值界的估计.

    定义8.1.3 设.令:

    (1)

    (2)集合.称复平面上以为圆心,以为半径的所有圆盘为Gerschgorin圆盘.

    定理8.1.8Gerschgorin圆盘定理)(1)设,则的每一个特征值必属于下述某个圆盘之中

    或者说,的特征值都在复平面上个圆盘的并集中.

    (2)如果个圆盘组成一个连通的并集,且与余下个圆盘是分离的,则内恰包含的一个特征值.

    证明 只就(1)给出证明,设的特征值,即有使

    记 ,考虑的第个方程,即

    于是

    这说明,的每一个特征值必位于的一个圆盘中,并且相应的特征值一定位于第个圆盘中(其中是对应特征向量绝对值最大的分量的下标).

    利用相似矩阵性质,有地可以获得的特征值进一步的估计,即适当选取非奇异对角阵

    并做相似变换.适当选取可使某些圆盘半径及连通性发生变化.

    例2 估计矩阵

    特征值的范围.

     的3个圆盘为

    由定理8.1.8,可知的3个特征值位于3个圆盘的并集中,由于是弧立圆盘,所以内恰好包含的一个特征值(为实特征值),即

    的其他两个特征值包含在的并集中.

    现选取对角阵

    做相似变换

    的3个圆盘为

    显然,3个圆盘都是孤立圆盘,所以,每一个圆盘都包含的一个特征值(为实特征值)且有估计

    下面给出理论上有关通过酉相似变换及正交相似变换可以约化一般矩阵到什么程序的问题.

    定理8.1.9Schur定理) 设,则存在酉阵使

    (上三角阵)

    其中的特征值.

    时,如果限制用正交相似变换,由于有复的特征值,不能用正交相似变换约化为上三角阵.用正交相似变换能将约化到什么程度呢?

    定理8.1.10(实Schur分解) 设,则存在正交矩阵使

    其中对角块为一阶或二阶方阵,且每一个一阶的实特征值,每个二阶对角块的两个特征值是的两个共轭复特征值.

    我们转向实Schur型的实际计算.

    定义8.1.4 设阶实对称矩阵,对于任一非零向量,称

    为对应于向量瑞利Rayleigh)商.

    定理8.1.11 设为对称矩阵(其特征值次序记为),则

    1 (对任何非零)

    2

    3

    证明 只证1,关于2,3留作习题.

    由于为实对称矩阵,可将对应的特征向量正交规范化,则有.设中任一向量,则有展开式

    于是

    从而1成立,结论1说明瑞利商必位于之间.

    关于计算矩阵的特征值问题,当时,我们还可按行列式展开的办法求的根.但当较大时,如果按展开行列式的办法,首先求出的系数,再求的根,工作量就非常大,用这种办法求矩阵特征值是不切实际的,由此需要研究求的特征值及特征向量的数值解法.

    本章将介绍一些计算机上常用的两类方法,一类是幂法及反幂法(迭代法),另一类是正交相似变换的方法(变换法).

    8.2 幂法及反幂法

    8.2.1 幂法

    幂法是一种计算矩阵主特征值(矩阵按模最大的特征值)及对应特征向量的迭代方法,特别适用于大型稀疏矩阵。反幂法是计算海森伯格或三对角阵的对应一个给定近似特征值的特征向量的有效方法之一。

    设实矩阵有一个完全的特征向量组,其特征值为,相应的特征向量为.已知的主特征值是实的,且满足条件

    ,  (8.2.1)

    现讨论求的方法.

    幂法的基本思想是任取一个非零的初始向量,由矩阵构造一向量序列

    (8.2.2)

    称为迭代向量.由假设,可表示为

    ,    (8.2.3)

    于是

    其中,显然

    由假设知,故

    (8.2.4)

    这说明序列 越来越接近的对应于的特征向量,或者说当充分大时

    (8.2.5)

    即迭代向量的特征向量的近似向量.

    下面再考虑主特征值的计算,用表示的第个分量,则

    ,        (8.2.6)

    因此

    ,           (8.2.7)

    也就是说两相邻迭代向量分量的比值收敛到主特征值.

    这种由已知非零向量及矩阵的幂构造向量序列{}以计算的主特征值(利用(8.2.7)式)及相应特征向量(利用(8.2.5)式)的方法称为幂法

    (8.2.6)知,的收敛速度由比值来确定,越小收敛越快,当接近于1时收敛会很慢.

    定理8.2.1 设个线性无关的特征向量,主特征值满足

    则对任何非零初始向量(8.2.4), (8.2.7)式成立.

    如果的主特征值为实的重根,即,且

    又设个线性无关的特征向量,对应的个线性无关特征向量为,则由(8.2.2)

    这说明当的主特征值是实的重根时,定理5的结论还是正确的.

    应用幂法计算的主特征值及对应的特征向量时,如果>1(或<1),迭代向量的各个不等于零的分量将随而趋向于无穷(或趋向于零),这样在计算机实现时就可能 “溢出”.为了克服这个缺点,就需要将迭代向量加以规范化.

    设有一向量,将其规范化得到向量

    其中表示向量的绝对值最大的分量,即如果有

    ,且为所有绝对值最大的分量中的最小下标.

    在定理8.2.1的条件下幂法可这样进行:任取一初始向量,构造向量序列

    (8.2.3)

      (8.2.8)

    这说明规范化向量序列收敛到主特征值对应的特征向量.

    同理,可得到

    收敛速度由比值确定.总结上述讨论,有

    定理8.2.2 设个线性无关的特征向量,主特征值满足,则对任何非零初始向量,按下述方法构造的向量序列

           (8.2.9)

    则有

    (1)

    (2)

    例3 用幂法计算

    的主特征值和相应的特征向量.计算过程如表8-1.

    [解]下述结果是用Visual C++6.0在微机上进行运算得到的,的分量值是舍入值.于是得到2.5365326,及相应的特征向量(0.74822631,0.64966657,1).和相应的特征向量的真值(8位数字)为

    表8-1 幂法计算结果

    K

    (规范化向量)

    0

    (1,

    1,

    1)

     

    1

    (0.90909091,

    0.81818182,

    1)

    2.7500000

    2

    (0.83760684,

    0.74358974,

    1)

    2.6590909

    3

    (0.79901559,

    0.70303527,

    1)

    2.6047009

    4

    (0.77741499,

    0.68033766,

    1)

    2.5752666

    5

    (0.76510819,

    0.66740583,

    1)

    2.5587919

    6

    (0.75802535,

    0.65996327,

    1)

    2.5494056

    7

    (0.75392531,

    0.65565500,

    1)

    2.5440035

    8

    (0.75154396,

    0.65315271,

    1)

    2.5408764

    9

    (0.75015815,

    0.65169652,

    1)

    2.5390602

    10

    (0.74935078,

    0.65084814,

    1)

    2.5380032

    11

    (0.74888009,

    0.65035355,

    1)

    2.5373874

    12

    (0.74860558,

    0.65006510,

    1)

    2.5370284

    13

    (0.74844545,

    0.64989683,

    1)

    2.5368191

    14

    (0.74835202,

    0.64979867,

    1)

    2.5366969

    15

    (0.74829751,

    0.64974139,

    1)

    2.5366257

    16

    (0.74826571,

    0.64970797,

    1)

    2.5365841

    17

    (0.74824715,

    0.64968847,

    1)

    2.5365598

    18

    (0.74823632,

    0.64967709,

    1)

    2.5365457

    19

    (0.74823000,

    0.64967045,

    1)

    2.5365374

    20

    (0.74822631,

    0.64966657,

    1)

    2.5365326

    8.2.2 加速方法

    1 原点平移法

    由前面讨论知道,应用幂法计算的主特征值的收敛速度主要由比值来决定,但当接近于1时,收敛可能很慢.这时一个补救的办法就是采用加速收敛的方法.

    引进矩阵

    其中为选择参数.设的特征值为,则的相应特征值为,而且的特征向量相同.

    如果需要计算的主特征值,就要适当选择,使仍然是的主特征值,且使

    应用幂法,使得在计算的主特征值的过程中得到加速.这种方法通常称为原点平移法.对于的特征值的某种分布,它是十分有效的.

    例4 设有特征值

    比值,作变换

    12,

    的特征值为

    应用幂法计算的主特征值的收敛速度的比值为

    虽然常常能够选择有利的值,使幂法得到加速,但设计一个自动选择适当参数的过程是困难的.

    下面考虑当的特征值是实数时,怎样选择使采用幂法计算得到加速.

    的特征值满足

       (8.2.10)

    则不管如何,的主特征值为.当我们希望计算时,首先应选择使

    且使收敛速度的比值

    显然,当,即为最小,这时收敛速度的比值为

    的特征值满足(8.2.10)能初步估计时,我们就能确定的近似值.

    当希望计算时,应选择

    使得应用幂法计算得到加速.

    例5 计算例3中矩阵的主特征值.

    作变换,设0.75,则

    应用幂法,计算结果如表8.2.

    8-2 幂法计算结果

    K

    (规范化向量)

    1

    2

    3

    4

    5

    6

    7

    8

    9

    10

    0.875

    0.78333333

    0.76834862

    0.75382166

    0.7516

    0.74913041

    0.7487941

    0.74836879

    0.74831858

    0.74824505

    0.75

    0.7

    0.66513761

    0.65796178

    0.65217778

    0.65105965

    0.65007315

    0.64989782

    0.64972854

    0.64970127

    1

    1

    1

    1

    1

    1

    1

    1

    1

    1

    2.0000000

    1.8750000

    1.8166667

    1.8004587

    1.7914013

    1.7888444

    1.7873301

    1.7869153

    1.7866589

    1.7865914

    由此得的主特征值为1.7865914,的主特征值

    2.5365914.

    与例3结果比较,上述结果比例3迭代15次还好.

    原点位移的加速方法,是一个矩阵变换方法.这种变换容易计算,又不破坏矩阵的稀疏性,但的选择依赖于对的特征值分布的大致了解.

    2 瑞利商加速

    由定理8.1.11知,对称矩阵可用瑞利商的极值来表示.下面我们将把瑞利商应用到用幂法计算实对称矩阵的主特征值的加速收敛上来.

    定理8.2.3 设为对称矩阵,特征值满足

    对应的特征向量满足,应用幂法(公式(8.2.9))计算的主特征值,则规范化向量瑞利商给出

    证明 由(8.2.8)式及

    .  (8.2.11)

    8.2.3 反幂法

    反幂法用来计算矩阵按模最小的特征值及其特征向量,也可用来计算对应与一个给定近似特征值的特征向量。

    为非奇异矩阵,的特征值依次记为

    相应的特征向量为的特征值为

    相应的特征向量为

    所以计算的按模最小的特征值的问题就是计算 的按模最大的特征值问题。

    对于应用幂法迭代(称为反幂法),可求得矩阵的主特征值1/,从而求得的按模最小的特征值

    反幂法迭代公式为:

    任取初始向量,构造向量序列

     迭代向量可以通过解方程组求得.

    定理8.2.4 设为非奇异矩阵且有个线性无关的特征向量,其对应的特征值满足

    则对任何初始非零向量,由反幂法构造的向量序列满足

    (1)

    (2)

    收敛速度的比值为

    在反幂法中也可以用原点平移法来加速迭代过程或求其他特征值及特征向量。

    如果矩阵存在, 对其应用幂法,得反幂法的迭代公式

                  2.12

    如果p的特征值的一个近似值,且设与其他特征值分离,

    (),

    就是说的主特征值,可用反幂法(2.12)计算特征值及特征向量。这样我们有下面的定理

        定理8.2.5 个线性无关的特征向量, 的特征值及对应的特征向量分别记为(i=1,2,…,n),的近似值, 存在,且

    ()

    则对任意的初始非零向量,有反幂法迭代公式(2.12)构造的向量序列满足

    (1) ;

    (2) ,

    且收敛速度由比值

    确定。

    由该定理知,对(其中)应用反幂法,可用来计算特征向量. 只要选择的的一个较好的近似且特征值分离情况较好,一般很小,常常只要迭代一二次就可完成特征向量的计算。

    反幂法迭代公式中的是通过解方程组

    求得的。为节省工作量,可先将进行三角分解 

    其中P为某个排列阵,于是求相当于解两个三角形方程组

    实验表明,按下述方法选择较好的:选使

    (2.13)

    用回代求解(2.13)即得,然后再按公式(2.12)进行迭代。

    反幂法迭代公式:

          1. 分解计算

    , 且保存信息。

          2. 反幂法迭代 

         (1)

     

    (2) k=2,3,…

    1) ,求得;
    ,求得

    2)

    3) 计算

    6 用反幂法求

    的对应于特征值λ=1.2679(精确特征值为)的特征向量.

    [] 用部分选主元的三角分解将 (其中)分解为

    其中

    ,得

    ,得

    对应的特征向量是

    由此看出的相当好的近似。

    特征值的真值为

    8.3 豪斯霍8.4 尔德方法

    8.4.1 引言

    本节讨论下面问题:

    (1) 用初等反射阵作正交相似变换约化一般实矩阵为上海森伯格阵(次对角元素下方的元素均为零)

    (2) 用初等反射阵作正交相似变换约化对称矩阵为对称三对角阵。

    于是,求原矩阵特征值问题,就转化为求上海森伯格阵或对称三对角阵的特征值问题。

    8.4.2 用正交相似变换约化一般矩阵为上海森伯格阵。

       。下面说明可选择初等反射阵使经正交相似变换约化一个上海森伯格阵。

      1)设

    其中,

    时进行下一步。当时,对维列向量待定,令

    这是阶初等反射阵(),使其满足,维向量)  

    易见

    现在令

    ,

    其中

        (2) 步约化:重复上述过程,设对已完成第1步,,第1步正交相似变换,即有

    其中阶上海森伯格阵,

    ,于是可选择初等反射阵使, 其中的计算公式为

    (8.3.2)

    ,

    (8.3.3)

    其中+1阶上海森伯格阵。第k步约化只需计算 (为对称阵时可以不算)

    (3) 重复上述过程,有

    定理8.3.1 (矩阵的上Hessenberg) , 则存在初等反射阵 使 

    (上海森伯格阵)

    算法1(Householder约化上H-矩阵) , 本算 法计算 (上海森伯格型),其中为初等反射阵的乘积。

        1.

        2. 对于

    计算初等反射阵使,即取

     2)约化计算

    3

    算法分析:  

    本算法约需要次乘法运算,要明显形成还需加次乘法。

    在计算的过程,的计算可能会引起溢出现象.为避免溢出现象,可以考虑向量的规范化,即,由计算得到.但应注意不能将保存到中.

    在实际计算时,不必要得出矩阵形式,只需要得出算法中交待的相关参数及向量即可.

        7 用豪斯霍尔德方法将

    矩阵约化为上Hessenberg阵。

       选取初等反射阵使 R1 , 其中

    计算: (规范化)

    则有 

     (2) 约化计算:令

    8.4.3 用正交相似变换约化对称阵为对称三对角阵

    定理8.3.2 (豪斯霍尔德约化对称阵为对称三对角阵) , 则存在初等反射阵 使 

    证明 由定理8.3.1,存在初等反射阵使,注意到,显然(对称矩阵),因此若对称,则必有对称,从而是对称三对角阵.

    定理证毕

    根据定理17的结论,我们容易对算法1加以改进得之如下算法2.

    算法2(豪斯霍尔得德约化对称阵为对称三对角阵)为对称矩阵,本算法确定初等反射阵使为对称三对角阵.C的对角元存放在数组内,C的次对角元存放在数组内,数组的初始值可用来存放,确定中向量的分量存放在A的相应位置.冲掉,约化A的结果冲掉A,数组A的上部分元素不变,如果第不需要变换则置为零.

    1. 对于

    (1)

    (2) 确定变换

    1) 计算

    2) 如果,则,转4

    3) 计算

    4)

    5)

    6)

    7)

    (3) 应用变换

    1)

    2) 计算

    对于

    (a)

    (b)

    3) 计算

    4) 计算对角线以下部分

    对于

    5) 继续循环

    2.

    对对称矩阵用初等反射阵正交相似约化为对称三对角阵大约需要次乘法。

    用正交矩阵进行相似约化有一些特点,如构造的容易求逆,且的元素数量级不大,这个算法是十分稳定的。

    8.5 QR

    QR方法是一种变换方法,是计算一般矩阵(中小型矩阵)全部特征值问题的最有效方法之一。目前QR方法主要用来计算:(1) 上海森伯格阵的全部特征值问题;(2) 计算对称三对角矩阵的全部特征值问题。QR方法具有收效快,算法稳定等特点。

    对于一般矩阵 (或对称矩阵), 则首先用豪斯霍尔德方法将化为上海森伯格阵(或对称三对角阵), 然后再用QR方法计算的全部特征值。

    8.5.1  基本QR算法

    ,基本QR迭代算法如下:

           8.4.1

    其中为上三角阵,为正交阵。

    定理8.4.1 基本QR算法的性质定理)记,则有

    1

    2

    3

    [证明] (1);

    (2) ;

    (3) 事实上,

    由(2)知,代入上式得

    从而便有

    [证毕]

    定理????? 设阶矩阵个特征值满足,并设阶非奇异矩阵的第行是对应于的左特征向量.如果有LU分解,则由(8.4.1)产生的矩阵的对角线以下的元素趋向于零,同时对角元素趋向于

    [证明]令:,则有.假定LU分解为,其中是单位下三角矩阵,是上三角矩阵.这样,我们有

    ???????

    其中.由于是单位下三角矩阵,而,故必有

    ???????

    现令的QR分解为:.由的非奇异,可要求的对角元素均为正负.将这一分解代入???????可得

    ???????

    充分大时,是非奇异的,故它有如下的QR分解

    ???????

    其中的对角元素均为正数.从(8.4.6)(8.4.8)

    于是有

    (8.4.9)

    ???????代入??????,得

    即已经找到了的一个QR分解.为保证分解中的上三角矩阵的对角元素均为正数,令

    其中的第个对角元素.于是我们有

    将上式与(8.4.4)比较,并注意到QR分解的唯一性,就有

    将上式代入(8.4.3)即有

    再注意到

    代入上式便得

    由于,故

    是上三角矩阵,加之是对角矩阵,因此,主对角线以下的元素趋于0.

    定理证毕

    8.5.2 Hessenberg

    实际计算时,为了减少每次迭代所需的运算量,总是先将原矩阵A经相似变换约化为一个上准三角矩阵——上Hessenberg矩阵,然后对约化后的矩阵进行QR迭代.

    定义8.4.1 ,当.则称Hessenberg矩阵上准三角矩阵.

    下面我们希望讨论这样一个问题:对,存在上Hessenberg矩阵与之相似.也就是说,总存在一个非奇异矩阵使得

    对于给定的,我们经次正交变换,就可将其约化为上Hessenberg矩阵.

    首先,我们将矩阵分块写成

    其中

    第一步:构造Householder变换,使其满足

    其中阶单位矩阵的第一列.再令

    (8.4.10)

    于是

    显然

    也就是说经第一步正交变换后,我们得到与相似的矩阵形如

    第二步:先对,进行类似于第一步的正交变换,将写成

    构造Householder变换,使得阶单位矩阵的第一列.令

    使得

    显然

    从而令

    于是

    如此进行步,就可找到Householder变换,使得

    其中

    再令

    则有

    (8.4.11)

    通常称该式为A的上Hessenberg分解。

    注意:根据以上Hessenberg约化过程可以发现,正交矩阵的第一列为,即

    算法8.4.1(计算上Hessenberg分解:Householder变换)

    这一算法计算出的上Hessenberg矩阵就存放在A所对应的存储单元内.运算量为;如果需要累积,则还需要再增加运算量

    此外,上述算法计算得到的上Hessenberg矩阵满足

    其中是正交矩阵,,这里是一常数,是机器精度.

    当然,我们亦可用Givens变换将Hessenberg化,一般所需要的运算量大约是算法8.4.1的二倍.但是,如果有较多的零元素,则适当安排Givens变换的次序,可使运算量大为减少.另外,为了节省运算量,也可采用列主元的Gauss消去法将约化为上Hessenberg矩阵.不过这样做,虽然运算量少,但数值稳定性差.

    尽管一般来讲上Hessenberg分解是不唯一的,然而我们可以证明

    定理8.4.3 设有如下两个上Hessenberg分解:

    8.4.12

    其中阶正交矩阵,是上Hessenberg矩阵.若,而且的次对角元素均不为零,则存在对角元素均为1或-1的对角矩阵使得

    (8.4.13)

    [证明] 假定对某个已证

    (8.4.14)

    其中.下面来证明存在使得

    从(8.4.12)可得

    分别比较上面两个矩阵等式的第列,可得

    8.4.15

    8.4.16

    分别在(8.4.15)(8.4.16)两边左乘,可得

    再利用(8.4.16)就有

    (8.4.17)

    (8.4.17)代入(8.4.15),并利用 (6.4.16) (8.4.18),可得

    (8.4.18)

    注意到,则

    ,故(8.4.18)蕴含着

    其中

    定理证毕

    定义8.4.2 一个上Hessenberg矩阵,如果其次对角元素均不为零,即,则称它是不可约的

    上述定理表明:如果这不可约的上Hessenberg矩阵,其中为正交矩阵,则完全由的第一列确定(为里是在相差一个正负号的意义下的唯一).

    现在假定是上Hessenberg矩阵,我们来考虑对进行一次QR迭代的具体实现问题.第一步是计算的QR分解.由于的特殊性,这一步可用个平面旋转变换来完成,其具体计算过程如下:

    第1次约化:作Givens变换

    其中仅有前两行的元素与的元素不同,其余各行元素与均相同,并且

    , 

    直至第k-1次约化毕,我们得到

    k次约化:令

    Givens变换

    那么

    其中仅第k行和第k+1的元素与的元素不同,其余各行元素与均相同,并且

    , 

    这样,我们便知

    是上三角矩阵.令,则,即这样就已完成了的QR分解.要完成一次QR迭代,我们还须计算

    由于是(1,2)坐标平面内的旋转变换,因此仅有前两列与不同,而的前两列由的前两列的线性组合构成.又是上三角矩阵,故必有如下形状:

    同样,是(2,3)坐标平面内的旋转变换,仅有第二列和第三列与不同,它们是的和第二列和第三列的线性组合,故有如下形状:

    如此进行下去,最后我们得到的仍是上Hessenberg矩阵.而且不难算出,这样进行的一次QR迭代的运算量是.注意,对一般方阵进行一次QR迭代的运算量是

    8.5.3 带原点位移的QR迭代

    从定理8.4.2已经知道,基本的QR算法是线性收敛的,其收敛速度取决于特征值之间的分离程度。为了加速其收敛速度,类似于反幂法,可引进原点位移。设第步迭代的位移为,则带原点位移的QR迭代如下:

    这里是给定的上Hessenberg矩阵。

    现在来讨论位移的选取。由于为上Hessenberg矩阵,故其最后一行仅有两个非零元素,若QR算法收敛,则当充分大时,就很小,因而就接近于的一个特征值。这样,根据从反幂法所得到的经验,我们可以选取位移为。事实上,对于这样选取的位移,可以证明,若很小的话,则经一次带原点位移QR迭代后,成立

    8.4.19

    显然,这只需考虑右下角的子矩阵的变化即可。由前面的讨论可知,将约化成上Hessenberg阵有步,现假定前面步已完成,此时的变为

    这是因为前步不改变的最后一行。约化的第步就是要消去,即确定使得

    从平面旋转变换的确定方法易知,此处

    这样,通过简单的计算可知

    (6.4.24)成立.我们看到通过原点位移,特征值的渐近收敛速度从线性收敛加速面变成了二次收敛.

    8.5.4 双重步位移的QR迭代

    上面所讨论的带原点位移的QR迭代,存在严重的缺点:若具有复共轭特征值,则实位移一般并不能起到加速的作用.为了克服这一缺点,我们下面来介绍双重步位移的QR迭代,其基本思想是将带原点位移的QR迭代合并为一步,以避免复数运算.

    ,考虑如下的迭代:

    8.4.20

    不失一般性,我们可以假定迭代(8.4.20)中出现的上Hessenberg矩阵都是不可约的.因若不然,迭代的某一步,已有

    则我们可以分别对进行QR迭代即可.

    大家已经知道,在一定条件下取位移可起到加速收敛的作用;然而,大家也熟知实矩阵可以有复特征值.这样,假如的尾部子矩阵

    有一对复共轭特征值时,我们就不能期望最终收敛于A的某个特征值,因而此种情形再取位移为就完全起不到加速收敛的作用.为了加速收敛,此时我们自然应该取作位移.但这样一来就必须涉及复数运算,而这又是我们所不希望的.为了避免复运算的出现,人们想用连续作两次位移,即进行

    这里我们记

    但是,这里的都是复数,要避免复数运算肯定不能按照上面给出的迭代直接进行,它只能是一个算法思想。

    非常巧妙的是,上面的算法可以通过实运算得到.为了证实这一点,我们令

    (8.4.21)

    , (8.4.22)

    另外

    即有

    (8.4.23)

    , (8.4.24)

    另一方面,由(8.4.23)可得

    (8.4.25)

    其中

    因此是一个实矩阵;而且如果均不是的特征值,并假定在迭代过程中选取的对角元素均为正数,则由(8.4.23)可推知,Q亦是实的;从而由(8.4.24)也是实的.这也就是说,在没有误差的情况下,用连续作两次位移进行QR迭代产生的仍是上Hessenberg矩阵.但是,实际计算时,由于舍入误差的影响,如此得到的一般并不一定是实的.这样为了确保计算得到的仍是实的,根据(8.4.23)(8.4.25),我们自然想到按如下的步骤来计算

    计算QR分解:

    计算

    然而,如此计算的第一步形成的运算量就是.当然,这是我们不希望的.幸运的是,定理8.4.3告诉我们:不论采用什么样的方法去求正交矩阵使是上Hessenberg矩阵,只要保证的第一列与的第一列一样,则就与本质上是一样的(所有元素的绝对值都相等).当然,这需要是不可约来加以保证的.因此是不可约的,则我们就可有很大的自由度去寻求更有效的方法来实现由的变换.下面的定理给出不可约的条件.

    定理8.4.4 是不可约的上Hessenberg矩阵,且均非的特征值,则也是不可约上Hessenberg矩阵.

    [证明]用反证法.记,并假定存在使得,而.比较等式的两边矩阵的前列,得

    由此可得

    (8.4.26)

    其中

    将代入(8.4.26),并注意到也是的多项式,就有

    (8.4.27)

    其中

    ,并注意到是不可约的上Hessenberg矩阵,直接计算可知的第个分量为

    这也就是说(8.4.27)有非零解,而这与均非的特征值蕴含着非奇异矛盾.

    定理证毕

    基于定理8.4.26.4.4,我们可以从另外的途径来实现的变换.首先,我们从(8.4.21)知,的第一列与的第一列共线(其实的第一列就相当于由的第一列单位化而得到的).而由(8.4.25)容易算出

    其中

    其次,如果Householder变换使:,从而有,即的第一列就与共线,从而的第一列就可作为的第一列,即,而由关于Householder变换的理论知,可以按如下方式确定:

    其中

    现令

    则我们只要能够找到第一列为的正交矩阵使为上Hessenberg矩阵,那么就是我们希望得到的.由本节所介绍的约化一个矩阵为上Hessenberg矩阵,的方法可知,这是容易办到的.这只需确定Householder变换使

    为上Hessenberg矩阵,则的第一列就为.而且由于所具有的特殊性,实现这一约化过程所需的运算量仅为

    事实上,由于用相似变换为只改变了的前三行和前三列,故有如下形状

    仅比上Hessenberg形多三个可能的非零元素"+".由的这种特殊性易知,用来约化为上Hessenberg形的第一个Householder变换具有如下形状

    其中为3阶Householder变换,而且具有如下形状

    一般地,第k次约化所用的Householder变换具有如下形状

    其中为3阶Householder变换,,而且具有如下形状

    因此,最后一次约化所用的Householder变换具有如下形状

    其中为2阶Householder变换.

    综述上面的讨论,就得到了著名的Francis双重步位移的QR迭代算法:

    算法8.4.2(双重步位移QR迭代)

    m=n-1

    s=H(m,m)+H(n,n)

    t=H(m,m)H(n,n)-H(m,n)H(n,m)

    x=H(1,1)H(1,1)+H(1,2)H(2,1)-sH(1,1)+t

    y=H(2,1)(H(1,1)+H(2,2)-s)

    z=H(2,1)H(3,2)

    for k=0:n-3

    q=max{1,k}

    H(k+1:k+3,q:n)= H(k+1:k+3,q:n)

    R=min{k+4,n}

    H(1:r,k+1:k+3)=H(1:r,k+1:k+3)

    x=H(k+2,k+1)

    y=H(k+3,k+1)

    if k<n-3

    z=H(k+4,k+1)

    end

    end

    H(n-1:n,n-2:n)= H(n-1:n,n-2:n)

    H(1:n,n-1:n)= H(1:n,n-1:n)

    该算法的运算量为;如果需要累积正交变换投影还需要增加运算量

    8.5.5 隐式QR算法

    前面的讨论已经解决了用QR方法求一个给定的实矩阵的Schur标准形的几个关键的问题.然而,作为一种实用的算法,还需给出一种有效的判定准则,来判定迭代过程中所产生的上Hessenberg矩阵的次对角元何时可以忽略不计.一种简单而实用的准则是,当

    时,就将看作零.这样做的理由是,在前面约化为上Hessenberg矩阵时就已经引进了量级为的误差.

    综合上面的讨论,就得到如下的隐式QR算法.该算法是计算一个给定的阶实矩阵的实Schur分解:,其中是正交矩阵,为拟上三角矩阵,即对角块为方阵的块上三角矩阵,而且每个的对角块必为有一对复共轭特征值.

    算法8.4.3(计算实矩阵的实Schur分解:隐式QR算法)

    输入

    Hessenberg化:用算法8.4.1计算的上Hessenberg分解,得

    收敛性判定:

    (i)把所有满足条件

    置零;

    (ii)确定最大的非负整数和最小的非负整数,使

    其中为拟上三角形,而为不可约的上Hessenberg形;

    iii)如果,则输出有关信息,结束;否则进行下一步.

    QR迭代:对用算法8.4.2迭代一次得

    计算

    然后转步(3).

    实际计算的统计表明,这一算法每分离出一个子矩阵平均约需2次QR迭代.因此,如果只计算特征值,则运算量平均约为;如果都需要,则运算量平均约为

    关闭提示 关闭

    确 认 取 消
    展开全文
  • 二、矩阵二阶偏导数  1. X轴方向上二阶偏导  2. Y轴方向上二阶偏导  3.第一次在X轴方向上求偏导,第二次在Y轴方向上求偏导  4.第一次在Y轴方向上求偏导,第二次在X轴方向上求偏导  5.结论 三、拉普拉斯...
  • “数据和特征决定了机器学习的上限,而模型和算法只是逼近这个上限而已”。线性模型是在统计机器学习中常用的模型,我们假设解释变量和响应变量的关系是线性的。真实情况未必如此。如果想仿造一段曲线,那么首先应该...
  • 推荐系统CTR预估学习路线:从LR到FM/FFM探索二阶特征的高效实现 推荐系统CTR预估学习路线:利用树模型自动化特征工程 推荐系统CTR预估学习路线:深度模型 推荐系统CTR预估学习路线:引入注意力机制作者:Eric陈健锋 ...
  • 利用QR算法求解矩阵特征值和特征向量

    万次阅读 多人点赞 2015-03-05 22:02:55
    利用QR算法求解矩阵特征值和特征向量为了求解一般矩阵(不是那种幼稚到shi的2 x 2矩阵)的特征值. 根据定义的话,很可能需要求解高阶方程... 这明显是个坑...高阶方程你肿么破...折腾了好久 1.我要求...
  • I'm trying to take a second derivative in python with two numpy arrays of data.For example, the arrays in question look like this:import numpy as npx = np.array([ 120. , 121.5, 122....
  • 在数列、微分方程、矩阵三个不同的领域中都见到了“特征值”这个名字,难道仅仅是因为他们都刻画了问题的特征吗?当时我想,一定不是的。或许这些特征值是同一个特征值。因为本人水平较低,这个...
  • 从机器学习、量子计算、物理到许多数学和工程的问题,都可以通过找到一个矩阵特征值和特征向量来解决。根据定义(标量λ、向量v是特征值、特征向量A):视觉上,Av与特征向量v位于同一直线上。这里有些例子。然而,Ax...
  • C++——多项式拟合

    千次阅读 多人点赞 2018-07-20 15:37:32
    C++——多项式拟合 目标:利用C++对txt或者xml中的数据,进行高阶或低阶多项式拟合  为方便以后查找,代码以及详细资料已打包,并上传至云盘(链接:https://pan.baidu.com/s/1bvUBIoxv7Avxeq_Cz6xOZQ 密码:u9qe...
  • 矩阵特征向量和特征值1. 特征值与特征向量定义2. 特征值的几何意义3. eigshow4. 绘图 1. 特征值与特征向量定义 设A是n阶方阵,如果存在常数λ和n维非零列向量x,使得等式Ax=λx,则称λ为A的特征值,x是对应特征值...
  • 多项式回归

    2018-08-22 21:50:00
    机器学习中一种常见的模式,是使用线性模型训练数据的非线性函数。这种方法保持了一般快速的线性方法的性能,同时允许它们适应更广泛的...如果我们想把抛物面拟合成数据而不是平面,我们可以结合二阶多项式...
  • /*********************************************************************************  该程序实现对n维输入x,n维输出y,给出n阶多项式拟合系数,功能和matlab
  • 多项式插值

    千次阅读 2015-10-09 15:12:25
    给定n+1n+1个点{(xi,yi)}ni=0\{(x_i,y_i)\}^{n}_{i=0}(称为插值点),所谓多项式插值就是找到一个多项式(称为插值多项式) y=P(x)=akxk+ak−1xk−1+⋯+a1x+a0y=P(x)=a_kx^k+a_{k-1}x^{k-1}+\cdots+a_1x+a_0 使得...
  • 多项式拟合正弦函数

    2018-10-14 00:11:54
    实验要求: 1. 生成数据,加入噪声; 2. 用高阶多项式函数...求解解析解时可以利用现成的矩阵求逆。梯度下降,共轭梯度要求自己求梯度,迭代优化自己写。不许用现成的平台,例如pytorch,tensorflow的自动微分工具。
  • 全局多项式插值法可以根据整个表面拟合多项式,而局部多项式插值法可以对位于指定重叠邻域内的多个多项式进行拟合。通过使用大小和形状、邻域数量和部分配置,可以对搜索邻域进行定义。...二阶全局多项式
  • 本篇博客针对三种联系十分紧密的矩阵分解(Schur分解、特征值分解、奇异值分解)依次介绍,它们的关系是Schur→EVD→SVDSchur\rightarrow{}EVD\rightarrow{}SVDSchur→EVD→SVD,也就是说由Schur分解可以推导出EVD,...
  • 矩阵A无逆矩阵,则称A为奇异矩阵。若A有逆矩阵,则称A是非奇异矩阵,简称非异阵 N阶行列式的超平行多面体的几何图形是由行(或列)向量张成的,而且这个n维超平行多面体与一个n维超长方体等体积。 一、代数意义 ...
  • 常见的曲线拟合方法   1.使偏差绝对值之和最小    2....   3....按偏差平方和最小的原则选取拟合曲线,并且采取二项式方程为拟合曲线的方法,称为...多项式拟合介绍 多项式拟合公式 多项式拟合问题描述 假定
  • 特征提取(二)Hessian矩阵

    千次阅读 2020-03-15 21:21:32
    黑塞矩阵(Hessian Matrix),是一个多元函数的二阶偏导数构成的方阵,描述了函数的局部曲率。黑塞矩阵常用于牛顿法解决优化问题,利用黑塞矩阵可判定多元函数的极值问题。在工程实际问题的优化设计中,所列的目标...
  • 将U PMNS矩阵视为SU(3)组的元素,并为其构造二阶矩阵多项式。 解决了构造混合矩阵的对数的反问题。 这样,标准参数化与指数参数化正好相关。 指数形式允许轻松分解和分别分析旋转和CP违规。 利用有关中微子混合的...

空空如也

空空如也

1 2 3 4 5 ... 20
收藏数 4,243
精华内容 1,697
关键字:

二阶矩阵的特征多项式