How to interpolate values less than 1? (interp1)
3 views (last 30 days)
Show older comments
Hello, I'm working with some data that I need to introduce in an equation, I've done this with interp1 before, but now the data values are way smaller that the passed examples I did, values goes around zero but not greater than one. At first I thought that it was the interpolation method I was using (I tried nearest, linear and spline), then I scaled the vector data and it was the same, here's the code:
function [t,x]=Cobelli2007_Daily_T2DM_Offline9E(CIx1, CIx2, CIx3, CIx4, CIx5, CIx6, CIx7, CIx8, CIx9, ke1, kp2, alfa, betha, tf, Ra)
[Ra]=Cobelli2006_T2DM_Daily(tf);
Ra=ans;
tspan = [0;tf];
options = odeset('MaxStep',0.1);
x0 = [CIx1; CIx2; CIx3; CIx4; CIx5; CIx6; CIx7; CIx8; CIx9];
uRa = @(t) interp1(1:length(Ra), Ra, t, 'spline');
[t,x] = ode45(@(t,x) f(t,x, ke1, kp2, alfa, betha,tf,uRa), tspan,x0,options);
assignin('base','t',t);
assignin('base','x',x);
function dxdt = f(t,x, ke1, kp2, alfa, betha,tf,uRa)
H = [8;12;20];
M = [31; 1; 1;];
H = (H.*60) + M;
tpw = 0.1;
n = 3;
VG=1.49;
k1=0.042;
k2= 0.071;
VI= 0.04;
m1=0.379;
m2=0.673;
m4=0.269;
m5=0.0526;
m6=0.8118;
HEb=0.6;
a=0.00006;
b=0.68;
c=0.00023;
d=0.09;
kp1=3.09;
kp3=0.005;
kp4=0.0786;
ki=0.0066;
Fcns=1;
Vm0=4.65;
if ((0<=t) & (t<=H(n)))
Vmx=0.034;
else if ((H(n)<=t) & (t<=(tpw+H(n))))
Vmx=0.034*0.75;
else if ((tpw+H(n)<=t) & (t<=tf))
Vmx=0.034;
end
end
end
Km0=466.21;
Kmx=0;
p2u=0.0840;
ke2=269;
K=099;
if ((0<=t) & (t<=H(n-1)))
btha=betha
else if (((H(n-1)<=t) & (t<=(tpw+H(n-1)))) | ((H(n)<=t) & (t<=(tpw+H(n)))))
btha=betha*0.75;
else if ((((tpw+H(n-1)<=t) & (t<=H(n)))) | ((tpw+H(n)<=t) & (t<=tf)))
btha=betha;
end
end
end
gamma=0.5;
EGP=kp1-kp2*x(1)-kp3*x(6)-kp4*x(8);
Uii= Fcns;
if x(1)>ke2
E=ke1*(x(1)-ke2);
else
E=0;
end
Uid= ((Vm0+Vmx*x(7))*x(2))/((Km0+Kmx*x(7))+x(2));
S=gamma*x(8);
HE=-m5*S + m6;
m3= (HE*m1)/(1-HE);
G=x(1)/VG;
I=x(4)/VI;
Ib=50;
dG=(EGP+uRa(t)-Uii-E-k1*x(1)+k2*x(2))/VG;
Sb= 2;
if dG > 0
Spo=x(9)+K*dG+Sb;
else
Spo=x(9)+Sb;
end
Gb= 100;
h=Gb;
if btha*(G-h) >= -Sb
X12= -alfa*(x(9)-btha*(G-h));
else
X12=-alfa*x(9)-alfa*Sb;
end
dxdt = [
EGP+uRa(t)-Uii-E-k1*x(1)+k2*x(2); %x(1)
-Uid+k1*x(1)-k2*x(2); %x(2)
-(m1+m3)*x(3)+m2*x(4) + S; %x(3)
-(m2+m4)*x(4)+m1*x(3); %x(4)
-ki*(x(5)-I); %x(5)
-ki*(x(6)-x(5)); %x(6)
-p2u*x(7)+p2u*(I-Ib); %X(t) %x(10)
-gamma*x(8)+Spo; %x(11)
X12 %x(12)
];
where Ra is attached, along with initial conditions and missing parameters, and tf=1800, the result I get after the "interpolation" is everything zero. What can I do to make it work? Thanks
0 Comments
Answers (1)
Walter Roberson
on 14 Feb 2018
In a few different places, your function is discontinuous. ode45() cannot handle discontinuous functions. If you are lucky ode45() will just end up thinking about it for a while and then complaining that it has hit a discontinuity, thus giving you a hint to review and fix your code; if you are not lucky then ode45 will merely give you the wrong answer.
Your code needs to stop the integration at each discontinuity, and restart ode45 there.
Your obvious discontinuities are at t = 12 and t = 20, but you also have discontinuities at t - 12.1 and t = 20.1. These ones are relatively easy to deal with by adjusting your tspan boundaries. But you have a discontinuity based upon btha*( x(1)/1.49 - 2 ) where btha has discontinuities at the mentioned times, and because that deals with the current value of one of the input x, you are going to need to use an event function to detect those crossings and tell it to terminate the integration at that point.
With regards to your interpolation question:
You have
uRa = @(t) interp1(1:length(Ra), Ra, t, 'spline');
As you are interpolating at t, and t is time, the implication is that Ra is a vector of values corresponding to times 1, 2, 3, 4 ,... length(Ra). But t starts at 0, and MATLAB will for sure try a number of time steps near 0, so you are asking for interpolation before the declared initial data point "1". MATLAB cannot interpolate before the given data coordinates or after the given data coordinates. To get values before the given data coordinates or after the given data coordinates, you need to specify 'extrap', preferably explicitly saying how you want extrapolation to take place.
See Also
Community Treasure Hunt
Find the treasures in MATLAB Central and discover how the community can help you!
Start Hunting!