function [opt, fall] = htopy(fu, svars, prex0, upperS, numchecks, box, Hchoice, ODEchoice) % htopy2 computes an approximate solution to the system of equations fu(x)=0. % It uses an homotopy represented by H(x,t)=0, with known solution for x at t=0, and % with H(x,1)=fu(x)=0, so any x satisfying (approx) H(x,1) is an (approx) solution to % fu(x), the system of interest. % This implementation continously deforms H(x,0) into a sequence of H(x,t) by making % x and t functions of another variable 's' (thus, x becomes x(s) and t becomes t(s)), % and solves a sequence of systems of ODE in order to 'track' the path defined by H(x,t). % For more explanations about this procedure and formulas involved, see Judd 1998, pages % 179 to 187) % % OUTPUT: % % fall: fall = 1 only if the function succeeds. Success means: 1) a candidate solution was found % (i.e. a t close enough to 1 is found), AND 2) the candidate solution lies inside the % bounds for x specified by the user (in 'box'). In any other case the function fails % and fall = 0. This is useful to instruct other functions not to use the value in opt % if fall is not 1. % % opt: if the function succeeds, opt is a vector of values x % that makes fu(x) approx. zero. If fails because the candidate solution lies % outside the box, then it returns and empty array []. If it fails because no % candidate solution was found it ret)urns the values of x and t for the last % iteration. % % INPUTS: % % fu: a cell array of strings representing a system of equations % (it can be just a single equation) in the form f(x)=0 % % svars : a cell array of strings containig the names of the variables included % in fu (in particular, do not include t or s here) % % prex0 : a vector with the initial conditions for variables in svars % % upperS : this scalar is the upper bound used in the ODE part for s. You can set upperS % relatively large and insert many check points between 0 and upperS. The function % will exit as soon as one candidate solution has been found in any of the % intervals defined by the checkpoints [0,s1,s2....upperS]. That can happen well % before s is close to upperS. % % numchecks: A scalar. Is the number of intermediate checks s1,s2,s3... between 0 % and upperS. The fucntion will divide [0,upperS] in (numchecks+1) intervals of equal % length % % box: an nx2 matrix. Each rows represent one variable in svars. The first column % holds the lower bound for that variable and the second column contains the upper bound. % % Hchoice: a string. If it takes the value 'FP' then the (1-t)(x-x0) is % done. If it takes the value 'N', it implements a Newton type homotopy (see page 179) % % ODEchoice: a string. If it takes the value 'ode45' then the solver ode45 % (runge-kutta, recommended choice for nonstiff problems) is used. If it takes the % value 'ode15s', a BDF-type solver is used. This appropiate for stiff % problem and when precision is less of a concern %%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%% %%% %%% Step 0: set tolerance parameters and default values %%% %%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%% %%% this is our max acceptable deviation from H(x,t) = 0, we make this check %%% at the end of each subinterval [s_{i} , s_{i+1}] Htol = 1e-4; %%% we accept t as close enough to 1 if t \in (1-epsilont , 1+epsilont) epsilont=1e-1; %% set the failure flag to 0, do not change to success unless on sol is found fall=0; %% set a void value for the solution as the default in case of failure opt=[]; %%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%% %%% %%% Step 1: symbolize relevant variables and set the homotopy H %%% %%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%% %%% simbolize s and t and some otehr symbols we may need syms t s pi e a; %%% augment the vector of initial conditions -prex0-, by adding the initial condition for t %%% (zero) at the end. The augmented vector is x0 x0= [prex0, 0]; numvar = size(svars,2); numfun = size(fu,2); %% check your system is square if numfun ~= numvar 'the number of variables is not equal to the number of equations' return; end %% initialize these to be symbolic but empty allvars = sym([]); allfun = sym([]); for i=1:numvar allvars = [allvars, sym(svars{i})]; allfun = [allfun, sym(fu{i})]; end %% now allvars contains all variables in symbolic form %% and allfun has all equations in fu but with symbolic status %%% Define your Homotopy if Hchoice=='FP' % this is the fixed point homotopy H = (1-t)*(allvars-prex0) + t*allfun elseif Hchoice=='N' % Newton's homotopy f0=allfun; allfun x0 for a=1:size(allvars,2) f0=subs(f0,{allvars(a)},{x0(a)}); end f0 H=allfun-(1-t)*f0 else display('Error: no Homotopy was selected or you typed an incorrect name') return; end %%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%% %%% %%% Step 2: Construct the system of ODEs (see eq. 5.9.4 , page 182) %%% %%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%% %%%%%% Step 2.1: obtaing the RHS of the system %% append symbolic t to end of list of symbolic x allvars = [allvars, t]; %% compute the symbolic Jacobian of H w.r.t. x and t J=jacobian(H,allvars); %% form the Right Hand Side of the ODE system allRHS = []; for i = 1:numvar+1 subJ = J; %% exclude the i-th column of the Jacobian before taking Determinant subJ(:,i)=[]; %% concatenate each eq vertically, do the sign switch allRHS=[allRHS; (-1)^i * (det(subJ))]; end %% display the RHS of the ODE system if mod(numvar,2)==0 allRHS=(-1)*allRHS; end %% now we change the notation of the symbolic expressions that appear %% in allRHS. The reason is to get an expression that, once transformed %% into a string (see step 2.2) will be usable as an inline function %% accepted as input by fsolve (see step 2.3). %% the first step to change the notation of allRHS is to create a row %% vector y=[y(1), y(2), ..., y(numvar+1)], where numvar was defined in %% step zero and corresponds to allvars-1. yy=sym([]); for j=1:size(allvars,2) sj = num2str(j); name = strcat('y(',sj,')'); yy(j)=name; end y=yy; %% the second step in changing the notation is to sobstitute the original %% variables given in allvars with the y(i) ones we constructed above. [m1,n1]=size(allRHS); B=jacobian(allRHS,allvars); for j=1:m1 for i=1:m1 if B(j,i)==0 allRHSs(j)=allRHS(j); else allRHSs(j)=subs(allRHS(j),{allvars(i)},{y(i)}); end allRHS(j)=allRHSs(j); end end allRHS %%%%%%% Step 2.2 %%%%%%% change the symbolic object to string object, so you can use inline functions newstring = '['; for i=1:size(allRHS,1) for j=1:1 newstring = [newstring char(allRHS(i,j)) ' ']; end newstring = [newstring ' ; ']; end newstring = newstring(1:end-3); %remove trailing punctuation newstring = strcat(newstring, ']'); %%%%%%% Step 2.3 : %%%%%%% create an inline function to make RHS usable with Matlab's ODE-functions f=inline(newstring,'s','y'); %%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%% %%% %%% Step 3: Do the path-tracking via ODE-solving %%% %%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%% %% these points define a uniform grid on [0,upperS] checkpoints = linspace(0,1,numchecks+1)*upperS; %% container for solutions passed by the ODE solver pastedvalues = []; pastedtime =[]; %% Htest will be used to asses value of H at the checkpoints. For now %% it is a symbolic object but it will be converted to a number once we %% have computed the values of t and the svars at which we want to compute %% the value of H. Htest = H; k=0; for j=1:numchecks j if (ODEchoice=='ode45') [timepoints, varvalues] = ode45(f,[checkpoints(j),checkpoints(j+1)],x0); elseif (ODEchoice=='ode15') [timepoints, varvalues] = ode15s(f,[checkpoints(j),checkpoints(j+1)],x0); else H=h; display('Error: no ode-solver is being called or you typed an incorrect name') return; end %% store the last chunk of x's t and s pastedvalues = [pastedvalues; varvalues(2:end,:)]; pastedtime = [pastedtime; timepoints(2:end)]; %% update the initial conditions for x and t, to be used in the next iteration x0 = varvalues(end,:); %% find the value of t closest to 1 for this [s_i,s_{i+1}] interval testvec =abs(varvalues(:,end) -1); mintestvec= min(testvec); k = k+size(testvec,1)-1; %% test if our 'best' t is close enough to 1, if it's not, just go on if (abs(mintestvec) < epsilont) %% if t is accepted then report a candidate solution z = min(testvec); i=1; while (testvec(i) ~= z) i = i+1; end %% opt does not include tH=h; %% set opt equal to your candidate solution opt = varvalues(i,1:end-1); %% test if your candidate solution lies inside the box for x or not if (opt >= box(:,2)') | (opt <= box(:,1)') display ('Out of bounds') out_of_bound_guess=opt %% your solution is out of the permisible region %% exit the function, do not report fall=1 return; end %% your solution passed the test, report success fall=1; %% in the case of success the program will create plots for each of %% the original variables. Each plot has t on the vertical axis and %% one of the independent variables on the orizontal axis. if m1 == 3 X1=pastedvalues(1:(k-size(testvec,1)+i),1); Y1=pastedvalues(1:(k-size(testvec,1)+i),2); Z1=pastedvalues(1:(k-size(testvec,1)+i),3); plot3(X1,Y1,Z1) grid on xlabel('x') ylabel('q') zlabel('t') title('3D Plot of the solution') else for l=1:m1-1 h(l)=subplot(m1-1,1,l); plot(pastedvalues(1:(k-size(testvec,1)+i),l)',pastedvalues(1:(k-size(testvec,1)+i),end)') xlabel('independent variable') ylabel('t') axis([h(l)],'tight') end end return end %% A cnadidate solution was not found on this subinterval. Before %% checking the next subinterval (that is, going to the next j) we %% assess if our new guess is such to satisfy the requirement that %% H is equal to zero along all the path. %% Tah is, we check taht H(x0,t)=0, where x0 is the last value of %% varvalues. for a=1:size(allvars,2) Htest=subs(Htest,{allvars(a)},{x0(a)}); end Htest if (abs(Htest)/(1+norm(x0)) > Htol) %% if we are off-path then use a solver to correct x, fixing t display('Warning: off-the path. The Htol loop has been activated to correct the initial guess for the next subinterval'); % Now the program uses fsolve to produce a new vector of initial % guesses for the original independent variables, given of t which % we get at the end of the subinterval. % To produce the new guess we proceed in three parts: %% PART I: H is evaluated at the last value of t and then the %% notation is changed so to get an expression that, once %% transformed into a string (see PART II) will be usable as an %% inline function accepted as input by fsolve (see PART III). [m2,n2]=size(H); y(n2+1)=x0(n2+1); B2=jacobian(H,allvars); H2=H; for u=1:n2 for v=1:m1 if J(u,v)==0 newHs(u)=H2(u); else newHs(u)=subs(H2(u),{allvars(v)},{y(v)}); end H2(u)=newHs(u); end end newH=sym([]); r0=zeros(n2,1); for r=1:n2 newH=[newH;H2(r)]; r0(r)=x0(r); end newH % shows H(y(1),...,y(numvar),t*). Where t* is the last value % of t computed by the ODE solver for the subinterval at hand. % PART II: then we need to make newH usable with fsolve. We make it a % string. newstring2 = '['; for c=1:size(newH,1) for w=1:1 newstring2 = [newstring2 char(newH(c,w)) ' ']; end newstring2 = [newstring2 ' ; ']; end newstring2 = newstring2(1:end-3); %remove trailing punctuation newstring2 = strcat(newstring2, ']'); %% PART III: we use the above string into an inline function %% that becomes the input of fsolve together with an initial %% guess for the original independent variables given by the last %% values generated for them by the ODE solver. %% The new guess we generate through fsolve is called "newx". f2=inline(newstring2,'y') newx=fsolve(f2,r0) %% Then we also implement a check to see how well we are on the path %% with this new values of the independent variables and the value %% of that ends the subinterval. That is we evaluate H at %% (newx,t*). We call "newHtest" the value that we compute and that %% we will compare to Htol. x02=[(newx)',x0(n2+1)] newHtest=H; for a=1:size(allvars,2) newHtest=subs(newHtest,{allvars(a)},{x02(a)}); end newHtest if abs(Htest) Htol) display('Matlab equation solver gave a result such that H is far from zero. Check Htol') return; else x0=x02 end % so either we start the check on the next subinterval with an % initial guess so that H is close enough to zero or we get one % of the two error messages above and the program stops. end %% this complete the loop for a single subinterval end %% all the subintervals in [0, upperS] have been examined %% if you made it up to this point without exit the function, it means that you never found %% a candidate solution display ('Warning: t did not get close enough to 1: increase upperS or refine stepsize for t')