variable coefficients-help needed please!!!

Question

Rose 2011-4-25

0
链接

此问题的直接链接

https://ww2.mathworks.cn/matlabcentral/answers/6150-variable-coefficients-help-needed-please

编辑： Rose 2016-2-22

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 个评论
显示无隐藏无

Oleg Komarov 2011-4-25

Please, format the code.

Rose 2011-4-25

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 2011-4-25

0
链接

此回答的直接链接

https://ww2.mathworks.cn/matlabcentral/answers/6150-variable-coefficients-help-needed-please#answer_8539

在 MATLAB Online 中打开

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:

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).
Now I get an error on d(1,2) because I don't have code for K, k, gamma etc..
No need to use globals, you can use nested functions which can see parent's workspace.

6 个评论
显示 4更早的评论隐藏 4更早的评论

Rose 2011-4-25

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

%%%%%%%%%%%%%%%%%%%%%%%%%%%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 2011-4-25

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

Rose 2011-4-25

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

请先登录，再进行评论。

Answer 2

Walter Roberson 2011-4-25

0
链接

此回答的直接链接

https://ww2.mathworks.cn/matlabcentral/answers/6150-variable-coefficients-help-needed-please#answer_8538

在 MATLAB Online 中打开

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 个评论
显示 -1更早的评论隐藏 -1更早的评论

Rose 2011-4-25

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

请先登录，再进行评论。

variable coefficients-help needed please!!!

2 个评论
显示无隐藏无

采纳的回答

6 个评论
显示 4更早的评论隐藏 4更早的评论

更多回答（1 个）

1 个评论
显示 -1更早的评论隐藏 -1更早的评论

另请参阅

类别

标签

Community Treasure Hunt

variable coefficients-help needed please!!!

2 个评论 显示 无隐藏 无

采纳的回答

6 个评论 显示 4更早的评论隐藏 4更早的评论

更多回答（1 个）

1 个评论 显示 -1更早的评论隐藏 -1更早的评论

另请参阅

类别

标签

Community Treasure Hunt

2 个评论
显示无隐藏无

6 个评论
显示 4更早的评论隐藏 4更早的评论

1 个评论
显示 -1更早的评论隐藏 -1更早的评论