[Date Prev][Date Next][Thread Prev][Thread Next][Date Index][Thread Index]
Re: Problem with Lsode..help
From: |
Carlo de Falco |
Subject: |
Re: Problem with Lsode..help |
Date: |
Mon, 27 Oct 2008 23:45:49 +0100 |
On 27/ott/08, at 21:28, genehacker wrote:
Hi Bryan, Here is the script
**********************
clf;
clc;
clear all;
1;
This command makes no sense, why is it there??
global x4;
this does not assign a value to x4
x5 = @(x1) k1max- ((k1max- k1min)/(1 + (f1 * x1/IC2)));
x6 = @(x1) k2max- ((k2max- k2min)/(1 + (0.01*f1 * x1/IC2)));
x7 = @(x1) k23min- ((k23min- k23max)/(1 + (0.01*f1 * x1/IC2)));
t0 = 8;
x8 = @(t,x2) ((t+0.001) .* (t<t0)) + (t0 .*(t==t0)) +
(t0*exp(0.1*t0)*exp(-0.1*t) .*(t > t0));
xdot = @(x,t) [ ( (x8(t,x(2)) * x4/pre_OC) - x5(x(1)) * x(1));
(x6(x(1)) * pre_OB - x7(x(1)) * x(2));
(x7(x(1)) * x(2) - k3 * x(3)) ];
the definition of xdot contains x4, so x4 should be initialized before
xdot is used
x = [0,0,0; 0,0,0];
len = 100;
time = linspace(0,len,len+1);
matrix = zeros(1,10);
matrix(1,10) = bmd = 100; % for bone initial state
a = 0.06; % OB weight
b = 0.1; % OC weight
for i=1:(length(time)-1)
seedtime = time(i:i+1);
[x,istate,msg] = lsode(xdot,[x(2,1),x(2,2),x(2,3)],seedtime);
lsode will repeatedly call the function xdot but, as x4 is not yet
defined, xdot(1)
will evaluate to the empty matrix [] so xdot(x,t) will be a 2-vector.
bmd = bmd + a*x(2,3) - b*x(2,1);
x4 = (vmax * pre_OC) / ( x8(time(i+1),x(2,2)) * (pre_OC + (x8(time(i
+1),x(2,2)) * (1 + (10 * x(2,2) / kI)) )) );
you won't get to this point as your script crashes on calling lsode
above, but this is also wrong. lsode as you called it above will
return a 3-vector which will overwrite the 2x3 matrix x
so x(2,2) will not exist after lsode returns. Plus vmax, pre_OC are
undefined.
matrix(i+1,:) = [time(i+1) x8(time(i+1),x(2,2)) x(2,1) x(2,2) x(2,3)
x5(x(2,1)) x6(x(2,1)) x7(x(2,1)) x4 bmd];
endfor
timerescale = 1;
plot(time/timerescale,10*matrix(:,3),"r"); hold on ;
plot(time/timerescale,10*matrix(:,5),"m"); hold on;
plot(time/timerescale,matrix(:,10),"b"); hold off;
xlabel("DAYS");
ylabel('Concentration');
axis([0 len]);
legend('OC','OB','bmd');
print -djpg foo.jpg;
*************************
c.
`