How to interpolate values less than 1? (interp1)

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

Answer 1

[Walter Roberson]

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.

[Annie]

Thanks for your answer, I tried with:

uRa = @(t) interp1(1:length(Ra), Ra, t, 'linear','extrap');

but it still does the same, I'm just testing it by plotting:

plot(t,uRa(t))

but it still gives me zero, and I tried it with a shorter version of the Ra data but it somehow not right:

What I need is that the interpolated vector, as called in the image, be the same as the original vector. What am I doing wrong?

Thanks for your help