MATLAB特征值的计算之eig()函数存在的问题

      本人这段时间一直在研究特征值计算问题，当然从理论上来说，这个问题很简单。甚至，我们自己可以通过公式，用C语言和FORTRAN语言编制相应的代码来实现，或者用PYTHON也行。结果，我选择了编程效率很高的MATLAB,直接调用了里面的eig()函数，发现了比较大的问题。

弹性系数矩阵可以通过voigit notation后，利用christoffel equation得到质点偏振 polarization vector. 通常对于各向同性二维介质来说，christoffel矩阵的形式很简单，元素值为Gij, i,j=1,2:

       kristoffel矩阵{Gij}的特征值就是上面矩阵对角线上要减去的值\rho V^2, V是相速度，特征向量为质点偏振方向。对于各向同性介质而言，知道lamda和mu可以算出其偏振方向。因此，编制MATLAB脚本，准备运算。lambda=6,mu=3;二维christoffel矩阵各元素为Gij：

       G11=C11*kx^2+C44*kz^2;
       G12=(C13+C44)*kx*kz;
       G22=C44*kx^2+C33*kz^2;

         二维各向同性介质弹性系数矩阵元素为：C11=C33=lamda+2*mu;   C12=lamda;  C44=mu ;

>>[vec,lam]=eig([12*sin(1.7*3.14)^2+3*cos(1.7*3.14)^2 9*sin(1.7*3.14)*cos(1.7*3.14);9*sin(1.7*3.14)*cos(1.7*3.14) 12*cos(1.7*3.14)^2+3*sin(1.7*3.14)^2])

vec =

   -0.5856   -0.8106
   -0.8106    0.5856                                                                                                      (1)

后面验证，该结果与理论计算不符合。

但是用符号函数计算，过程如下：

>>syms t;

>> [vec,lam]=eig([12*sin(t)^2+3*cos(t)^2 9*sin(t)*cos(t);9*sin(t)*cos(t) 12*cos(t)^2+3*sin(t)^2])
 
vec =
 
[ -cos(t)/sin(t), sin(t)/cos(t)]
[              1,             1]

明显看到符号函数计算出来的P波polarization vector为[sint(t),cos(t)]；

调用

>>m1y=matlabFunction(vec);

>>m1y(1.7*pi)

ans =

    0.7265   -1.3764
    1.0000    1.0000                                                                                        (2)

对比(1)和（2）其结果显然不同，(2)中第二行是经过了归一化。

（1）和（2）计算结果和极化方向不一致，第四象限和第二象限相同，也就是第四象限刚好与理论结果相反。

需要特别注意。实际应用中将sin(t)替换成nx(kx),cos(t)替换成nz(kz).

>> [vec,lam]=eig([12*kx^2+3*kz^2 9*kx*kz;9*kx*kz 12*kz^2+3*kx^2])

特征向量：

vec =
[ -kz/kx, kx/kz]
[      1,     1]                                                                                                                         （3）

特征值：

lam =
[ 3*kx^2 + 3*kz^2,                 0]
[               0, 12*kx^2 + 12*kz^2]                                               

所以对(3)进行必要的修改，具体根据象限来修改。

vec=[-kz kx;kx kz];这个结果应该对了。

利用PYTHON编程，对上述计算实现后，可以对比结果，作为另外一种参考。

import numpy as np
def frange(start,stop, step=1.0):
       while start < stop:
            yield start
            start +=step
for x in frange(1.7*3.14,1.7*3.14+0.1,0.1):
          G11=12*np.sin(x)**2+3*np.cos(x)**2
          G12=9*np.sin(x)*np.cos(x)
          G22=12*np.cos(x)**2+3*np.sin(x)**2
          G=np.array([[G11,G12],[G12,G22]])
          t=np.linalg.eig(G)
          print(t)

该脚本运行的结果为：                            （3）

(array([12.,  3.]), array([[ 0.81060546,  0.58559268],
       [-0.58559268,  0.81060546]]))                        

显然3和12是特征值，后面是特征向量。显然（3）的运算结果与MATLAB结果相比较而言，它们的特征向量的方向不同。不过python自己的第二象限和第四象限的特征向量也是等同的。这与极化方向理论不符，需要修改。

同理，我们还可以使用MAPLE的符号计算功能，求christoffel矩阵的特征值和特征向量，

先要加载线性代数库：

with(linalg)

接着

A := matrix(2, 2, [12*kx^2 + 3*kz^2, 9*kx*kz, 9*kx*kz, 3*kx^2 + 12*kz^2]);
              

 

 

 MAPLE和MATLAB的结果一致，说明MATLAB调用的是MAPLE里面的模块。

上面弹性参数代表的各向同性介质，可以看到P波偏振向量（特征向量）就是(kx,kz)；S波的偏振向量就是(-kz,kx)。然而对于各向异性介质，如二维VTI介质来说，稍微有些复杂，

上式为弹性参数矩阵得到的kristoffel矩阵A，特征值和特征向量

 

MATLAB 代码:

>> syms C11 C44 C12 C33 kx kz

>> cg0 = [C11 * kx ^ 2 + C44 * kz ^ 2 (C12 + C44) * kx * kz; (C12 + C44) * kx * kz C33 * kz ^ 2 + C44 * kx ^ 2];%MAPLE公式可以转换成MATLAB代码。

>> [vec,lam]=eig(cg0)
 
vec =
 
[ ((C11*kx^2)/2 + (C44*kx^2)/2 + (C33*kz^2)/2 + (C44*kz^2)/2 - (C11^2*kx^4 - 2*C11*C33*kx^2*kz^2 - 2*C11*C44*kx^4 + 2*C11*C44*kx^2*kz^2 + 4*C12^2*kx^2*kz^2 + 8*C12*C44*kx^2*kz^2 + C33^2*kz^4 + 2*C33*C44*kx^2*kz^2 - 2*C33*C44*kz^4 + C44^2*kx^4 + 2*C44^2*kx^2*kz^2 + C44^2*kz^4)^(1/2)/2)/(kx*kz*(C12 + C44)) - (C44*kx^2 + C33*kz^2)/(kx*kz*(C12 + C44)),

((C11*kx^2)/2 + (C44*kx^2)/2 + (C33*kz^2)/2 + (C44*kz^2)/2 + (C11^2*kx^4 - 2*C11*C33*kx^2*kz^2 - 2*C11*C44*kx^4 + 2*C11*C44*kx^2*kz^2 + 4*C12^2*kx^2*kz^2 + 8*C12*C44*kx^2*kz^2 + C33^2*kz^4 + 2*C33*C44*kx^2*kz^2 - 2*C33*C44*kz^4 + C44^2*kx^4 + 2*C44^2*kx^2*kz^2 + C44^2*kz^4)^(1/2)/2)/(kx*kz*(C12 + C44)) - (C44*kx^2 + C33*kz^2)/(kx*kz*(C12 + C44))]
[1, 1]
 

最终发现，无论各向同性介质还是各向异性介质，用matlab的eig()函数得到的是已经归一化的特征向量。例如本来向量是(kx,kz)，而matlab得到的是(kx/sqrt(kx^2+kz^2,kz/sqrt(kx^2+kz^2)))。因此，对于各向同性介质得到的特征向量（ux，uz）需要乘以系数sqrt(kx^2+kz^2)。