How to use else if with ODE45 solver?

%define ODEs
%DOX conc in blood
if t<Tiv 
    diffeqs(1,1) = (-k_clr_DOX*C_DOX_Bl)+ R_DOX
else 
    R_DOX = 1 
    R_DOX = 0
end
 
%BEV conc in blood
if t<Tiv
    diffeqs(2,1) = (-k_clr_BEV*C_BEV_Bl)+ R_BEV
else
    R_BEV = 1 
    R_BEV = 0
end

This is what I have and it is showing an error message when I try to solve this?

Answer 1

[Walter Roberson]

The mathematics used in the RK methods such as ode45(), require that for any one call the function be continuous and its first two derivatives are continuous. If the second derivative is not continuous, the most likely result is that the function will return incorrect results without noticing. If the first derivative is not continuous, sometimes it will notice and and abort the integration and other times it will not notice and will continue producing wrong results without noticing.

If for any one call you cross the time boundary between < Tiv and >= Tiv, then your function is not continuous, and if you are lucky you will get a message about it; if you are unlucky then it will return an answer without telling you that the answer is wrong.

We notice that if in the case where t >= Tiv, you do not assign anything to diffeqs(1,1) or diffeqs(2,1), which could potentially be a problem (depending what else you do.)

What you should do is run the code once with time span 0 to Tiv, and then take the final conditions from that and use them as the input boundary conditions for running a second call with time span Tiv to whatever.

[Skye Cameron]

Would you be able to give me an example as to how I would write this code? I feel slightly lost!

[Skye Cameron]

Full code attached below, I have diffeqs specified earlier in the code

%Run simulation
%define the range of time
range = (0:3600*1); %24 hour
 
%define inital conditions of ODEs
ICs = [0,0,0,1,0];
 
% call the user-defined functions
% pass the range of independent variables
% pass the intitial conditions
[tsol,varsol] = ode45(@ode_sys, range, ICs);
 
figure 
hold on 
box on
plot (tsol/3600, varsol(:,1),"LineWidth",2)
plot (tsol/3600, varsol(:,2),"LineWidth",2)
plot (tsol/3600, varsol(:,3),"LineWidth",2)
plot (tsol/3600, varsol(:,5), "LineWidth",2)
title ('Concentration as a function of time')
legend ('Concentration DOX In Blood','Concentration BEV In Blood','Concentration DOX In Tumour','Normalised Concentration BEV in Blood')
set(gca,'xtick', (0:3:24)) %set scale on x-axis
xlim([0,24])
xlabel ('Time (h)')
ylabel ('Concentration (kg/m^3)')
hold off
 

Define the functions

function diffeqs = ode_sys(t,var)
%define variables
C_DOX_Bl = var(1)
C_BEV_Bl = var(2)
C_DOX_T = var(3)
phi = var(4)
Cs_BEV_Bl = var(5)
 
%define the parameters FOR DOX
a = -4.10e-6 %alpha in 1/s
B = 1.22e-5 %beta in 1/s
y = -8.10e-6 %gamma in 1/s
k_ak = 5.00e-7 % in 1/s
sv_BT = 20000 % s/v in brain tumour no unit
sv_NT = 7000 %s/v in normal tissue no unit
k_e_DOX = 5.8e-4 %elimination factor in 1/s
P_DOX = 3.00e-6 %transvascular permeability in tumour no unit
t_half_DOX = 1194.8 %half life in s
k_clr_DOX = 5.801e-4 %clearance in 1/s
Vt = 4.19 %tumour volume in cm2
Ve = 1.4665 %tumour extracellular from table in cm2
Tiv_DOX = 7200 %2 hours in s
D_DOX = 204 %dose of DOX in mg/m2
Vdist_DOX = 2.125 %L using 25L/kg for 85kg man?????
R_DOX = 10
 
%define the parameters for BEV
k_e_BEV = 1.2e-5 %elimination factor in 1/s
t_half_BEV = 57750 %in s
k_clr_BEV = (log(2))/t_half_BEV
Tiv_BEV = 7200 %2 hours in s
D_BEV = 850 %dose of BEV on average about 10mg/kg bu
Vdist_BEV = 3.29
R_BEV = 5
 
%define ODEs
%DOX conc in blood
if t<Tiv 
    diffeqs(1,1) = (-k_clr_DOX*C_DOX_Bl)+ R_DOX
else 
    R_DOX = 1 
else

Illegal use of reserved keyword "else".

Error in connector.internal.fevalMatlab

Error in connector.internal.fevalJSON

    R_DOX = 0
end
 
%BEV conc in blood
if t<Tiv
    diffeqs(2,1) = (-k_clr_BEV*C_BEV_Bl)+ R_BEV
else
    R_BEV = 1 
else
    R_BEV = 0
end
 
 
diffeqs(3,1) = ((P_DOX*phi*sv_BT*Vt)*(C_DOX_Bl - C_DOX_T) - (k_e_DOX*C_DOX_T*Ve))/Ve
diffeqs(4,1) = ((phi*(a+(B*phi)+(y*phi^2)))-(phi*k_ak*Cs_BEV_Bl))
diffeqs(5,1) = C_BEV_Bl/(D_BEV/Vdist_BEV)
end

[Walter Roberson]

What is your intention for

if something
else
    code
else
    more code
end

??

[Skye Cameron]

It is to find drug concentration in blood over a 24 hour period, I was advised to use an if loop, with the equation (-k_clr_DOX*C_DOX_Bl)+ R_DOX to help control thr 'infusion'

[Walter Roberson]

Suppose I wrote

if a < 5
   disp('outcome 1')
else
    disp('outcome 2')
else
    disp('outcome 3')
end

then under what circumstances would you expect outcome 3 to be what is displayed?

[Walter Roberson]

to help control thr 'infusion'

If you are adding more drug within the time span of any one ode45() call, then you will have problems, unless you can match derivatives. The easiest approach is to have any one ode45() call cover only the intervals between injections. For example if you inject every 6 hours then

curvar = var0;
hours_6 = 6*60*60;
for K = 1 : 4
  [t{K},var{K}] = ode45(@ode_sys, (K-1:K)*hours_6, curvar);
  curvar = var{K}(end,:) + amounts_to_inject;
end