Re: Problem with Lsode..help

[Carlo de Falco]

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.
`