Unable to meet integration tolerances without reducing the step size below the smallest value allowed (1.136868e-13) at time t.

3 views (last 30 days)
for tf=1:1:8
tspan = t0:step:tf*dose_length;
lt=length(tspan);
sx=size(xx1,1);
x0=xx1(sx,:);
% t0=t0+dose_length;
options = odeset('RelTol',10,'AbsTol',[10 10 10 10 ]);
[tt,xx1]= ode15s (@odebladder_exp,tspan,x0,options,muB,muE,alpha,pA,pT,pR,pK,pD,gRate,0, 0,K,g_T);
t0=t0+dose_length;
sx=size(xx1,1);
x1((tf-1)*lt+1:lt*tf,:)=xx1;
t1((tf-1)*lt+1:tf*lt)=tt;
end
  5 Comments
Sveta
Sveta on 4 Dec 2023
Edited: Torsten on 4 Dec 2023
thanks a lot for your help!
the program works when I change only one parameter in ParamBCG file:alpha = 0.7*10^(-4) => alpha = 3.7*10^(-4)
Wait for your answer.
The best,
Svetlana.
% model for Bladder :
echo off
close all
clear
% 1 - bacilla, 2 - immune respfonse, 3 - tumor with bacilla, 4 - tumor without bacilla
% mu1 - death of bacilla
% mu2 - death of factor
% mu3 - death of cells of immune response (ther experiments show that the average lifetime of a lymphocyte is 24.25 days,
% giving a death rate of 1/24.25 = 0.412 Philis )
% beta - proportion rate of infection of tumor cells with bacilla
% par1 - death bacilla by macrophag
% par2 - tumor cells stimulate growth of effector cells
%
%global mu1 mu3 beta alpha par1 par2 par3 gRate b_incj count ;
ParamBCG;
MAINT = 1;
global muE muB;
global pK pR pA pT;
global r;
global alpha pD K g_T IFN
% TU_i = 7.8*10^6;%100;%2*10^7;
DOSE = 2.2*10^6;
t0=0;
xx1 = [0 20 0 8*10^6 ];
%x1=zeros(12100,4);
step = 0.1;
dose_length = 7;
xxt=xx1(4);
YINDEX = 500;
urot_State(1:YINDEX) = struct('growth_rate',0,'init_size',0,'result',zeros(2100),'t1',zeros(2100),'color','r', 'lnwdth',1.5);
sd=0;%stable desease
pr=0;nr=0;cr=0;
T0_ini= 10^6;%
T0_fin= 7.6*10^7;
T0_step= 10; %WAS 5; %determines how many elements in the vector
T0=linspace(T0_ini,T0_fin,T0_step);
DOSE = 2.2*10^6;%was 2.2
r0_ini= 0.016;%MICHAEL 0.008
r0_fin= 0.024;
% r0_ini= 0.0025;
% r0_fin= 0.0038;%%3.7.22- was 0.0015-0.0025 ; was 98
r0_step = 5; %WAS10;%detemines how many elelments in the vector
r0=linspace(r0_ini,r0_fin,r0_step);
kk=size(r0);
f_r=randperm(kk(2));
kk=size(T0);
f_T = randperm(kk(2));
k=1;
for i_T=1:T0_step;
for i_r=1:r0_step;
gNum = f_r(i_r);
gRate = r0(gNum);
t0=0; tf=7;
tspan = [t0 tf];
gNumT=f_T(i_T);
t_0 = T0(gNumT)
% xx1 = [b_incj APC_i AT_i AB_i EB_i ET_i IL2_i TI_i t_0];
xx1 = [0 20 0 t_0];
%before treatment
for tf=1:1:8
tspan = t0:step:tf*dose_length;
lt=length(tspan);
sx=size(xx1,1);
x0=xx1(sx,:);
% t0=t0+dose_length;
options = odeset('RelTol',0.001,'AbsTol',[0.001 0.001 0.001 0.001 ]);
[tt,xx1]= ode15s (@odebladder_exp,tspan,x0,options,muB,muE,alpha,pA,pT,pR,pK,pD,gRate,0, 0,K,g_T);
t0=t0+dose_length;
sx=size(xx1,1);
x1((tf-1)*lt+1:lt*tf,:)=xx1;
t1((tf-1)*lt+1:tf*lt)=tt;
end
b_incj = 2.2*10^6; %DOSE ;
for tf=9:1:14 %was 14
tspan = t0:step:tf*dose_length;
lt=length(tspan);
sx=size(xx1,1);
xx1(sx,1)=xx1(sx,1)+b_incj;
x1((tf-1)*lt+1:lt*tf,:)=xx1;
disp(size(xx1));
x0=xx1(sx,:);
disp(size(xx1));
%tic
[tt,xx1]= ode15s (@odebladder_exp,tspan,x0,options,muB,muE,alpha,pA,pT,pR,pK,pD,gRate,0, IFN,K,g_T);
%toc
t0=t0+dose_length;
disp(size(t1));
disp(size(tt));
t1((tf-1)*lt+1:tf*lt)=tt;
disp(size(xx1));
end
%without treatment
for tf=15:1:27
tspan = t0:step:tf*dose_length;
lt=length(tspan);
sx=size(xx1,1);
x0=xx1(sx,:);
%tic
[tt,xx1]= ode15s (@odebladder_exp,tspan,x0,options,muB,muE,alpha,pA,pT,pR,pK,pD,gRate,0, 0,K,g_T);
%toc
t0=t0+dose_length;
sx=size(xx1,1);
x1((tf-1)*lt+1:lt*tf,:)=xx1;
t1((tf-1)*lt+1:tf*lt)=tt;
end
%maintanance treatment
b_incj=0;
kmaint=0;
for tf=28:1:140
tspan = t0:step:tf*dose_length;
lt=length(tspan);
sx=size(xx1,1);
x0=xx1(sx,:);
%tic
[tt,xx1]= ode15s (@odebladder_exp,tspan,x0,options,muB,muE,alpha,pA,pT,pR,pK,pD,gRate,0, IFN,K,g_T);
%toc
t0=t0+dose_length;
sx=size(xx1,1);
if MAINT ==1
if (tf==28 || tf==40 || tf==64 || tf==88 ||tf==112 || tf==136 || kmaint>0 && kmaint<3)
% if (tf==27 || tf==42 || tf==60 || tf==76 || tf==94 || tf==108 || kmaint>0 && kmaint<3)
b_incj=DOSE ;
kmaint=kmaint+1;
else
b_incj = 0;
kmaint=0;
end
end
xx1(sx,1)=xx1(sx,1)+b_incj;
disp(size(xx1));
x1((tf-1)*lt+1:lt*tf,:)=xx1;
t1((tf-1)*lt+1:tf*lt)=tt;
end
%after treatment
b_incj = 0;
for tf=141:1:170
tspan = t0:step:tf*dose_length;
lt=length(tspan);
sx=size(xx1,1);
x0=xx1(sx,:);
[tt,xx1]= ode15s (@odebladder_exp,tspan,x0,options,muB,muE,alpha,pA,pT,pR,pK,pD,gRate,0, IFN,K,g_T);
t0=t0+dose_length;
sx=size(xx1,1);
x1((tf-1)*lt+1:lt*tf,:)=xx1;
t1((tf-1)*lt+1:tf*lt)=tt;
end
zz = 1:length(t1);
t1=t1(1:length(t1));
adTw1=x1(zz,4);
if adTw1(end)<0.1
cr=cr+1;%complete response
col='-b';
lnwdth = 2.3;
else
if adTw1(end)<=(0.7*adTw1(1))%0.5
lnwdth = 4;
col='--g';
pr=pr+1;%partial response
else
if (adTw1(end)>0.7*adTw1(1))&&(adTw1(end)<=1.2*adTw1(1))
sd=sd+1;%stable desease
col=':k';
lnwdth = 1.5;
else
col='-r';
lnwdth = 0.8;
nr=nr+1; %progressive desease
end
end
end
hold on
figure (1);
%subplot(2,1,2);
plot( t1,adTw1,col,'LineWidth',lnwdth);
urot_State(k).color = col;
urot_State(k).growth_rate=gRate;
urot_State(k).init_size=T0(gNumT);
urot_State(k).result=x1;
urot_State(k).t1=t1;
urot_State(k).lnwdth=lnwdth;
k=k+1;
end;
end
t_0 = 26000000
71 4
71 4
1 568
71 1
71 4
71 4
71 4
1 639
71 1
71 4
71 4
71 4
1 710
71 1
71 4
71 4
71 4
1 781
71 1
71 4
71 4
71 4
1 852
71 1
71 4
71 4
71 4
1 923
71 1
71 4
71 4 71 4 71 4 71 4 71 4 71 4 71 4 71 4 71 4 71 4 71 4 71 4 71 4 71 4 71 4 71 4 71 4 71 4 71 4 71 4 71 4 71 4 71 4 71 4 71 4 71 4 71 4 71 4 71 4 71 4 71 4 71 4 71 4 71 4 71 4 71 4 71 4 71 4 71 4 71 4 71 4 71 4 71 4 71 4 71 4 71 4 71 4 71 4 71 4 71 4 71 4 71 4 71 4 71 4 71 4 71 4 71 4 71 4 71 4 71 4 71 4 71 4 71 4 71 4 71 4 71 4 71 4 71 4 71 4 71 4 71 4 71 4 71 4 71 4 71 4 71 4 71 4 71 4 71 4 71 4 71 4 71 4 71 4 71 4 71 4 71 4 71 4 71 4 71 4 71 4 71 4 71 4 71 4 71 4 71 4 71 4 71 4 71 4 71 4 71 4 71 4 71 4 71 4 71 4 71 4 71 4 71 4 71 4 71 4 71 4 71 4 71 4 71 4
t_0 = 26000000
71 4
71 4
1 12070
71 1
71 4
71 4
71 4
1 12070
71 1
71 4
71 4
71 4
1 12070
71 1
71 4
71 4
71 4
1 12070
71 1
71 4
71 4
71 4
1 12070
71 1
71 4
71 4
71 4
1 12070
71 1
71 4
71 4 71 4 71 4 71 4 71 4 71 4 71 4 71 4 71 4 71 4 71 4 71 4 71 4 71 4 71 4 71 4 71 4 71 4 71 4 71 4 71 4 71 4 71 4 71 4 71 4 71 4 71 4 71 4 71 4 71 4 71 4 71 4 71 4 71 4 71 4 71 4 71 4 71 4 71 4 71 4 71 4 71 4 71 4 71 4 71 4 71 4 71 4 71 4 71 4 71 4 71 4 71 4 71 4 71 4 71 4 71 4 71 4 71 4 71 4 71 4 71 4 71 4 71 4 71 4 71 4 71 4 71 4 71 4 71 4 71 4 71 4 71 4 71 4 71 4 71 4 71 4 71 4 71 4 71 4 71 4 71 4 71 4 71 4 71 4 71 4 71 4 71 4 71 4 71 4 71 4 71 4 71 4 71 4 71 4 71 4 71 4 71 4 71 4 71 4 71 4 71 4 71 4 71 4 71 4 71 4 71 4 71 4 71 4 71 4 71 4 71 4 71 4 71 4
t_0 = 26000000
71 4
71 4
1 12070
71 1
71 4
71 4
71 4
1 12070
71 1
71 4
71 4
71 4
1 12070
71 1
71 4
71 4
71 4
1 12070
71 1
71 4
71 4
71 4
1 12070
71 1
71 4
71 4
71 4
1 12070
71 1
71 4
71 4 71 4 71 4 71 4 71 4 71 4 71 4 71 4 71 4 71 4 71 4 71 4 71 4 71 4 71 4 71 4 71 4 71 4 71 4 71 4 71 4 71 4 71 4 71 4 71 4 71 4 71 4 71 4 71 4 71 4 71 4 71 4 71 4 71 4 71 4 71 4 71 4 71 4 71 4 71 4 71 4 71 4 71 4 71 4 71 4 71 4 71 4 71 4 71 4 71 4 71 4 71 4 71 4 71 4 71 4 71 4 71 4 71 4 71 4 71 4 71 4
Warning: Failure at t=6.160957e+02. Unable to meet integration tolerances without reducing the step size below the smallest value allowed (1.818989e-12) at time t.
1 4
Unable to perform assignment because the size of the left side is 71-by-4 and the size of the right side is 1-by-4.
text(350,120e7,[' CR=',num2str(cr)],'Color','b','FontSize',10);
text(350,100e7,[' PR=',num2str(pr)],'Color','g','FontSize',10);
text(350,80e7,[' SD=',num2str(sd)],'Color','k','FontSize',10);
text(350,60e7,[' NR=',num2str(nr)],'Color','r','FontSize',10);
disp([cr ,sd , pr , nr ])
xlabel('days');
ylabel('tumor cells');
title(' BCG induction&maintenance for 50 patients')
text (-45, 23*10^8, 'B ');
% num_k=randperm(500);
% figure(2);
% hold on
% for i_T=1:100;%length(num_k)/2;
% k=num_k(i_T);
% t1=urot_State(k).t1;
% x1=urot_State(k).result;
% col=char(urot_State(k).color);
% t1=t1(1:length(t1));
% zz = 1:length(t1);
% lnwdth = urot_State(k).lnwdth;
%
% adTw1=x1(zz,9);
% plot( t1,adTw1,col,'LineWidth',lnwdth);
% end;
% hold off
% zz = 1:length(t1);
% t1=t1(1:length(t1));
% adB=x1(zz,1);
% adE=x1(zz,2);
% adTb=x1(zz,3);
% adTw=x1(zz,4);
%
% %%plotyy
% hold on
% figure (1);
% subplot(2,1,1);
% %plot( t1,adTw,'-r','LineWidth',3.1);
% % plot( t1,adB,'-k','LineWidth',1.9);
%
% plot( t1,adB,'--k',t1,adTw,'-r','LineWidth',2.1);
% legend('BCG','T');
%
% subplot(2,1,2);
% %plot( t1,adTw,'-r','LineWidth',3.1);
% plot( t1,adB,'-k','LineWidth',1.9);
%
% xlabel('time(days)');
% ylabel('number of cells ');
function xdot = odebladder_exp(t,x,muB,muE,alpha,pA,pT,pR,pK,pD,r,b_injc,IFN,K,g_T)
% ode23s ode45- good for eps = 0
%global mu1 mu3 beta alpha par1 par2 par3 gRate b_incj count ;
%%%%xdot = zeros (4,1);
% 1 - bacilla, 2 - immune response, 3 - tumor with bacilla, 4 - tumor without bacilla
% this time we try without INF
%n=3;
if IFN == 1
Q1= 10; Q2 = 7.82;% Q1= 30; Q2 = 21.82; %Q1=3*10; Q2 = 3*7.3;
else
Q1=1; Q2 = 1;
end
xdot(1) = b_injc+ x(1)*(- pT*x(4) - pA*x(2) -muB) ;%- pT*x(4)
xdot(2) = ((pR*x(1)*x(2))+ (alpha-pD)*x(3)*x(2)+alpha*x(2)*x(3))*Q2 -x(2)*muE;
xdot(3) = pT*x(1)*x(4) - Q1*pK*x(2)*x(3);
%xdot(4) = r*x(4)-pT*x(1)*(1-x(4)/K) - Q1*pK*x(4)*x(2)/(1+5.2*10^3);% WORK!
xdot(4) = r*x(4)*(1-x(4)/K)-pT*x(1)*x(4) - Q1*pK*x(4)*x(2)/(1+g_T);%;*x(4)RIGHT
%xdot(4) = r*x(4)*(1-x(4)/K) - Q1*pK*x(4)*x(2)/(1+5.2*10^3) ;%;*x(4)-pT*x(1)*(1-x(4)/K)
xdot=xdot';
end
function ParamBCG()
global muE muB;
muE = 0.19;
% muE1=0.19;muE2=0.034;
muB = 0.1;
global pA pT pD pK pR %pE = pR
%pA= 2.16*10^(-6); this value do not allow eradication of the tumor
pA= 0.1925*2.16*10^(-6);
pT= 0.028*10^(-2);%0.028*10^(-3);%MICHAEL
pD = 1.03*10^(-10);
pK=0.4*10^(-6);%10^(-6);%(-7);
pR = 0.223e-8;%WAS0.223e-8;
global g_T r;
r = 0.0015;%0.0038;
g_T=5.2*10^3;
global alpha b_injc
alpha = 0.7*10^(-4);%3.7*10^(-6);Michael
b_injc = 2.2e6;
global K;
K = 10^(12);
global IFN;
global MAINT;
IFN = 0;%1;
MAINT = 0;
end
Sam Chak
Sam Chak on 4 Dec 2023
@Sveta, I just want to check the stability of the original system with some random initial conditions. So far, I haven't encountered any integration failure messages in this relatively simple code. I prefer to include the fixed parameters inside the ODE function for easy reference. Which parameters do you intend to change during the integration?
tspan = [0 1000];
x0 = rand()*[1e7 2e4 3e4 4e4];
[t, x] = ode45(@odebladder_exp, tspan, x0);
plot(t, x), grid on
function xdot = odebladder_exp(t, x)
xdot = zeros (4, 1);
%% parameters (some are unused)
muB = 0.1;
muE = 0.19;
alpha = 0.7e-4;
pA = 0.1925*2.16e-6;
pT = 0.028e-2;
pR = 0.223e-8;
pK = 0.4e-6;
pD = 1.03e-10;
r = 0.0015;
b_injc = 2.2e6;
IFN = 0;
K = 1e12;
g_T = 5.2e3;
Q1 = 1;
Q2 = 1;
%% differential equations
xdot(1) = b_injc + x(1)*(- pT*x(4) - pA*x(2) - muB);
xdot(2) = (pR*x(1)*x(2) + (alpha - pD)*x(3)*x(2) + alpha*x(2)*x(3))*Q2 - x(2)*muE;
xdot(3) = pT*x(1)*x(4) - Q1*pK*x(2)*x(3);
xdot(4) = r*x(4)*(1 - x(4)/K) - pT*x(1)*x(4) - Q1*pK*x(4)*x(2)/(1 + g_T);
end

Sign in to comment.

Answers (1)

Torsten
Torsten on 4 Dec 2023
Edited: Torsten on 4 Dec 2023
I had to add the line
global alpha pD K g_T IFN
at the beginning of your program to make the code work because the above global variables were not available.
In some phase of your code run, the integrator is no longer able to integrate your system of differential equations and gives up with the error message
Warning: Failure at t=6.160957e+02. Unable to meet integration tolerances without reducing the step size below the smallest value allowed (1.818989e-12) at time t.
After each call of the integrator, you expect an array of size (71x4) as result and you want to save it into a different array. But since the integrator stops already at time 616,.. , the (71x4) array is not available and you make an assignment that would only have been possible if the integrator had succeeded. This gives the error message
Unable to perform assignment because the size of the left side is 71-by-4 and the size of the right side is 1-by-4.
and MATLAB stops.
I cannot help you in this respect since the error stems from the integration of your model equation that I don't know.
  1 Comment
Sveta
Sveta on 5 Dec 2023
thanks for your message. I sow matLAB has settled on integration. it is a problem of the function of Matlab - ode45 or ode15s. Maybe is it possible to change them to another function or to change set of parameters in options?
I try to solve it too much time and cannot to find something good...
Please, help me.

Sign in to comment.

Categories

Find more on Symbolic Math Toolbox in Help Center and File Exchange

Products


Release

R2021a

Community Treasure Hunt

Find the treasures in MATLAB Central and discover how the community can help you!

Start Hunting!