Problems with oder45, for or if

5 visualizzazioni (ultimi 30 giorni)
Ewa Kozub
Ewa Kozub il 18 Lug 2013
I have problem with my program. Seemingly everything is working well but when i run ode45 I get wrong results. My program doesn't see change of Ginf. I really don't know why. Ginf in two first iterations is very high and i should see this on graph. but it doesn't see
function F=model_magda(t,y)
F=zeros(2,1);
global Vg Vi Rinf Cinf Gg gamma lambda alfa C1 a1 Rm Ki Kg Tinf UFR Gd Vglucose Vinsulin Ginf
for t=1:61
if t-1<Tinf
Ginf(t)=Rinf*Cinf;
else
Ginf(t)=0;
end;
Vglucose(t)=Vg-UFR*(t-1);
Vinsulin(t)=Vi-UFR*(t-1);
F(1)=(Gg-gamma*y(1)*y(2)-lambda*y(1)+Ginf(t)-Kg*(y(1)-Gd)+y(1)*UFR)/Vglucose(t);
F(2)=(Rm/(1+exp((C1-y(1))/a1))-alfa*y(2)-Ki*y(2))/Vinsulin(t);
t=t+1;
end
-------------------------------------------------
global Vg Vi Rinf Cinf Gg gamma lambda alfa C1 a1 Rm Ki Kg Tinf UFR m Gd Vglucose Vinsulin yo meani meang time to Ginf
UFR=3/240; %0.004166667; %[L/min]
m=78; %[kg]
Gd=5.615292712;
Rinf=60/1000 %[L/min]
Vinf=0.5*m*100/(33*1000); %[L]
Cinf=33/180.16*1000*1000/100; %[mmol/L]
Tinf=Vinf/Rinf %[min]
Minf=Tinf*Rinf*Cinf;
Minfg=Minf*180.16/1000; %Ginfa=Rinf*Cinf
%----------srednia---------------
Go=6.63;
Io=8.89;
Gg=0.621129599; %[mmol/min]
gamma=0.002993997; %{L^2/mU*min]
lambda=0.101455909; %[L/min]
alfa=0.095577674; %[L/min]
C1=25.51490467; %[mmol/L]
a1=1.194742573; %[mmol/L]
Rm=79.44266365; %[mU/min]
Vg=8.764594059; %[L]
Vi=6.673505545; %[L]
time=[ 0 10 20 30 40 60];
meang=[ 6.63 18.67 14.19 11.59 10.9 9.08];
meani=[ 8.89 25.07 17.54 13.05 10.99 8.76];
Ki=63/1000; %[L/min]
Kg=223.2/1000; %[L/min]
[to yo]=ode45('model_magda',[0:60],[Go Io]);
[to yo(:,1) yo(:,2)]
figure(1)
plot(to,yo(:,1),'b')
hold on
plot(time,meang,'b--x')
grid on;
xlabel('Time[min]');
ylabel('glucose concetration in blood [mmol/l]');
legend('simulation','dane');
figure(2)
plot(to,yo(:,2),'r')
hold on
plot(time,meani,'r--x')
grid on;
xlabel('Time [min]');
ylabel('insulin concetration in blood [mU/l]');
legend('simulation','dane');
figure(3)
plot(yo(:,1),yo(:,2),'k')
hold on
plot(meang,meani,'k--x')
legend('simulation','dane');

Risposte (0)

Categorie

Scopri di più su Frequently-used Algorithms in Help Center e File Exchange

Tag

Community Treasure Hunt

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

Start Hunting!

Translated by