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

1. clear;clc;close all
2. tic
3. %并行设置%
4. delete(gcp('nocreate'));%关闭之前的并行
5. c=parcluster('local');
6. c.NumWorkers=50;%手动设置线程数(否则默认最大线程为12)
7. parpool(c, c.NumWorkers);
8. %并行设置%
10. w=1;
11. u=2.5;
13. N=30;
16. valsto=zeros(2*N+2,100);
17. v10=linspace(-3,3,100).';
20. parfor jj=1:length(v10)
23.     v1=v10(jj);
25.     [~,val]=eig(vpa(H(u,v1,w,N),30));%设置30位高精度计算
26.     %[~,val]=eig(H(u,v1,w,N));
27.     tval=diag(val);
28.     [~,index]=sort(real(tval));
29.     val=tval(index);
30.     valsto(:,jj)=val;
32. end
34. figure
35. hold on
36. plot(linspace(-3,3,100),abs(valsto(:,:)),'k.')
38. %%%%%%%%%%=========%%%%%%%%%%%%%%
39. %%%%%%%%%%=========%%%%%%%%%%%%%%
43. clearvars -except N
44. syms u v w real
47. val2=eig(H(u,v,w,N));
48. val2=subs(val2,{u,w},{2.5,1});
50. valsto2=zeros(2*N+2,100);
51. v10=linspace(-3,3,100).';
54. parfor jj=1:length(v10)
57.     v1=v10(jj);
58.     tval=val2;
59.     tval=double(subs(tval,v,v1));
60.     [~,index]=sort(real(tval));
61.     tval=tval(index);
62.     valsto2(:,jj)=tval;
64. end
66. plot(linspace(-3,3,100),abs(valsto2(:,:)),'bo')
67. hold off
68. saveas(gcf,"VAPtest0.fig")
69. toc
71. function [H0]=H(u,v,w,N)
73. T1=[];
74. T2=[];
75. for hh=1:N
76.     T1=[T1;v+u;w+u];
77.     T2=[T2;v-u;w-u];
78. end
79. T1=[T1;v+u];
80. T2=[T2;v-u];
81. H0=diag(T1,1)+diag(T2,-1);
84. end

复制代码

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

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

免责声明：如果侵犯了您的权益，请联系站长，我们会及时删除侵权内容，谢谢合作！更多信息从访问主页：qidao123.com:ToB企服之家，中国第一个企服评测及商务社交产业平台。