I am running ode15s solver but I get an error " Failure at t=1.576194e-08. Unable to meet integration tolerances without reducing the step size " How can I solve this error?

6 views (last 30 days)
ABHIRUP SAMADDAR
ABHIRUP SAMADDAR on 14 Jan 2022
Commented: ABHIRUP SAMADDAR on 17 Jan 2022
Following is my code and the function file named as R_of_tau called inside the script-file as R_of_tau
clear all;
tic
format long e
global p0 c cigma kapa ow mu PA0 PApulse0 tt PAWave R0 ii pg0 R_eqb n0 n fre ttmax rho maxtt prob_t_end delta_n dtt
%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%
%%%%%%%%%%%%%%%%%%%%%%%% CONSTANTS AND VALUES %%%%%%%%%%%%%%%%%%%%%%%%%%%%%
%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%
radius=380e-6;% size of the absorber (radius)
Rinitial=100e-9;
NR=5000;
cspeed=1500; %sound speed
Rstep=radius/NR; %the total step for time integral over one radial distance
delta_n=0;
fre=5e5; %Frequency of ultrasound
Tb=1/fre; %Time period for bubble oscillation
PA0=2e3; %in Pa for sine wave amplitude HIFU
PApulse0=6.98e6; %in Pa for peak negative PA pulse
ow=2*pi*fre; %Omega (Angular Velocity) of ultrasound
mua=25000; %optical absoprtion in /m
cspeed=1500; %sound speed
Fluence=20; %fluence in mJ/cm^2
H0=Fluence*10; %convert fluence to J/m^2
Gama=0.2; %Gruneisen constant
NR=5000;
Rstep=radius/NR; %the total step for time integral over one radial distance
dt=Rstep/cspeed;
p0=1.01e5; %Pressure at Infinity %%%%%%%%%%%%%%%%%%%%%%%%%%%%% why not blood pressure?
c=1500; %Speed of Sound
cigma=0.0725; %Surface Tension Coefficient (N/m)
% radius=150e-6 %Size of the blood vessel
kapa=1.4; %Polytropic Index for the gas inside the bubble
rho=1000; %Density of liquid (kg/m^3)
mu=0.005; %Viscosity
Rg=8.3145; %Universal Gas Constant (J/mol K)
M=28.966e-3; %Kg/mol (28.966 kg/kmol) Molecular Weight of the Gas Inside the Bubble
Tinf=293.15; %Fluid Temperature in Infinity in Kelvin
g=9.81; %acceleration due to gravity%
h=radius - Rinitial; %Depth of bubble from the blood vessel wall
ph=rho*g*h; % Hydrostatic pressure on the bubble surface
Henry_ref=1.3e-3; % Henry constant for oxygen in water at 298.15 K
Tref=298.15; % Reference temperature in K
delta_sol=1700;
H = Henry_ref*exp(delta_sol*(1/Tinf - 1/Tref)); % Henry's constant at Tinf
c_0 = p0/H; %%%%%%%%%%%%%%%%%% Saturation Gas Concentration in the Liquid mol/m^3
c_i=1.5*c_0;
Diffusivity=2e-9; %Diffusivity Constant in m^2/s
%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%
%%%%%%%%%%%%%% SETTING UP PRESSURE SOLVER AND MAIN FOR LOOP %%%%%%%%%%%%%%%
%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%
ndelay=400;
A=load('PA3ns760um.mat');
A.Finalpressure = A.Finalpressure(4500:19499);
A.tt0=A.tt0';
A.tt0=A.tt0(1:15000);
tmax=max(A.tt0);
NN = length(A.tt0);
A.tt0 = [A.tt0;(tmax+dt:dt:(tmax+ndelay*dt))'];
A.Finalpressure = [zeros(ndelay,1);A.Finalpressure];
tt=A.tt0; %A.ttnewnew; %Time increments for PA
maxtt=max(tt);
dtt=tt(2)-tt(1);
PAWave= A.Finalpressure; %A.Finalpressurenewnew; %PA data points based on tt
%plot(tt,PAWave)
lk=10; %Number of cycles to be calculated
ttmax=2*pi*lk/ow; %Maximum value for time
options1=odeset('RelTol',1e-6,'AbsTol',1e-7); %Options for the solver of keller
%Expanding time and PA points to have tmax seconds
%deltatt=min(diff(tt));
%tt=(0:deltatt:ttmax)'; %new tt replaced
%PAWave(length(PAWave)+1:length(tt))=0; % PA points modified to account for new tt
deltatt=2e-8;
ttadd=(maxtt+deltatt):deltatt:ttmax;
tt=[tt;ttadd']'; %new tt replaced
PAWave(length(PAWave)+1:length(tt))=0; % PA points modified to account for new tt
% figure (1e9+1)
% plot(tt,PAWave,'LineWidth',2);
% title('Convoluted Pressure from Paltauf with tstart=0 in expaned time','FontSize',16);
% xlabel('time','FontSize',14,'FontWeight','bold');
% ylabel('Convoluted Pressure','FontSize',14,'FontWeight','bold');
%saveas(figure(9e12+1),'ConvPressurePaltaufTstart0.png')
prob_t_domain=[linspace(0,maxtt,100),linspace((maxtt+dtt),ttmax,lk)]; % !!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!
R0=zeros(length(prob_t_domain),1);
R_eqb=zeros(length(prob_t_domain),1);
Rdot0=zeros(length(prob_t_domain),1);
R_eqb(1)=Rinitial; %Initial equilibrium bubble radius
R0(1)=Rinitial; %Initial bubble radius for calculations
Rdot0(1)=0; %Initial Velocity of bubble Radius
% to calculate n0 we use ideal gas law
n0=(p0+2*cigma/R_eqb(1))*(4/3*pi*R_eqb(1)^3)/(Rg*Tinf);
n=n0;
% to calculate the initial gas concentration
tau_end=0;
tau_prime_last=0;
DataTime=[];
DataRR0=[];
DataReqbR0eqb=[];
for ii=1:1:length(prob_t_domain)-1
% disp(['starting for loop for iteration ' num2str(ii)]);
prob_t_start=prob_t_domain(ii);
prob_t_end=prob_t_domain(ii+1);
%for ll=1:1
pg0=p0+2*cigma/R_eqb(1);
y0=[R0(ii) Rdot0(ii)]; % y0 has initial values for R and R'
[t,y]=ode15s(@SYNCHRONIZE,[prob_t_start prob_t_end],y0,options1);
% y(:,1) gives R and y(:,2) gives R'
tau_start=tau_end; % it equals the last tau_end for the new loop
[R_tau,tau_end] = R_of_tau(t,y(:,1),tau_start);
% R in R_tau is now a function of tau instead of t
syms c_s(R)
c_s=c_0*pg0/p0*(R_eqb(1)/R)^(3*kapa)*(n/n0)*(R_eqb(ii)/R_eqb(1))^(3*(kapa-1));
F(R)=c_s-c_i;%(10)
% Notice that tau is always between 0 and tau_end for solving.
tau_prime_begin = tau_prime_last;
%% Creation of tau_prime array %%
num_of_tau_prime=50;
tau_prime=linspace(tau_start,tau_end,num_of_tau_prime);
tau_prime_last=tau_prime(end);
delta_tau_prime=tau_prime(2)-tau_prime(1);
%% CREATION OF tau_double_prime array %%
num_of_tau_double_prime=50;
tau_double_prime=linspace(tau_prime_begin,tau_prime_last,num_of_tau_double_prime);
delta_tau_double_prime=tau_double_prime(2)-tau_double_prime(1);
%% CREATION OF MY FUNCTIONS TO BE differentiated within the integral %%
z = ((tau_prime_last-tau_double_prime).^0.5)./(tau_end-tau_double_prime);
%% Numerical differentiation of z with respect to tau_double_prime %%
diff_z=[];
for ss = 2:num_of_tau_double_prime-1
cen_diff=(z(ss+1)-z(ss-1))/...
(tau_double_prime(ss+1)-tau_double_prime(ss-1));
diff_z=[diff_z cen_diff];
end
%% Second Integration in Eller and Flynn paper %%
inT2=0;
for jj=1:num_of_tau_double_prime-4
inT2=inT2+(delta_tau_double_prime/2)*( (double(F(R_tau(tau_double_prime(jj+1)))*diff_z(jj+1))...
+ double(F(R_tau(tau_double_prime(jj)))*diff_z(jj))) );
end
%% First Integration in Eller and Flynn paper %%
inT1=0;
for mm=1:num_of_tau_prime-2
inT1=inT1+(delta_tau_prime/2)*( (tau_end-tau_prime(mm+1))^0.5/(R_tau(tau_prime(mm+1)))^3 +...
(tau_end-tau_prime(mm))^0.5/(R_tau(tau_prime(mm)))^3 );
end
%% production of integrals %%
inT_prod=inT2*inT1;
%% CALCULATION OF CHANGE IN MOLES BY FIRST ORDER SOLUTION %%
delta_n=32*Diffusivity*inT_prod;
n=n+delta_n;
coef_to_solve_for_R_eqb=[p0 2*cigma 0 -3*n*Rg*Tinf/(4*pi)];
R_eqb_new=roots(coef_to_solve_for_R_eqb);%(12)
R_eqb_new(find(imag(R_eqb_new)))=0;
%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%
for aaa=1:1:length(R_eqb_new) %
if R_eqb_new(aaa)>0; %
R_eqb(ii)=R_eqb_new(aaa); % = % R_eqb(ii)=real(R_eqb_new(real(R_eqb_new)>0));
end %
end %
%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%
R_eqb(ii+1)=R_eqb(ii);
DataTime=[DataTime; t];
DataRR0=[DataRR0; y(:,1)/R0(1)];
DataReqbR0eqb=[DataReqbR0eqb; R_eqb(ii)/R_eqb(1)];
R0(ii+1)=y(end,1);
Rdot0(ii+1)=y(end,2);
end
save 100nm-us1pa0.24-2.mat
plot(DataTime*1e3,DataRR0)
xlabel('Time(ms)')
ylabel('R/R0')
xlim([0 0.03])
% ylim([0 1000])
% xlim([9.96 10])
% xlim([0 0.08])
set(gca,'FontSize',14)
title('9.97-10 ms')
% set(gcf,'color','none'); %????????
% set(gca,'color','none');
% title('100 nm, Ultrasound only (0.82 MPa)')
title('100 nm, Ultrasound(0.82 MPa)+ Laser (PA0 = 0.24 MPa)')
%??????????????????????????????????
% figure (2)
% plot(DataTime,DataRR0)
% %title('PA = 1MPa US = 1MPa for 80 nm bubble')
% title('US = 1MPa for 80 nm bubble')
% xlabel('Time(s)')
% ylabel('R/R0')
% xlim([0 1e-5])
toc
function [R_tau_func,tau_end] = R_of_tau(t,R,tau_start)
% This function maps t to tau. first we do a change of reference in time to
% change t such that it starts at zero and not at t1
% so (t1,t2) changes to (0,t2-t1). In this way integration and solving for
% tau becomes easier.
R_func=fit(t,R,'linearinterp');
R_func2=fit(t,R.^4,'linearinterp');
time=(linspace(t(1),t(end),500))';
tau=zeros(length(time),1);
Rtau=zeros(length(time),1);
for ii=1:length(time)
tau(ii)=integrate(R_func2,time(ii),t(1))+tau_start;
Rtau(ii)=R_func(time(ii));
end
R_tau_func=fit(tau,Rtau,'linearinterp');
%q is now a function of tau
tau_end=tau(end);
end

Answers (1)

Pratyush Roy
Pratyush Roy on 17 Jan 2022
Hi Abhirup,
It is useful to set the absolute and relative tolerances to a higher value to avoid this warning. In case of sharp discontinuities in the reference input such as when tracking a reference signal, it is possible that this warning is generated at intermediate time steps.
MATLAB is trying to reduce the time step to a really small amount in order to counter the abrupt change due to the discontinuity in the reference signal.
For sharp discontinuities,it might not be possible to avoid this warning. However for non discontinous inputs we can set relative and absolute tolerance to a higher number than the default setting.
For more information on relative and absolute tolerance please look at the link below:
We can set the tolerances to a higher value for example by using the following commands:
options = odeset('RelTol',1e-4,'AbsTol',1e-6);
Hope this helps!
  1 Comment
ABHIRUP SAMADDAR
ABHIRUP SAMADDAR on 17 Jan 2022
Dear Pratyush Roy,
Yes, I increased my relative and absolute tolerance values as 1e-4,and 1e-6 respectively. It worked for me. Thank you for helping me to fix out this error.
Sincerely,
Abhirup Samaddar

Sign in to comment.

Community Treasure Hunt

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

Start Hunting!