non_trivial.m

利用电磁场的源激发方法来计算光子晶体波导例如光子晶体光纤

% find a non_trivial solution to A*x = 0;
function [err,sol] = non_trivial(mZ, mBNegate, varargin)
if ~isempty(varargin{:})
        opts.v0 = varargin{:};
end
opts.p = 5;
opts.tol = 1e-2;
opts.maxit = 10;
opts.issym=1;
opts.isreal = 0;
opts.disp = 0;

A = mZ'*mZ;
%[rows, cols] = size(A);
%sA = sparse(rows, cols);
%ind = abs(A)>rms(A)/1e4;
%sA(ind) = A(ind);
mZ(mBNegate) = -mZ(mBNegate);
B = mZ'*mZ;
%sB = sparse(rows, cols);
%ind = abs(B)>rms(B)/1e4;
%sB(ind) = B(ind);

%A = Z'*Z;
%[r,c] = size(Z);
%Z1 = [Z(:,1:round(c/2)) -Z(:,round(c/2)+1:end)];
%B = Z1'*Z1;
B = (B+B')/2;
A = (A+A')/2;
try
        %gsvd
        [sol ,err] = eigs(A, B, 1, 'sm', opts);
catch
         sol = 1;
         err = -1;
%         d = eig(A); min(d)
%         d = eig(B); min(d)
%         return
%         disp('Arnoldi algorithm failed - using gsvd...')
%         [r,c] = size(mZ);
%         mZ1 = mZ;
%         mZ(mBNegate) = -mZ(mBNegate);
%         [U,V,X,C,S] = gsvd(mZ,mZ1);
%         [err, ind] =  min(abs(diag(C)./diag(S)));
%         X1 = X'\eye(size(X'));
%         sol = X1(:,ind);
end

err= sqrt(abs(err));
%[u,s,out,flag] = svds(sparse(A),1,0);
%flag
%[rows, cols] = size(A);
%[u,s,v] = svd(A);
%out = v(:,cols-num+1:cols);
%return
%Z = A(1:rows-1,1:cols-1);
%v = -A(1:rows-1,cols);
%I = Z\v;
%out = [I; 1];