ODE45 specific time point settings

Hello!

I'm trying to stimulate BIG model in different situations (Beta-cells , insulin , glucose).

The following function is very simple. I'm able to change the parameter 's' at t>100.

%%% BIG Model %%%

function yprime = big(t,y) 


s=1;
delta=1;
q=1;


if t>100
  s=0.2;
end

    
yprime = zeros(3,1);
yprime(1) = m - s*y(2)*y(1);%dG/d
yprime(2) = (q * y(3) * (y(1))^2) - delta * y(2); %dI/dt
yprime(3) = y(3)*0.01*(y(1)-5);%dB/dt

The second thing that I need to do is to change m=2 at t=90, t=110, t=200 (rest of the time m=1)

The problem is that the ODE is choosing the time points by itself and it skips these specific time point. How can I overcome this?

This is the code that i'm using to solve this.

clear all
clc

[t,y] = ode45(@(t,y) big(t,y) , [0:1:400] ,[5,1/5,1/125]);



subplot(1,3,1)
plot(t,y(:,1))
title('Glucose')
xlabel('Time')
subplot(1,3,2)
plot(t,y(:,2))
title('Insulin')
xlabel('Time')
subplot(1,3,3)
plot(t,y(:,3))
title('Beta-Cells')
xlabel('Time')

Answer 1

[Star Strider]

I verified that this version works, however I do not see any significant change in the integrated results.

These added lines:

tcnd = (((t>=90) & (t<91)) | ((t>=110) & (t<111)) | ((t>=200) & (t<201)));
m = tcnd*2 + (~tcnd)*1;

‘toggle’ ‘m’ between ‘1’ and ‘2’ at the appropriate times. An additional fprintf call allows verification.

The (Slightly Revised) Code —

function yprime = big(t,y) 
s=1;
delta=1;
q=1;
if t>100
  s=0.2;
end

tcnd = (((t>=90) & (t<91)) | ((t>=110) & (t<111)) | ((t>=200) & (t<201)));
m = tcnd*2 + (~tcnd)*1;
% fprintf('t = %f\t\tm = %d\n',t,m)                   % Un-Comment To See dm’ At Each Time
yprime = zeros(3,1);
yprime(1) = m - s*y(2)*y(1);%dG/d
yprime(2) = (q * y(3) * (y(1))^2) - delta * y(2); %dI/dt
yprime(3) = y(3)*0.01*(y(1)-5);%dB/dt
end

The Plot —

The only discernable difference are the ‘spikes’.

[Eran Sandman]

I need to see 3 spikes.

So complicated :/

[Star Strider]

They are all present.

The spike at 90 does not appear because of the if block defining ‘s’ keeps the ‘s’ at 1 until t=100. (I verified that they all exist at the correct locations.) The ‘Beta Cells’ plot will not show any spikes.

Plot this additional figure (after the subplot array) to show all of them (or at least where they should be):

spikes = @(t) (((t>=90) & (t<91)) | ((t>=110) & (t<111)) | ((t>=200) & (t<201)));
mv = spikes(t);
mplot = mv*2 + (~mv)*1;

figure
plot(t,y(:,1))
hold on
plot(t, mplot*5, '-r')
hold off
title('Glucose')
xlabel('Time')

producing:

Ths ‘spikes’ function (and the ‘mv’ and ‘mplot’ assignments) use the same logical array as I used in the differential equation code.

So my code creates all of them correctly, and they would appear in the plot if allowed. The spike at 90 is not allowed because of the way ‘s’ is defined.

.

[Eran Sandman]

Thanks for you answer.

I still need to see spike like it's suppose to be. How can we overcome the situation with s?

[Star Strider]

My pleasure!

I have absolutely no idea about ‘s’ because I don’t recognise the model you’re using.

That aside, my code does exactly what you requested.

If my Answer helped you solve your problem, please Accept it!

.

[Eran Sandman]

Evenutally i changed the conditions of 't' to nearby points which likley the intergration goes through them. So it worked well.

[Star Strider]

I’m happy you got it sorted.

I still don’t understand what the problem is with ‘s’. In spite of my experience with metabolic models (specifically glucose regulatory models), I don’t recognise the one you’re using.

Answer 2

[darova]

Try tolerance

function yprime = big(t,y) 
m = 1;
if any( abs([90 110 200]-t)<0.5 )
    m = 2;
end

Answer 3

[Boxn Hen]

Do you change m at one point like t=90 ?Ode steps small size ,if the tspan is [0,1] and you set m=1 when t=0,while m=2 in other time,it will be meaningless.What I experience is effective if you change m for a period.

function yprime = big2(t,y) 
s=1;
delta=1;
q=1;
m=1;
if t>100
  s=0.2;
end
 if t>90
     m=2;
 end
yprime = zeros(3,1);
yprime(1) = m - s*y(2)*y(1);%dG/d
yprime(2) = (q * y(3) * (y(1))^2) - delta * y(2); %dI/dt
yprime(3) = y(3)*0.01*(y(1)-5);%dB/dt
end