How do I solve the problem the eigenvalue is very slow to be solved analytically？

Question

I'm very upset this code doen't run a result,it just running all the time.I think the reason is the analytic derivation process which need a analytic solution of matrix W. I hope someone would give me some help or point the problems,I will be very grateful.

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

Please explain the purpose of this line: `syms JH J B T de m1i m1j m2i m2j`. You might want to start by explaining what is the expression you're trying to code, where it comes from, etc. Are you trying to find a general expression for the result, or do you want to perform the computation once and be done with it? In case of the latter, I'll give you a friendly tip: **DO NOT USE `syms`**. :) — Dev-iL, Sep 25 '19 at 09:18
Thank you for your advice first, but if I don't use 'syms', how can I solve the derivation of 'B'? — Mengr, Sep 25 '19 at 11:39
It appears from your post that _using_ `syms` doesn't give you an answer either, so I'm not sure what the argument is about. As mentioned, please give some background about your problem so that we could see the bigger picture and perhaps suggest more suitable approaches. — Dev-iL, Sep 25 '19 at 12:18
Just by looking at the code, I am most worried about the line `E1=eig(W);` where you're trying to compute the eigenvalues of the complicated-looking `W` matrix symbolically. Did you try debugging this code, stepping line by line, and determining where it gets stuck? — Dev-iL, Sep 25 '19 at 12:24

Cris Luengo · Accepted Answer · 2019-09-25T17:12:11.867

As suggested by @Dev-iL in a comment, the problem with the code is that it tries to compute very complicated expressions symbolically. It is much simpler to evaluate the expression numerically instead. The numerical approximation to the derivative will be the only place where you don't get exact results, but as we'll see later, the error is very, very small.

This is how I rewrote your code. First I took all your code up to the diff line, and made it into a function (simplified slightly for clarity). This function does not use symbolic math at all, it simply computes the value of f1 at the given value of B. Note that it uses the fa function in the question unchanged.

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

Next, we estimate the derivative numerically by computing the value of f1 at two locations really close to the point B (but not too close, as we'll get into numerical rounding errors):

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

We change the value of delta to verify we have a good approximation. With 1e-6 and 1e-4 and 1e-3 I get the exact same value: 1.6303. This indicates that the function is very smooth around B, and our estimate of the derivative is correct.

Yes,it's really a good method. Then what about the second order derivative? — Mengr, Sep 26 '19 at 07:22
@KarryMa: You can get that by computing the first derivative twice. See [here](https://en.wikipedia.org/wiki/Finite_difference#Higher-order_differences) under “second order central”. — Cris Luengo, Sep 26 '19 at 13:07

How do I solve the problem the eigenvalue is very slow to be solved analytically？

1 Answers1