Как мне решить проблему, собственное значение очень медленно решается аналитически analy - PullRequest
Новая
       83

Как мне решить проблему, собственное значение очень медленно решается аналитически analy

0 голосов
/ 25 сентября 2019

Я очень расстроен, этот код не запускает результат, он просто запускается все время. Я думаю, что причина в аналитическом процессе деривации, который требует аналитического решения матрицы 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

1 Ответ

1 голос
/ 25 сентября 2019

Как предложено @ Dev-iL в комментарии , проблема с кодом заключается в том, что он пытается вычислять очень сложные выражения символически.Вместо этого гораздо проще оценить выражение численно.Числовое приближение к производной будет единственным местом, где вы не получите точных результатов, но, как мы увидим позже, ошибка очень и очень мала.

Вот как я переписал ваш код.Сначала я взял весь ваш код до строки diff и превратил его в функцию (немного упрощенную для ясности).Эта функция вообще не использует символическую математику, она просто вычисляет значение f1 при заданном значении B.Обратите внимание, что он использует функцию 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

Далее мы оцениваем производную численно, вычисляя значение f1 в двух точках, действительно близких к точке B (но не слишком близко, так как мы попадем в ошибки округления чисел): 

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

Мы изменим значение delta, чтобы убедиться, что у нас хорошее приближение.С 1e-6 и 1e-4 и 1e-3 я получаю точно такое же значение: 1.6303.Это указывает на то, что функция очень гладкая в районе B, и наша оценка производной верна.

...