【问题标题】:How do I solve the problem the eigenvalue is very slow to be solved analytically?如何解决特征值解析求解很慢的问题?
【发布时间】:2019-09-25 08:51:53
【问题描述】:

我很沮丧这段代码没有运行结果,它一直在运行。我认为原因是解析推导过程需要矩阵 W 的解析解。我希望有人能给我一些帮助或指出问题,我将不胜感激。

tic;
clc;
clear;
syms B
W11=fa(B,1,1,1,1);W12=fa(B,1,1,1,0);W13=fa(B,1,1,0,1);W14=fa(B,1,1,0,0);
W21=fa(B,1,0,1,1);W22=fa(B,1,0,1,0);W23=fa(B,1,0,0,1);W24=fa(B,1,0,0,0);
W31=fa(B,0,1,1,1);W32=fa(B,0,1,1,0);W33=fa(B,0,1,0,1);W34=fa(B,0,1,0,0);
W41=fa(B,0,0,1,1);W42=fa(B,0,0,1,0);W43=fa(B,0,0,0,1);W44=fa(B,0,0,0,0);
W=simplify([W11,W12,W13,W14;W21,W22,W23,W24;W31,W32,W33,W34,W41,W42,W43,W44]);
E1=eig(W);
E2=sort(E1);
E3=E2(4);
T=1;
be=1/T;
f1=-1/be*log(E3);
ma=-diff(f1);
m1=subs(ma,{B},2)
toc;

function s1=fa(B,m1i,m1j,m2i,m2j)
JH=2;J=1;T=1;de=0.2;
sx=1/2*[0 1;1 0];
sy=1/2*[0 -1i;1i 0];
sz=1/2*[1 0;0 -1];
u=eye(2);
be=1/T;
H11=kron(sx,kron(sx,kron(u,u)))+kron(sy,kron(sy,kron(u,u)))...
+de*kron(sz,kron(sz,kron(u,u)));
H12=kron(u,kron(sx,kron(sx,u)))+kron(u,kron(sy,kron(sy,u)))...
+de*kron(u,kron(sz,kron(sz,u)));
H2=kron(sz,kron(u,kron(u,u)))+kron(u,kron(sz,kron(u,u)));
H3=kron(u,kron(u,kron(u,u)));
H=-JH*(H11+H12)+J*H2*(m1i+m1j)+1/2*J*H3*(m1i+m2i+m1j+m2j)...
    -B*1/2*H3*(m1i+m2i+m1j+m2j);
va=eig(H);
s1=sum(exp(-be*va));
end

【问题讨论】:

  • 请说明此行的用途:syms JH J B T de m1i m1j m2i m2j。您可能想首先解释您要编写的表达式是什么,它来自哪里等。您是要为结果找到一个通用表达式,还是要执行一次计算并完成它?如果是后者,我会给你一个友好的提示:不要使用syms。 :)
  • 先谢谢你的建议,但是如果我不使用'syms',如何解决'B'的推导?
  • 从你的帖子看来,using syms 也没有给你答案,所以我不确定这个论点是关于什么的。如前所述,请提供有关您的问题的一些背景信息,以便我们能够看到更大的图景,并可能提出更合适的方法。
  • 仅通过查看代码,我最担心的是 E1=eig(W); 行,您在其中尝试象征性地计算看起来复杂的 W 矩阵的特征值。您是否尝试过调试这段代码,逐行执行,并确定卡在哪里?
  • 是的,你是对的,它卡在 'E1=eig(W)'。

标签: matlab derivative eigenvalue


【解决方案1】:

正如@Dev-iL 在a comment 中所建议的那样,代码的问题在于它试图以符号方式计算非常复杂的表达式。相反,以数字方式评估表达式要简单得多。导数的数值近似是唯一不能得到精确结果的地方,但正如我们稍后会看到的,误差非常非常小。

这就是我重写您的代码的方式。首先,我将您的所有代码添加到 diff 行,并将其制成一个函数(为清楚起见,稍微简化了一些)。这个函数根本不使用符号数学,它只是在B 的给定值下计算f1 的值。请注意,它使用问题中的fa 函数不变。

function out = f1(B)
W = [fa(B,1,1,1,1), fa(B,1,1,1,0), fa(B,1,1,0,1), fa(B,1,1,0,0);...
     fa(B,1,0,1,1), fa(B,1,0,1,0), fa(B,1,0,0,1), fa(B,1,0,0,0);...
     fa(B,0,1,1,1), fa(B,0,1,1,0), fa(B,0,1,0,1), fa(B,0,1,0,0);...
     fa(B,0,0,1,1), fa(B,0,0,1,0), fa(B,0,0,0,1), fa(B,0,0,0,0)];
E = sort(eig(W));
E = E(4);
T = 1;
be = 1/T;
out = -1/be*log(E);
end

接下来,我们通过计算两个非常接近点 B 的位置处的 f1 的值来估算导数(但不要太接近,因为我们会遇到数值舍入误差):

B = 2;
delta = 1e-6;
(f1(B-delta) - f1(B+delta))/(2*delta)

我们改变delta 的值来验证我们有一个好的近似值。使用1e-6 和1e-4 和1e-3 我得到完全相同的值:1.6303。这表明函数在B 附近非常平滑,我们对导数的估计是正确的。

【讨论】:

  • 是的,确实是个好方法。那么二阶导数呢?
  • @KarryMa:你可以通过计算一阶导数两次来得到它。请参阅“二阶中心”下的here。
猜你喜欢
  • 1970-01-01
  • 2014-09-05
  • 1970-01-01
  • 1970-01-01
  • 1970-01-01
  • 2016-02-05
  • 1970-01-01
  • 1970-01-01
相关资源
最近更新 更多