MATLAB was running for a long time before warning me there is an error

I was trying to solve an integrodifferential equaiton using numerical approach. My thought is that, I should use symbolic to evaluate the integral first, and then convert it into a double array to solve the equation.

The followings are my code:

%% EQUATIONS GOVERNING THE BUBBLE DYNAMICS
% Yihan Zhang
clc;clear;close all;
%% parameters
tau_y=8.6; % Pa
gamma=0.07; % Pa m
rho=1000; % kg/m^3
G=81.5; % Pa
R0=100e-6; % m
P0=1.13e5; % Pa
eta_s=0.001; % Pa s
w=((3*P0+4*gamma/R0+4*G)/(rho*R0^2))^0.5;
%% deltaP=0.1kPa
deltaP=100; % Pa
[t,R]=ode45(@BD,[0,0.001],[100e-6;0;deltaP]);
figure(1)
plot(2*pi.*t.*w,R(:,1)./R0)
xlabel('2 \pi t \omega');ylabel('R/R_0');
hold on
%% deltaP=0.2kPa
deltaP=200; % Pa
[t,R]=ode45(@BD,[0,0.001],[100e-6;0;deltaP]);
figure(1)
plot(2*pi.*t.*w,R(:,1)./R0)
xlabel('2 \pi t \omega');ylabel('R/R_0');
%% deltaP=2.5kPa
deltaP=2500; % Pa
[t,R]=ode45(@BD,[0,0.001],[100e-6;0;deltaP]);
figure(1)
plot(2*pi.*t.*w,R(:,1)./R0)
xlabel('2 \pi t \omega');ylabel('R/R_0');
axis([595,600,0.6,1.8]);
legend('\Delta p = 0.1 kPa','\Delta p = 0.2 kPa','\Delta p = 2.5 kPa')
%% EQUATIONS GOVERNING THE BUBBLE DYNAMICS
function dydt=BD(t,R)
%% parameters
gamma=0.07; % Pa m
rho=1000; % kg/m^3
G=81.5; % Pa
R0=100e-6; % m
P0=1.13e5; % Pa
deltaP=R(3);
eta_s=0.001; % Pa s
w=((3*P0+4*gamma/R0+4*G)/(rho*R0^2))^0.5;
%% try to evaluate the integral in integrodifferential equation
syms r
tau_rr=-4*G.*(R(1).^3-R0^3)./(3.*r.^3);
tau_tt=-tau_rr./2;
INT=(tau_rr-tau_tt)./r;
I=double(int(INT,R(1),inf));
%% other parameters
Pg=(P0+2*gamma/R0).*(R0./R(1)).^3;
P_inf=P0+deltaP.*sin(w.*t);
tau_rrR=double(subs(tau_rr,R(1)));
P=Pg+tau_rrR-eta_s.*4.*R(2)./R(1)-2.*gamma./R(1);
%% integrodifferential equation
dydt=[R(2);((P-P_inf-tau_rrR+2.*I)./rho-3/2.*R(2).^2)./R(1);0];
end

The result should be a diagram with three different curves.

At first, it runs successfully and gave me one diagram, but then it told me that there is an error shown by the error message below:

Error using symengine
Unable to convert expression into double array.

Error in sym/double (line 692)
        Xstr = mupadmex('symobj::double', S.s, 0);

Error in Bubble_dynamics5>BD (line 55)
P=Pg+tau_rrR-eta_s.*4.*R(2)./R(1)-2.*gamma./R(1);

Error in ode45 (line 299)
    f2 = odeFcn_main(t2, y2);

Error in Bubble_dynamics5 (line 24)
plot(2*pi.*t.*w,R(:,1)./R0)

I am quite confussed with this. Could someone help me please?

Sorry for my bad English.

Answer 1

[Sai Bhargav Avula]

Hi,

In the code you have written for deltaP = 200 is encountering the case of INF in the expression.

The issue can be because the integral that is being computed is not convergent enough. You can debug by using the command.

dbstop if naninf  

and then running your ode45 call. It will stop when the nan or infinity shows up and you can then examine the values of variables to see how it happened.