Dear Chris,

First of all, I would like to make two points:

- in my entire career as a lecturer and university professor (as well as in my life), I have never used a red herring. I wanted to emphasise this.
- The first time I solved the Rayleigh-Plesset equation was in 1989 with a stiff routine in Fortran. I have never doubted that ODE solvers are efficient. I am successfully using some of the ones you contributed to SciML.

Regarding the mathematical and physical aspects associated with the Rayleigh-Plesset equation, I am not the only one to have obtained results with a bubble radius reaching 0 or négative values which stop computation. (see figure below).

Reference
This happens when the acoustic pressure is too high or when the radius of a bubble is not adapted. And this has a physical meaning. Indeed, the Rayleigh-PLesset equation is only valid for spherical bubbles. However, towards the end of the implosion, instabilities appear (sometimes Rayleigh-Taylor instabilities, sometimes for very small bubbles: an inward jet). This can lead to a bubble fragmenting before reaching its minimum radius.

In the spherical approximation, I have of course tried Rosenbrock23, KenCarp4, DPRKN6() with different reltol and abstol. I can reproduce well the evolution of the radius as a function of time or the evolution of the temperature as a function of time. Here are 2 figures for a bubble of initial radius R0 = 2 µm in a water/glycerine mixture. Already at this stage, I am well aware that I must modify the rpnnp() function. Indeed, a bubble at R0 = 2 µm contains ~1e9 N2 molecules (or Ar atoms). If even we could reach the density of a liquid at the end of the collapse, the radius should not be less than 180 nm. However, in the case of figure R-t, it already reaches 45 nm…

Calculated with KenCarp4.

So, be sure I know that the rpnnp() function must be modified. And not only from this point of view, the T-t curve should make you tick. The expansion phase is slow enough (for low ultrasonic frequencies: 20 kHz - 40 kHz) for heat transfer to occur, while the collapse is fast enough to be considered adiabatic. Secondly, the evaporation of the liquid at the interface and its condensation during collapse must be taken into account. Not only that, but due to a temperature increase of the order of 3_000 to 5_000 K (spectroscopic data), there are radical chemical reactions (pyrolysis) and chemiluminescence reactions.

Coming back to the mathematical problem, you are perfectly right, of course (I have to transmit complex u0). You must be aware that sometimes I don’t manage to formalise code well without a concrete example.

Finally, I don’t want to leave the impression that I said that the rpnnp problem was unsolvable. It turns out that I neglected it for several months following work on ODEs (chemical kinetics, NMR, 1D Schröginger, etc.), following a chemical kinetics study with Catalyst.jl, translating Think Julia’s book into French and writing a document for (bio)chemists about Flux.jl.

So I’ll be coming back to rpnnp in a more relentless way. Please keep in mind that when I ask a question, I am mostly looking for concrete help in the form of an example. I try never to disturb the forum for trivia but only when I really don’t understand something.

Sincerely,

Thierry