variable coefficients-help needed please!!!

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!!

[Oleg Komarov]

Please, format the code.

[Rose]

I have a system of 4 ODEs having variable coefficients, e.g k1, k2 etc. I am trying to solve the equations using ode15s. I defined the variable coefficients as piecewise constant functions and made sub-routines for all the non-constant coefficients. But I am getting an error.

Answer 1

[Oleg Komarov]

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.

[Rose]

thanks Oleg! If it includes 366th day and the 1st day into account, then i think its not good! because actually 366th day of the year is 1st day of the next year?

Anyway, I changed my code and it is:

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).^2)/(K(t)^2+y(2).^2)).*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

%Nested k3

function y = k3(t)

y = [0.03384,0.1984,0.1460,0.1900];

idx = logical(histc(t,[0,91.25,182.5,273.75,366]));

y = y(idx);

end

%Nested mu

function y = mu(t)

y = [0,500,1500,500];

idx = logical(histc(t,[0,91.25,182.5,273.75,366]));

y = y(idx);

end

%Nested d1

function y = d1(t)

y = [0.005263,0.02272,0.04,0.02272];

idx = logical(histc(t,[0,91.25,182.5,273.75,366]));

y = y(idx);

end

%Nested d2

function y = d2(t)

y = [0.005300,0.2,0.2,0.2];

idx = logical(histc(t,[0,91.25,182.5,273.75,366]));

y = y(idx);

end

%Nested r

function y = r(t)

y = [0.0045,0.0165,0.0165,0.0045];

idx = logical(histc(t,[0,91.25,182.5,273.75,366]));

y = y(idx);

end

%Nested K

function y = K(t)

y = [10000,20000,30000,20000];

idx = logical(histc(t,[0,91.25,182.5,273.75,366]));

y = y(idx);

end

%Nested k

function y = k(t)

y = [0.00003125,0.00003125,0.00003125,0.00003125];

idx = logical(histc(t,[0,91.25,182.5,273.75,366]));

y = y(idx);

end

%Nested alpha

function y = alpha(t)

y = [0.4784,0.4784,0.1321,0.1321];

idx = logical(histc(t,[0,91.25,182.5,273.75,366]));

y = y(idx);

end

%Nested gamma

function y = gamma(t)

y = [10^(-6),0.0002,0.0000002,0.0002];

idx = logical(histc(t,[0,91.25,182.5,273.75,366]));

y = y(idx);

end

end

%%%%%%%%%%%%%%%%%%%%%%%%%%%main.m%%%%%%%%

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);

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:

Error in ==> funl12>k at 73

idx = logical(histc(t,[0,91.25,182.5,273.75,366]));

Error in ==> funl12 at 16

dy(2,1) = (mu(t).*(y(2).^2)/(K(t)^2+y(2).^2)).*exp((-y(1)).*k)-(k3(t).*y(1).*y(2))./(y(2)+y(3)-1e-12)-(d1(t)+gamma(t).*y(4)).*y(2);

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);

[Rose]

i read your answer again and I understood about that 366th day thing. Thanks for that. Now, the problem is the error i am getting.

[Oleg Komarov]

in d(2,1) you call _k_ without supplying inputs. What do you want _k_ to return in that case?

function y = k(t)

if nargin == 0

y = ...something...

else

y = [0.00003125,0.00003125,0.00003125,0.00003125];

idx = logical(histc(t,[0,91.25,182.5,273.75,366]));

y = y(idx);

end

end

[Rose]

it works!!! thanks Oleg!! I am really grateful!! :-)

[Rose]

One more thing, if I want to run the code for years ans years, what should I do? Any suggestions?

[Rose]

I want to write these piecewise constant functions k1,k2,d1 etc into periodic functions. How can I do it?

Answer 2

[Walter Roberson]

if 0<t<=91.25;

is always true.

It is parsed as

if (0<t)<=91.25

and 0<t will have the value 0 (false) or 1 (true), both of which are less than or equal to 91.25

Recode to

if 0 < t && t <= 91.25

[Rose]

i made that change but still getting the same error!!