Re: matlab to MMA, need some help converting this code.
From: sean_incali (sean_incali_at_yahioo.com)
Date: 11/12/04
- Next message: Brad Cooper: "Re: MuPad 3.0 extraction of individual eigenvectors"
- Previous message: Peter Pein: "Re: Mathematica digit handling query"
- In reply to: sean kim: "matlab to MMA, need some help converting this code."
- Messages sorted by: [ date ] [ thread ]
Date: Fri, 12 Nov 2004 06:34:51 +0000 (UTC)
No takers?
I really could use soem startin goff point. I don't much about matlab.
Can someone suggests some ways to convert this code?
thanks in advance
sean
On 10 Nov 2004 09:37:47 -0800, sean kim wrote:
>does anyone know both mathematica AND matlab?
>
>if so can any of you help me out with this?
>
>i have received some codes written in matlab, and i would like to
>convert it to mma.
>
>anyone up for the challenge *and* have the time and energy? It will
be
>nice to see line by line explanation of the conversion.
>
>following is the code which I'm having problems understanding.
>
>
>thanks for any insights
>
>sean
>
>
>-------
>
>
>function result= gillespie(tf, v0)
>% function result= gillespie(tf, v0)
>%
>% runs the Gillespie algorithm for constitutive expression up to
>% a time tf. the value of v0 may be specified but has a default value
>% of 0.01
>
>if nargin == 1
> v0= 0.01;
>end;
>
>% initial amounts
>D= 1;
>M= 0;
>N= 0;
>t= 0;
>% initialize other variables
>j= 2;
>nr= 4; % number of reactions
>
>% create some large result matrix (this is not essential but
>% speeds up the program)
>result= zeros(100000, 4);
>
>% store data for first time point
>result(1,:)= [t D M N];
>
>while t < tf
>
> % conversion
> s(1)= D;
> s(2)= M;
> s(3)= N;
>
> % store data
> result(j,:)= [t D M N];
> j= j+1;
>
> % calculate propensities
> a= calpropensities(s, v0);
> a0= sum(a);
>
> % generate two random numbers
> r= rand(1,2);
>
> % calculate time of next reaction
> dt= log(1/r(1))/a0;
>
> % calculate which reaction is next
> for mu= 1:nr
> if sum(a(1:mu)) > r(2)*a0
> break;
> end;
> end;
>
> % carry out reaction
> t= t+dt;
> if mu == 1
> M= M+1;
> elseif mu == 2
> M= M-1;
> elseif mu == 3
> N= N+1;
> elseif mu == 4
> N= N-1;
> end;
>
>end;
>
>% store data for last time point
>result(j,:)= [t D M N];
>
>% remove any excess zeros from result
>result(j+1:end,:)= [];
>
>
>
>%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%
>function a= calpropensities(s, v0)
>
>% rates
>v1= 0.04;
>d0= log(2)/180;
>d1= log(2)/3600;
>
>% conversion
>D= s(1);
>M= s(2);
>N= s(3);
>
>% propensities
>a(1)= v0*D;
>a(2)= d0*M;
>a(3)= v1*M;
>a(4)= d1*N;
- Next message: Brad Cooper: "Re: MuPad 3.0 extraction of individual eigenvectors"
- Previous message: Peter Pein: "Re: Mathematica digit handling query"
- In reply to: sean kim: "matlab to MMA, need some help converting this code."
- Messages sorted by: [ date ] [ thread ]
Relevant Pages
|
