bvp_hom.asv

来自「计算动力学系统的分岔图」· ASV 代码 · 共 59 行

function f = BVP_Hom(x,YS,YU)
global hetds 

% extract ups
ups = reshape(x,hetds.nphase,hetds.npoints);
p = num2cell(p);
% -----------
% equilibrium condition
 J=hetds.niteration;
 x0=x(:,1);
 for i=1:J
 f(1:hetds.nphase,1)=feval(hetds.func,0,x0,p{:}); %f^i(x0)-x0
 end
 idx=hetds.nhase*(hetds.npoint-1);
 % f^J(x-i,\alpha)-x_{i+1}=0, I=2,3,...,npoints-1
 idx=hetds.nphase;
 for i=2:hetds.npoints-2
     xx=x(:,i);
    for i=1:J
     xx=feval(fpmds.func,0,xx,p{:});
    end
    f(idx+1:i*hetds.nphase,1)=x(:,i+1)-xx;%f^i(xi-1)-xi
    idx=idx+hetds.nphase;
 end
x1=x(:,hetds.npoints)
 for i=1:J
 f(idx+1:npoints*nphase,1)=feval(hetds.func,0,x1,p{:}); %f^i(x1)-x1
 end 

% Component 4
% -----------
% Ricatti equations
R(1:hetds.nu,1) = RicattiEval(x0,p,1,YU);
idx=hetds.nu
R(idx+1:hetds.ns,1) = RicattiEval(x1,p,1,YS);
% Component 5
% -----------
% Last and first vectors along stable and unstable eigenspaces
Q0U = hetds.Q0;
Q1S = hetds.Q1;
if hetds.nu
    Q1U = Q0U * [-YU'; eye(size(YU,1))];
    for i=1:hetds.nneg
        f(end+1,1) = (ups(:,1) - x0)' * Q1U(:,end-i+1);
    end
end
if hetds.nneg
    Q1S = Q0S * [-YS'; eye(size(YS,1))];
    for i=1:hetds.npos
        f(end+1,1) = (ups(:,end) - x0)' * Q1S(:,end-i+1);
    end
end

% Component 6
% -----------
% Distances from endpoints to equilibrium equal to epsilons
f(end+1,1) = norm(ups(:,1) - x0) - eps0;
f(end+1,1) = norm(ups(:,end) - x0) - eps1;