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

Question

Sveta 2023년 12월 3일

0
링크

이 질문에 대한 바로 가기 링크

https://kr.mathworks.com/matlabcentral/answers/2055674-unable-to-meet-integration-tolerances-without-reducing-the-step-size-below-the-smallest-value-allow

댓글: Sveta 2023년 12월 5일

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
이전 댓글 3개 표시이전 댓글 3개 숨기기

Sveta 2023년 12월 4일

편집: Torsten 2023년 12월 4일

MATLAB Online에서 열기

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

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

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

1 568

71 1

71 4

1 639

71 1

71 4

1 710

71 1

71 4

1 781

71 1

71 4

1 852

71 1

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

1 12070

71 1

71 4

1 12070

71 1

71 4

1 12070

71 1

71 4

1 12070

71 1

71 4

1 12070

71 1

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

1 12070

71 1

71 4

1 12070

71 1

71 4

1 12070

71 1

71 4

1 12070

71 1

71 4

1 12070

71 1

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 2023년 12월 4일

MATLAB Online에서 열기

@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

댓글을 달려면 로그인하십시오.

이 질문에 답변하려면 로그인하십시오.

Answer 1

Torsten 2023년 12월 4일

0
링크

이 답변에 대한 바로 가기 링크

https://kr.mathworks.com/matlabcentral/answers/2055674-unable-to-meet-integration-tolerances-without-reducing-the-step-size-below-the-smallest-value-allow#answer_1365299

편집: Torsten 2023년 12월 4일

MATLAB Online에서 열기

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
이전 댓글 -1개 표시이전 댓글 -1개 숨기기

Sveta 2023년 12월 5일

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.

댓글을 달려면 로그인하십시오.

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

댓글 수: 5
이전 댓글 3개 표시이전 댓글 3개 숨기기

답변 (1개)

댓글 수: 1
이전 댓글 -1개 표시이전 댓글 -1개 숨기기

참고 항목

카테고리

태그

제품

릴리스

Community Treasure Hunt

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

댓글 수: 5 이전 댓글 3개 표시이전 댓글 3개 숨기기

답변 (1개)

댓글 수: 1 이전 댓글 -1개 표시이전 댓글 -1개 숨기기

참고 항목

카테고리

태그

제품

릴리스

Community Treasure Hunt

댓글 수: 5
이전 댓글 3개 표시이전 댓글 3개 숨기기

댓글 수: 1
이전 댓글 -1개 표시이전 댓글 -1개 숨기기