matlab中高精度计算函数vpa与非厄米矩阵本征值的求解

clear;clc;close all
tic
%并行设置%
delete(gcp('nocreate'));%关闭之前的并行
c=parcluster('local');
c.NumWorkers=50;%手动设置线程数(否则默认最大线程为12)
parpool(c, c.NumWorkers);
%并行设置%

w=1;
u=2.5;

N=30;


valsto=zeros(2*N+2,100);
v10=linspace(-3,3,100).';


parfor jj=1:length(v10)


    v1=v10(jj);

    [~,val]=eig(vpa(H(u,v1,w,N),30));%设置30位高精度计算
    %[~,val]=eig(H(u,v1,w,N));
    tval=diag(val);
    [~,index]=sort(real(tval));
    val=tval(index);
    valsto(:,jj)=val;

end

figure
hold on
plot(linspace(-3,3,100),abs(valsto(:,:)),'k.')

%%%%%%%%%%=========%%%%%%%%%%%%%%
%%%%%%%%%%=========%%%%%%%%%%%%%%



clearvars -except N
syms u v w real


val2=eig(H(u,v,w,N));
val2=subs(val2,{u,w},{2.5,1});

valsto2=zeros(2*N+2,100);
v10=linspace(-3,3,100).';


parfor jj=1:length(v10)


    v1=v10(jj);
    tval=val2;
    tval=double(subs(tval,v,v1));
    [~,index]=sort(real(tval));
    tval=tval(index);
    valsto2(:,jj)=tval;

end

plot(linspace(-3,3,100),abs(valsto2(:,:)),'bo')
hold off
saveas(gcf,"VAPtest0.fig")
toc

function [H0]=H(u,v,w,N)

T1=[];
T2=[];
for hh=1:N
    T1=[T1;v+u;w+u];
    T2=[T2;v-u;w-u];
end
T1=[T1;v+u];
T2=[T2;v-u];
H0=diag(T1,1)+diag(T2,-1);


end

上面代码没用高精度计算函数vpa之前数值计算的特征值能谱有部分和符号计算(篮圈)是对不上的，如下：

使用了vpa高精度函数后，数值计算和符号计算的结果完全一致：

论文1、DOI: 10.1103/PhysRevLett.125.226402

2、DOI: 10.1103/PhysRevLett.122.250201有讨论非厄米矩阵如何计算出正确特征值的问题。

为了获得一个非厄米矩阵的正确本征值，当非厄米矩阵尺寸变大时，计算精度应该也要随之加大。(Ref：10.1103/PhysRevLett.125.226402)

有一个重要的观念，对于任意大小的非厄米矩阵，能正确计算其本征值所需的最小数值精度是未知的。