Unable to meet integration tolerances without reducing the step size below the smallest value allowed (1.136868e-13) at time t.
Mostra commenti meno recenti
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 Commenti
Sveta
il 3 Dic 2023
Torsten
il 3 Dic 2023
You only want a precision of 10 ? That means the true solution and the numerical solution are allowed to deviate by 10 ?
Further comments are not possible since we cannot execute your code (function odebladder_exp is missing).
Sveta
il 4 Dic 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
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
@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
Risposte (1)
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 Commento
Sveta
il 5 Dic 2023
Categorie
Scopri di più su Stochastic Differential Equation (SDE) Models in Centro assistenza e File Exchange
il 3 Dic 2023
il 5 Dic 2023
Community Treasure Hunt
Find the treasures in MATLAB Central and discover how the community can help you!
Start Hunting!
Translated by 
