% % CFD_path.m % % Run coupled finite differences on the model % % A -> B -> C % % with rates k1 = theta, k2 = 1. % % usage: CFD_path(tend,initial,theta,epsilon) % % tend = max time for the program to run % % initial = initial condition, ex: [30, 0, 0] % % theta = value of parameter at which we are computing the % derivative. % % epsilon = epsilon used in centered finite difference % % example: CFD_path(3,[300, 0, 0] , 0.5, 0.5*1/20) % function [x_plus,x_minus] = CFD_path(tend,initial,theta,epsilon) %rate constants k = [theta 1]; %maximum number of steps to allow. maxi = 10^4; % reaction vectors (there are three for each of the 2 reactions now) zeta = [-1 -1 0 0 0 0; 1 1 0 -1 -1 0; 0 0 0 1 1 0; -1 0 -1 0 0 0; 1 0 1 -1 0 -1; 0 0 0 1 0 1]; %initialize. T represents time and x is the state of the enlarged system. %x(1:3) gives the system with parameter theta + epsilon/2 whereas x(4:6) %gives the system with parameter theta - epsilon/2. T = 0; x = zeros([6,1]); x(1:3) = initial; x(4:6) = initial; %lambda1 are intensities for x(1:3). lambda2 are intensities for x(4:6) lambda1 = zeros([2,1]); lambda2 = zeros([2,1]); %As are intensities for actual representation. The first 3 %components represnets the first reaction, the second 3 represent % the second reaction. A = zeros([6,1]); %Set Poisson times and next Poisson Jumps. One for each of the 6 actual %reactions. Also set t_ks. Tk = zeros([6,1]); Pk = zeros([6,1]); t = zeros([6,1]); %set the first unit exponentials for the 6 Poisson processes. for j = 1:6 r = rand; Pk(j) = log(1/r); end for i=1:maxi-1, %update intensity functions lambda1(1) = (k(1)+epsilon/2)*x(1); lambda1(2) = k(2)*x(2); lambda2(1) = (k(1)-epsilon/2)*x(4); lambda2(2) = k(2)*x(5); %Now get the intensity functions for the 6 actual channels. A(1) = min(lambda1(1),lambda2(1)); A(2) = lambda1(1) - A(1); A(3) = lambda2(1) - A(1); A(4) = min(lambda1(2),lambda2(2)); A(5) = lambda1(2) - A(4); A(6) = lambda2(2) - A(4); %find which reaction happens next. for j = 1:6 if A(j) > 0 t(j) = (Pk(j) - Tk(j))/A(j); else t(j) = inf; end end %find the time of the next reaction. loc = 1; for ell = 2:6 if t(loc)>t(ell) loc = ell; end end %if all intensities are zero, we're done. quit. if t(loc) == inf; break end %update time T = T + t(loc); %update the state of the system. x = x + zeta(:,loc); %update the one Poisson process that fired. r = rand; Pk(loc) = Pk(loc) + log(1/r); %update local times. Tk = Tk + A*t(loc); if T >= tend break end end x_plus = x(1:3); x_minus = x(4:6); %diff = (x(2) - x(4))/h; %CPUTime = etime(clock,bt) end %The program %%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%% %%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%% %%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%