MATLAB: Variable coefficients-help needed please!!!

ode15surgent

let me write the exact code i am using:
%%%%%funl12.m%%%%%%%%%
function dy = funl12(t,y) 
global k1 k2 k3 mu d1 d2 r K k alpha gamma n;
dy(1,1) = (k1(t).*(y(4)-y(1)).*y(3))./(y(2)+y(3)-1e-12) - k2(t).*y(1).*(y(2)./(y(2)+y(3)-1e-12)) ;
dy(2,1) = (mu(t).*(y(2).^n)/(K(t)^n+y(2).^n)).*exp((-y(1)).*k)-(k3(t).*y(1).*y(2))./(y(2)+y(3)-1e-12)-(d1(t)+gamma(t).*y(4)).*y(2);
dy(3,1) = (k3(t).*y(1).*y(2))./(y(2)+y(3)-1e-12)-(d2(t)+gamma(t).*y(4)).*y(3);
dy(4,1) = r(t).*y(4).*(1-(y(4)./(alpha(t).*(y(2)+y(3)-1e-12)))); 
%%%%%%%%%%%%%%%%%%%k1.m%%%%%%%%%%%%
function y = k1(t) 
if 0<t<=91.25; 
y=0.04226; 
elseif 91.25<t<=182.5; 
y=0.1593; 
elseif 182.5<t<=273.75; 
y=0.1460; 
else 273.75<t<=365; 
y = 0.1489; 
end;
%%%%%%%%k2.m%%%%%%%%%%%%%%
function y = k2(t) 
if 0<t<=91.25; 
y=0.008460; 
elseif 91.25<t<=182.5; 
y=0.04959; 
elseif 182.5<t<=273.75; 
y=0.03721; 
else 273.75<t<=365; 
y=0.04750; 
end;
%%%%%%%%%%%%%
etc for all other parameters too.
%%%%%%%main.m%%%%%%%%%%%
n=2; 
k=0.00003125; 
finaltime = 365; 
tspan = [0 finaltime]; 
y0=[0;10000;0;0];%winter 
options = odeset('RelTol',1e-12,'AbsTol',1e-12);
sol = ode15s(@funl12,tspan,y0,options);
%[t,y] = ode15s('funpyo',tspan,y0);
X = sol.x; 
Y = (sol.y)';
plot(X,Y(:,1),'r-', X,Y(:,2),'g-', X,Y(:,3), X, Y(:,4)) 
xlabel('time') 
ylabel('Y') 
legend('m', 'X', 'y', 'M') 
figure
%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%
the error i am getting is:
??? Attempted to access k1(0); index must be a positive integer or logical.
Error in ==> funl12 at 4 
dy(1,1) = (k1(t).*(y(4)-y(1)).*y(3))./(y(2)+y(3)-1e-12) - k2(t).*y(1).*(y(2)./(y(2)+y(3)-1e-12)) ;
Error in ==> odearguments at 110 
f0 = feval(ode,t0,y0,args{:}); % ODE15I sets args{1} to yp0.
Error in ==> ode15s at 228 
[neq, tspan, ntspan, next, t0, tfinal, tdir, y0, f0, odeArgs, odeFcn, ...
Error in ==> simulate_lis2 at 66 
sol = ode15s(@funl12,tspan,y0,options);
%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%
i don't know what should i do!!

Best Answer

I redefined your code as:
function dy  =  funl12(t,y)
dy(1,1)  =  (k1(t).*(y(4)-y(1)).*y(3))./(y(2)+y(3)-1e-12) - k2(t).*y(1).*(y(2)./(y(2)+y(3)-1e-12)) ;
dy(2,1)  =  (mu(t).*(y(2).^n)/(K(t)^n+y(2).^n)).*exp((-y(1)).*k)-(k3(t).*y(1).*y(2))./(y(2)+y(3)-1e-12)-(d1(t)+gamma(t).*y(4)).*y(2);
dy(3,1)  =  (k3(t).*y(1).*y(2))./(y(2)+y(3)-1e-12)-(d2(t)+gamma(t).*y(4)).*y(3);
dy(4,1)  =  r(t).*y(4).*(1-(y(4)./(alpha(t).*(y(2)+y(3)-1e-12))));
      % Nested k1
      function y  =  k1(t)
          y   = [0.04226,0.1593,0.1460,0.1489];
          idx = logical(histc(t,[0,91.25,182.5,273.75,366]));
          y   = y(idx);
      end
      % Nested k2
      function y  =  k2(t)
          y   = [0.008460,0.04959,0.03721,0.04750];
          idx = logical(histc(t,[0,91.25,182.5,273.75,366]));
          y   = y(idx);
      end
  end
Remarks:
  1. Histc uses the following boundary behavior low <= t < up and this is consistent with k(0). In your code k(0) didn't return any value, thus the error you got. In this case it returns 0.04226 (if k1) or 0.008460 (if k2). Note that to include 365 I shifted it to 366 (assuming you'll never have values bigger equal than 366).
  2. Now I get an error on d(1,2) because I don't have code for K, k, gamma etc..
  3. No need to use globals, you can use nested functions which can see parent's workspace.
Related Question