i am trying to plot bubble radius and step time for rayleigh plesset equation using euler method( ode45 for code) , but i am not getting the results as required.

function dydt = bubbleDynamics(t, y)
sigma=0.072;
mu=0.798e-3;
rho=996;
P0=1e-6;
Pv=0.00424;
Pg0=1e-5;
R0=1e-5;
R = y(1);
dRdt = y(2);
P= Pv + Pg0*(R0 / R)^3;
d2Rdt2 = ((P - P0 - (2*sigma / y(1))) /(y(1)*rho)) - 1.5*(y(2)^2)/y(1) - (4*mu*y(2)/rho*y(1)^2);
dydt = [dRdt; d2Rdt2];
end
tspan = 0:2e-9:0.00035; %  time range
R0=1e-5;
dR0=0;
sigma=0.072;
[t, y] = ode45(@bubbleDynamics, tspan, [R0, dR0]);
R = y(:, 1); % Extract bubble radius
plot(t, R);
xlabel('Time');
ylabel('Bubble Radius');
axis ([0 4e-5 0 4e-5])

Answer 1

[Walter Roberson]

In the below, notice that the first y value is about 6e-10 by iteration #426 -- but you are dividing by y(1) and y(1)^2 in your calculation, so that is giving large outputs.

Your y(1) appears to related to radius, so it appears to me that this is saying that as the radius gets smaller that the collapse rate accelerates dramatically.

That sounds like plausible bubble dynamics to me. Collapsing bubbles are known to result in cavitation and flashes of light, and low-order effectis might possibly become very important.

format long g
tspan = 0:2e-9:0.00035; %  time range
R0=1e-5;
dR0=0;
sigma=0.072;
[t, y] = ode45(@bubbleDynamics, tspan, [R0, dR0]);

dydt values got large by iteration #426

dydt = 2×1

1.0e+00 * -228957.899481437 -1.18518514742992e+20

y = 2×1

1.0e+00 * 6.63463960113321e-10 -228957.899481437

Warning: Failure at t=2.298637e-06. Unable to meet integration tolerances without reducing the step size below the smallest value allowed (6.776264e-21) at time t.

R = y(:, 1); % Extract bubble radius
plot(t, R);
xlabel('Time');
ylabel('Bubble Radius');
axis ([0 4e-5 0 4e-5])

function dydt = bubbleDynamics(t, y)
persistent counter triggered
if isempty(counter); counter = 0; triggered = false; end
counter = counter + 1;

sigma=0.072;
mu=0.798e-3;
rho=996;
P0=1e-6;
Pv=0.00424;
Pg0=1e-5;
R0=1e-5;
R = y(1);
dRdt = y(2);
P= Pv + Pg0*(R0 / R)^3;
d2Rdt2 = ((P - P0 - (2*sigma / y(1))) /(y(1)*rho)) - 1.5*(y(2)^2)/y(1) - (4*mu*y(2)/rho*y(1)^2);
dydt = [dRdt; d2Rdt2];

if ~triggered && any(abs(dydt) > 1e20)
    triggered = true;
    fprintf('dydt values got large by iteration #%d\n', counter);
    dydt
    y
end
    
end

[Walter Roberson]

I ran the calculation using symbolic numbers, in case there was a lot of cancelation handling or something like that. The 1e20 boundary was hit slightly earlier, at 709 iterations. This tells us that the problem is likely inherent to the equations as coded, rather than being a problem of limited precision.

Answer 2

[Rifat Shahrear]

I have found this result after taking assistantance from the coding that was previously given by others. However, can anyone explain why does the buddle radius become infinite after being zero? is there any problem?

[Torsten]

However, can anyone explain why does the buddle radius become infinite after being zero? is there any problem?

Sure. You divide by the radius when you compute d2Rdt2.