MATLAB Matlab Sim Code for ODE & Gillespie for Reversible Reaction

AI Thread Summary
The discussion focuses on a MATLAB code implementation for simulating predator-prey dynamics using both ordinary differential equations (ODE) and the Gillespie algorithm. Initial population densities are defined for predators (A) and prey (B), with specific reaction rates set for the model. The ODE simulation calculates population changes over time, while the Gillespie method uses stochastic processes to model the reactions. The user expresses difficulty in uploading a PDF that explains the chemical reactions involved, citing file size issues. There is also a suggestion to create a code repository for better collaboration and version control among forum members, highlighting the usefulness of shared code contributions. The conversation includes attempts to resolve the file upload problem, with advice on zipping the file to reduce size.
DrWahoo
Messages
45
Reaction score
0
Here is the code you input into matlab. Aini etc, are the initial values of the population densities. A for predator, B for Prey.
Code:
% example for ODE and Gillespie

% one reversible reaction
 b1 = .033;
bo = .00047;
a1 = .0022;
ao = .00055;

Aini = 5168;
Bini = 34;

%% Basic ODE simulation
dt = 0.01;
Tmax = 1;

TV = (0:dt:Tmax)';
NT = length(TV);

ABVals = zeros(size(TV,1),2);

A=Aini;
B=Bini;
ABVals(1,:) = [A, B];

% loop over times
for it=2:NT
    dAdt = - b1*A*B + bo*A;
    dBdt = a1*B*A -ao*B 
    A = A + dAdt * dt;
    B = B + dBdt* dt;
    
    ABVals(it,:) = [A, B];
    
    
end

figure(1)
clf;
plot(TV,ABVals);
xlabel 'T';
legend('A','B');
hold on

return;
%% Gillespie

% loop goes over time steps

%Tmax = 1;
MaxSteps = 1000;

% define a large array to record everything
ABRecord = zeros(MaxSteps,2);
TRecord = zeros(MaxSteps,1);ABCount = [Aini,Bini];

t=0;TRecord(1) = t;
ABRecord(1,:) = ABCount;
istep=2;
while (t<= Tmax  && istep<=MaxSteps),
    
    A=ABCount(1);
    B=ABCount(2);
    
    % propensities ('gamma')
    PropElkEaten = b1 * A*B;
    PropWolvesBorn = a1*A*B;
    PropWolvesDie = ao*B;
    PropElkBorn = bo*A;
    
    % distribution of next reaction times follows exp( - gamma * t)
    NextFwdTime = exprnd(1/PropElkEaten,1);
    NextRevTime = exprnd(1/PropWolvesBorn,1);
    NextWolvesDie = exprnd(1/PropWolvesDie,1);
    NextElkBorn = exprnd(1/PropElkBorn,1);
    
    if NextFwdTime < NextRevTime,
        % forward reaction (A->B) happens
        t = t + NextFwdTime;
        ABCount = ABCount + [-1,1];
    else
        % reverse reaction (A->B) happens
        t = t + NextRevTime;
        ABCount = ABCount + [1,-1];
    end
    
    % record
    ABRecord(istep,:) = ABCount;
    TRecord(istep) = t;
    
    % done with this step, increment the step count
    istep=istep+1;
    
end

if istep-1<MaxSteps,
    TRecord = TRecord(1:istep-1);
    ABRecord = ABRecord(1:istep-1,:);
end

figure(1)
stairs(TRecord,ABRecord)
legend('A - ODE','B - ODE','A - Gillespie','B - Gillespie');
 
Physics news on Phys.org
Bit hard to read...
 
I am trying to upload a PDF of what the chemical reaction is doing to help the logic in the iterations. Maybe Mark can help me upload it. Its saying my file size is to large.
 
This was part of an exercise I did sometime last year I think. It is larger than it needs to be since I was solving the equations with a manually implemented RK4 scheme.

EDIT: Wow. The forum recognizes abbreviations, and I didn't even have to do anything!
 

Attachments

I think the code contributions in this and other recent threads are useful and unselfish, which is nice.

Would it be an idea to set up a small MHB code repository (perhaps using git or svn)? This would allow for easier forking, editing and keeping track of authors and changes. (In short, it would help to maintain an overview.)

Of course, I cannot oversee how feasible this is, but the possibility crossed my mind.
 
DrWahoo said:
I am trying to upload a PDF of what the chemical reaction is doing to help the logic in the iterations. Maybe Mark can help me upload it. Its saying my file size is to large.

Did you try zipping the file before uploading it?
 
MarkFL said:
Did you try zipping the file before uploading it?

Yeah, still won't upload. Its 3.1 Megs.
I will shoot you an email Mark.
 

Similar threads

Back
Top