% search_bifurcation.m % function [x_b,theta_b,lambda_b,iflag] = search_bifurcation(fun_name, ... % theta_0,theta_1,a,b,x0,num_lambda); % % This MATLAB routine searches for a bifurcation point % along a user defined path in parameter space, where % theta varies as a function of an order parameter lambda. % The inputs to the program are: % fun_name = a character string containing the name of the function % theta_0 = the value of the parameter vector at lambda = 0 % theta_1 = the value of the parameter vector at lambda = 1 % a = the initial lambda value at which one starts as an initial guess % b = the final lambda value that one uses as an initial guess % x0 = the initial guess of the solution % num_lambda = the number of lambda points to try as an initial guess % % The routine then attempts to search for a bifurcation point along % the line using values of the initial guesses of lambda between % a and b until a bifurcation point is found. The MATLAB routine % fsolve() is invoked to solve the augmented system. % % The program output is: % x_b = the solution at the bifurcation point % theta_b = the bifurcation point parameter vector % lambda_b = the lambda value at which the bifurcation point % was found % % K. J. Beers. MIT ChE. 9/23/03 function [x_b,theta_b,lambda_b,ifound_b] = search_bifurcation(fun_name, ... theta_0,theta_1,a,b,x0,num_lambda, ... max_Newton_iter,atol,max_line); % extract number of unknowns in system N = length(x0); Naug = N+1; % generate the vector of lambda initial guesses lambda_guess_vect = linspace(a,b,num_lambda); % now, iterate with different guesses of lambda until a % bifurcation point is found ifound_b = 0; for k=1:num_lambda lambda = lambda_guess_vect(k); x = x0; % get parameter vector theta = (1-lambda)*theta_0 + lambda*theta_1; % call fsolve with augmented set of equations Options = optimset('LargeScale','off','Display','off'); z0 = [x0; lambda]; [z_b,faug_b,exitflag] = fsolve(@sb_calc_f_aug, ... z0,Options,fun_name,theta_0,theta_1); x_b = z_b(1:N); lambda_b = z_b(N+1); if(exitflag > 0) % if sucessful ifound_b = 1; theta_b = (1-lambda_b)*theta_0 + lambda_b*theta_1; return; end end % if did not find a bifurcation point ifound_b = 0; x_b = x; theta_b = theta; lambda_b = lambda; return; %------------------------------------------------------------- % This routine computes the augmented function vector that is % to be zero at the bifurcation point. function f_aug = sb_calc_f_aug(z,fun_name,theta_0,theta_1); N = length(z) - 1; x = z(1:N); lambda = z(N+1); theta = (1-lambda)*theta_0 + lambda*theta_1; f_aug = zeros(N+1,1); f_aug(1:N) = feval(fun_name,x,theta); % now, estimate Jacobian from finite differences and % compute determinant Jac = sb_calc_Jac(fun_name,x,theta,f_aug(1:N)); det_Jac = det(Jac); f_aug(N+1) = det_Jac; return; %-------------------------------------------------------------- % This routine uses finite differences to evaluate the Jacobian % matrix of the function. function Jac = sb_calc_Jac(fun_name,x,theta,f); N = length(x); Jac = zeros(N,N); epsilon = sqrt(eps); for k=1:N x_new = x; x_new(k) = x(k) + epsilon; f_new = feval(fun_name,x_new,theta); Jac(:,k) = (f_new-f)/epsilon; end return;