Effacer les filtres
Effacer les filtres

variable coefficients-help needed please!!!

2 vues (au cours des 30 derniers jours)
Rose
Rose le 25 Avr 2011
Modifié(e) : Rose le 22 Fév 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 commentaires
Oleg Komarov
Oleg Komarov le 25 Avr 2011
Please, format the code.
Rose
Rose le 25 Avr 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.

Connectez-vous pour commenter.

Réponse acceptée

Oleg Komarov
Oleg Komarov le 25 Avr 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 commentaires
Rose
Rose le 25 Avr 2011
One more thing, if I want to run the code for years ans years, what should I do? Any suggestions?
Rose
Rose le 25 Avr 2011
I want to write these piecewise constant functions k1,k2,d1 etc into periodic functions. How can I do it?

Connectez-vous pour commenter.

Plus de réponses (1)

Walter Roberson
Walter Roberson le 25 Avr 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 commentaire
Rose
Rose le 25 Avr 2011
i made that change but still getting the same error!!

Connectez-vous pour commenter.

Catégories

En savoir plus sur Programming dans Help Center et File Exchange

Community Treasure Hunt

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

Start Hunting!

Translated by