# ODE45: doubts about the result. Correct or not?

I have to solve $$\left\{\begin{matrix} x\ddot{x}+\frac{3}{2}\dot{x}^2=\frac{1}{\rho}\left (P_v-a_1t^2-a_2t-P_0 \right )\\ x(0)=5\cdot10^{-6}\\ \dot{x}(0)=0 \end{matrix}\right.$$ I do not know if it may be helpful, but note that $a_1=2.488125\cdot 10^{14}$ is very big while the timespan is very small ($t_f=70\cdot 10^{-6}$). I get the result without MATLAB returning any error or warning, but this result does not convince me (I expected the peak to be around $t=4\cdot 10^{-5}$). I have already tried switching to ode15s, but the output is the same. What else can I do to see if the output changes or not and then be sure that this result is correct?

function vdot = no_thermal_effect_98(t,v)
rho = 959.78;
P0 = 101325;
Pvap = 94285.313;
a1 = (P0 - 1800) / (20e-6)^2;    %nondimensional
a2 = 2 * (1800 - P0) / (20e-6);  %nondimensional
%
vdot = zeros(2,1);
vdot(1) = v(2);
vdot(2) = -1.5 * v(2) * v(2) / v(1) + 1 / (v(1) * rho) *...
(Pvap - P0 - a1 * t^2 - a2 * t);

run

x0 = 5e-6; %meters
tf = 70e-6; %seconds
[t,v] = ode45(@no_thermal_effect_98,[0,tf],[x0,0]);
[t,v(:,1)];
plot(t,v(:,1))

Edit

Let us first consider the case where the pressure outside the bubble is constant: $x\ddot{x}+\frac{3}{2}\dot{x}^2=\frac{1}{\rho}\left (P_v-P_0\right )$.

If $P_0>P_v$ the bubble collapses (i.e. the radius decreases in a monotonic fashion); conversely, if $P_0<P_v$ the bubble grows.

Now let us consider a pressure variation outside the bubble described by the equation $P(t)=a_1t^2+a_2t+P_0$.

At $40 \mu s$ the pressure outside the bubble is greater than the pressure inside the bubble (here in the picture, indicated by the orange line and equal to $Pv=94285.313 Pa$). That is why I expected the peak, that is, the moment where the bubble started collpasing, would be $40 \mu s$.

• You've used the default (not very accurate) tolerance in ode45(). What happens when you tighten the error tolerance? – Brian Borchers Nov 28 '16 at 16:16
• @BrianBorchers I tried options = odeset('RelTol',1e-13,'AbsTol',1e-13); I got the same output. – Marco Nov 28 '16 at 16:20
• Why do you expect the output to be different from what you've got? – Brian Borchers Nov 28 '16 at 17:03
• Sundials also gives the same plot. I don't think this has anything to do with ode45. – Kirill Nov 28 '16 at 20:09
• I think the comparison with the constant-pressure case is wrong: if $\dot x$ is positive and $\ddot x$ given by the r.h.s. starts decreasing ($\ddot x<0$ for $t>24\mathrm{\mu s}$ for this ODE), it will take some time before $\dot x = 0$, there's no reason it would happen immediately. I think a better comparison might be with a piecewise-constant example, $P_v(t) = \mathrm{const} > P_0$ for $t\in[0,t_1)$ and $P_v(t) = \mathrm{const} < P_0$ for $t\in [t_1,t_f]$. – Kirill Nov 28 '16 at 21:06

I ran your same code in Python

from scipy.integrate import odeint
import numpy as np
import matplotlib.pyplot as plt

def func(v, t):
rho = 959.78
P0 = 101325
Pvap = 94285.313
a1 = (P0 - 1800) / (20e-6)**2
a2 = 2 * (1800 - P0) / (20e-6)
v1 = v[1];
v2 = -1.5 * v[1] * v[1] / v[0] + 1 / (v[0] * rho) *\
(Pvap - P0 - a1 * t**2 - a2 * t)
return [v1, v2]

y0 = 5e-6
tf = 70e-6
t = np.linspace(0, tf, 1000)
y = odeint(func, [y0, 0], t)
plt.plot(1e5*t, 1e5*y[:, 0])

and I obtained the same result

If you are worried about the numerical values of your constants, maybe you can change the units or non-dimensionalize your equation. For example you can change the time for $\tau = t/t_0$, with $t_0 = 20\times 10^{-6}$. And use the chain rule to make the changes to the derivatives in the left—hand side to obtain something like

$$y y''+\frac{3}{2}y'^2 = \frac{1}{\rho}\left [(P_v - P_0) - (P_0 - P_1)\tau^2 - 2(P_1 - P_0)\tau\right]$$

with $y = x/t_0$, $y' = dy/d\tau$, and $P_1 = 1800$. I tried this version and gave me the same answer.

Notes

• In both cases, I scaled the axes by $10^5$.
• It might be better to non-dimensionalize your equation using some knowledge of the physics in your problem. You can check reference 1 to get some ideas.

## References

1. Kudryashov, Nikolay A., and Dmitry I. Sinelshchikov. "Analytical solutions of the Rayleigh equation for empty and gas-filled bubble." Journal of Physics A: Mathematical and Theoretical 47.40 (2014): 405202.