variable coefficients-help needed please!!!

1 view (last 30 days)
Rose
Rose on 25 Apr 2011
Edited: Rose on 22 Feb 2016
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!!
  2 Comments
Rose
Rose on 25 Apr 2011
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.

Sign in to comment.

Accepted Answer

Oleg Komarov
Oleg Komarov on 25 Apr 2011
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.
  6 Comments
Rose
Rose on 25 Apr 2011
I want to write these piecewise constant functions k1,k2,d1 etc into periodic functions. How can I do it?

Sign in to comment.

More Answers (1)

Walter Roberson
Walter Roberson on 25 Apr 2011
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
  1 Comment
Rose
Rose on 25 Apr 2011
i made that change but still getting the same error!!

Sign in to comment.

Community Treasure Hunt

Find the treasures in MATLAB Central and discover how the community can help you!

Start Hunting!