Link to home
Start Free TrialLog in
Avatar of thebellman
thebellmanFlag for United States of America

asked on

Error when running ODE15s program.

I am using the ODE15s function in MATLAB and am getting an error.

 I am getting the following error (attached) when running the second part of the code.  Matlab is saying the varible T is not being used.  I'm a newbie to matlab and I can't figure out why it's complaining about T.

User generated image
Code for the function I used:
%Dua's function practice
function xdot = dua_glucose(t,x,u)

%t is time, x is the arrays, u is the insulin infusion rate

% (1) Input
% Exogenous Insulin infusion rate (mU/min)
% % u = 3;

% (2) States 
% Plasma Glucose Conc. (mg/dL)
G = x(1,1);
% Plasma Insulin Conc. (mU/L)
I = x(2,1);
% Plasma Insulin Conc. (1/min) in remote compartment
X = x(3,1);

%x is a matrix or (combination of arrays G,X,I). Consisting of a 3(variables) by n
%matrix which n is the time step size. 

% (3) Parameters From the original ex
%Assuming these conditions until literature
%is investigated further, For future will need
%Normal person values, and differing diabetic?
% Basal values of glucose and insulin conc. *Normal values
G_basal = 4.5; % mmol/L
%%X_basal = 15; % mU/L *Note there is not xb
I_basal = 15; % mU/L
% For a type-I diabetic These are parameters
P1 = 0.028735; % min-1
P2 = 0.028344; % min-1
P3 = 5.035e-5; % mU/L
%insulin distribution volume
V1 = 12; % L
%the fractional disappearance rate of insulin
n = 5/54; % 1/min
%rate of exogenously infused glucose
Fg = 0.5/t; %mg/min (divide by t?)
%glucose distribution space
Vg = 5; %dL 
% (4) Disturbances:
% Meal glucose disturbance (mg/dL min-1)
D = Fg/Vg;

%(5) Equations
Gdot = -P1*G - X*(G+G_basal) + D;
Xdot = -P2*(X) + P3*(I);
Idot = -n *(I+I_basal) + (u / V1);

% Vector the function produces
xdot = [Gdot; Xdot; Idot];

Open in new window


Code for the ODE portion:
clear;
clc;

% Steady State Initial Conditions for the States
% Basal values of glucose and insulin conc.
G_ss = 4.5; % mg/dL
I_ss = 15; % mU/L
X_ss = 15; % 1/min

%x_ss = [G_ss;I_ss;X_ss]; 
x_0 = [G_ss I_ss X_ss];

% Steady State Initial Condition for the Control
u_ss = 16.667; % mU/min

% Steady State for the Disturbance *Assuming
d_ss = 0; % mmol/L-min

% Final Time (min)*Leaving the same 
%as example in order to compare
tf = 100; t0 = 0;
t_step = 100; tt_step = 5;
dt = (tf - t0)/t_step;
t = 0;

% initialization
uu = zeros(1,t_step+1)';
G = zeros(1,t_step+1)';
X = zeros(1,t_step+1)';
I = zeros(1,t_step+1)';

% Control Input
%Why is this like this? 
for i=1:t_step+1
    uu(i) = 60*sin(pi*(i-1)/t_step)+10;
end
% u = 100;

%Why is tspan written like this? 
for i=1:t_step+1
    tspan = (t:(dt/tt_step):t+dt);
    t = t + dt; %where did you get this from?
    u = uu(i);
%[t,x] = ode15s('dua_glucose',[0 tf],x_ss);
[t,x] = ode15s(@(t,x) dua_glucose(t,x,u),tspan,x_0,u);

% Separate out the state values that are in the array 
G(i) = x(tt_step+1,1);     
X(i) = x(tt_step+1,2);     
I(i) = x(tt_step+1,3);     
x_0 = [G(i),X(i),I(i)];
end

t = 0:dt:tf;

Graphing

figure(1);
plot(t,G);
legend('Glucose');
xlabel('Time (min)');
ylabel('mg/L');

figure(2);
plotyy(t,I,t,X);
legend('Plasma Insulin','Remote Compartment');
xlabel('Time (min)');
ylabel('mU/L');

Open in new window

This question needs an answer!
Become an EE member today
7 DAY FREE TRIAL
Members can start a 7-Day Free trial then enjoy unlimited access to the platform.
View membership options
or
Learn why we charge membership fees
We get it - no one likes a content blocker. Take one extra minute and find out why we block content.