1 2 3 4 5 6 7 8 9 10 11 12 13 14 15 16 17 18 19 20 21 22 23 24 25 26 27 28 29 30 31 32 33 34 35 36 37 38 39 40 41 42 43 44 45 46 47 48 49 50 51 52 53 54 55 56 57 58 59 60 61 62 63 64 65 66 67 68 69 70 71 72 73 74 75 76 77 78 79 80 81 82 83 84 85 86 87 88 89 90 91 92 93 94 95 96 | %% main.m X=0.1; C=(16/X)+1; d1=zeros(7,C); e1=zeros(7,C); MM = repmat(1:7, 7, 1); for k=1:C for p = 1:7 eval(sprintf('x%i = 0;',p)); end K=0.1*k-8.1; A = arrayfun(@(i,j) f(i, j, K, x1, x2, x3, x4, x5, x6, x7), MM', MM); [ve, ei]=eig(A); S=zeros(1,7); for w=1:1:7; for v=1:1:7; for u=1:1:7; S=S+A(v,u)* ve(u,w) *conj(A(v,u)* ve(u,w)); end end end w=1:1:7; S(w)=S; sorted = sort(S); m1=sorted(1); index1=find( S <= m1); for l = 1:7 eval(sprintf('alpha%i = ei(%i, index1);',l,l)); eval(sprintf('x%i = ve(%i, index1);',l,l)); end alpha=alpha1+alpha2+alpha3+alpha4+alpha5+alpha6+alpha7; s1=real(alpha); t1=-imag(alpha); d1(:,k)=s1; e1(:,k)=t1'; end %% f.m function value = f(i, j, K, x1, x2, x3, x4, x5, x6, x7) % disp([i, j]); if i == j dd = i - 4; eval([sprintf('value = (0.5*(K+%i*pi)^2+0.5)', dd), ... '+(1-0.3i)*(2*x1*conj(x1)+2*x2*conj(x2)+2*x3*conj(x3)+2*x4*conj(x4)+x5*conj(x5)+2*x6*conj(x6)+2*x7*conj(x7));']); elseif i > j ii = 7; d = i - j; if d == 6 value = 0; elseif d < 6 && d > 0 value = 0; while 1 if ii - d <= 0 break; elseif ii == i && ii - 1 - d <= 0 break; elseif ii == i ii = ii - 1; end % disp(sprintf('value = value + x%i*conj(x%i);', ii, ii - d)); eval(sprintf('value = value + x%i*conj(x%i);', ii, ii - d)); ii = ii - 1; end value = (1-0.3i) * value; end if d == 1 value = -0.25 + value; end elseif i < j ii = 7; d = j - i; if d == 6 value = 0; elseif d < 6 && d > 0 value = 0; while 1 if ii - d <= 0 break; elseif ii == j && ii - 1 - d <= 0 break; elseif ii == j ii = ii - 1; end % disp(sprintf('value = value + x%i*conj(x%i);', ii - d, ii)); eval(sprintf('value = value + x%i*conj(x%i);', ii - d, ii)); ii = ii - 1; end value = (1-0.3i) * value; end if d == 1 value = -0.25 + value; end end % MM = repmat(1:7, 7, 1); % A = arrayfun(@(i,j) f(i, j, K, x1, x2, x3, x4, x5, x6, x7), MM', MM); |
Direct link: https://paste.plurk.com/show/1822816