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?

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

Answer 1

[Pratyush Roy]

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:

https://www.mathworks.com/matlabcentral/answers/26743-absolute-and-relative-tolerance-definitions

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!

[ABHIRUP SAMADDAR]

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