https://github.com/Data2Dynamics/d2d
Tip revision: 83d6eedd5f6c33ae94dbce6c4e620d2f1bb57b1f authored by Joep Vanlier on 09 April 2015, 08:41:55 UTC
Forgot to update C version number
Forgot to update C version number
Tip revision: 83d6eed
arMomentODE2.m
% calculate moment ODEs
%
% arMomentODE2(m)
%
% calculate up to 2nd moment
%
% A Moment Closure Method for Stochastic Reaction Networks
% Chang Hyeong Lee, Kyeong-Hun Kim, Pilwon Kim, February 25, 2009
function arMomentODE2(m)
global ar;
N = ar.model(m).N;
fv = sym(ar.model(m).fv);
x = cell(1,length(ar.model(m).x));
covar = cell(length(ar.model(m).x));
extra = cell(length(ar.model(m).x),length(ar.model(m).x),length(ar.model(m).x));
for jx = 1:length(x)
x{jx} = sprintf('%s', ar.model(m).x{jx});
for jx2 = 1:length(x)
covar{jx,jx2} = sprintf('COVAR_%s_%s', ar.model(m).x{jx}, ar.model(m).x{jx2});
for jx3 = 1:length(x)
extra{jx,jx2,jx3} = sprintf('EXTRA_%s_%s_%s', ar.model(m).x{jx}, ...
ar.model(m).x{jx2}, ar.model(m).x{jx3});
end
end
end
x = sym(x);
covar = sym(covar);
for jx = 1:length(x)
fprintf('%s\n', char(x(jx)));
end
for jx = 1:length(x)
for jx2 = 1:length(x)
fprintf('%s\n', char(covar(jx,jx2)));
end
end
dxdt = sym(zeros(size(x)));
for ji = 1:length(x)
for jk = 1:length(fv)
tmp = sym(0);
for jl = 1:length(x)
for jm = 1:length(x)
tmp = tmp + 0.5 * diff(diff(fv(jk), x(jl)), x(jl)) * covar(jl,jm);
end
end
dxdt(ji) = dxdt(ji) + N(ji,jk) * (fv(jk) + tmp);
end
fprintf('"%s"\n', char(dxdt(ji)));
end
dcovardt = sym(zeros(size(covar)));
for ji = 1:length(x)
for jj = 1:length(x)
for jk = 1:length(fv)
tmp1 = sym(0);
tmp2 = sym(0);
tmp3 = sym(0);
tmp4 = sym(0);
tmp5 = sym(0);
for jl = 1:length(x)
tmp1 = tmp1 + diff(fv(jk), x(jl)) * covar(jj,jl);
tmp2 = tmp2 + diff(fv(jk), x(jl)) * covar(ji,jl);
for jm = 1:length(x)
tmp3 = tmp3 + 0.5*diff(diff(fv(jk), x(jl)), x(jm)) * covar(jl,jm);
tmp4 = tmp4 + 0.5*diff(diff(fv(jk), x(jl)), x(jm)) * extra(jj,jl,jm);
tmp5 = tmp5 + 0.5*diff(diff(fv(jk), x(jl)), x(jm)) * extra(ji,jl,jm);
end
end
dcovardt(ji,jj) = dcovardt(ji,jj) + N(ji,jk)*tmp1 + ...
N(jj,jk)*tmp2 + ...
N(ji,jk)*N(jj,jk)*(fv(jk) + tmp3) + ...
N(ji,jk)*tmp4 + ...
N(jj,jk)*tmp5;
end
fprintf('"%s"\n', char(dcovardt(ji,jj)));
end
end