Implicit formulation of finite difference with ode15s

I'm solving a set of 3 PDEs, 4 ODEs, and 1 AE in MatLab. The model describes the dynamic behavior of an argon/water vapor spherical bubble in water undergoing an oscillatory pressure field (P = P0+PA*sin(wt)). The first 2 PDEs describe the mass and energy conservation along the bubble radius [0,R], while the third PDE is the energy conservation in the liquid [R,+inf]:

The 4 ODEs are for the functions R(t), R'(t), P(t) and m(t), i.e. the time dependence of bubble radius, bubble wall speed, pressure in the bubble and total mass in the bubble. An additional AE is required to find the value of overdot m which is the mass flux at the interface. In order to solve the system, I've discretized the PDEs using a straightforward finite difference method with second-order accurate central differences for the spatial derivatives (properly modified at the boundaries to one-sided differences of the same accuracy order):

Therefore, I obtain a system of only ODEs and AEs that I solve defining a mass matrix and passing it to the ode15s solver. I am quite sure the code is properly written as it gives reasonable outputs for certain range of the parameters' value. However, for other ranges, the spatial derivatives at r=R of the liquid temperature PDE become very sharp and the code gives unrealistic outputs (or it just crashes before finishing the integration). In literature, it is advised to use an implicit formulation for that particular PDE. (For whoever works in a university there are good possibilities to have access to that paper at the following link PAPER. The following picture is an example of the spatial profile of the liquid temperature at different times obtained for certain parameters values.

If the acoustic pressure is increased above a certain threshold, the gradients of those profile at the bubble wall (r=R) should be much sharper (as the temperature goes up to several thousand Kelvins in few nanoseconds. That's when my code fails.

My question is whether it is possible to implement the implicit method within the tools I am using and, if yes, how to do it. I would like to keep using ode15s because the time dependence is also really steep and discretizing time as well would require writing a variable time step algorithm (which I am not familiar with).

Thanks in advance.

[Bill Greene]

You really haven't provided enough information for anyone to help you. I suggest including your equations and also a plot of the solution that shows the behavior you describe.

The term "implicit" is usually used to describe a general class of ODE solution methods that require the solution of a system of equations at each time step; the ode15s function uses an implicit algorithm.

[Davide Masiello]

Hi Bill, thanks for replying. I've updated my question as you suggested.

Answer 1

[Bill Greene]

I am still unclear on many aspects of your problem. And I am unable to access the paper you reference so that is of no help to me.

But here are a few comments (guesses really) based on what I can see from your equations. Your three PDE have both first and second order spatial derivatives. When the first order terms become large relative to the second, there are difficulties in obtaining a numerical solution. There are many references describing numerical solution of convection-diffusion equations you can consult for more information.

One simple approach to avoid numerical difficulties is to discretize the first order derivatives using a forward difference formula rather than the central difference you have used. This is often referred to as "upwinding" in the literature.

Finally, I am still unclear about what you (or the paper) mean by an "implicit method." But, as I tried to hint in my previous comment, if you are using ode15s for solution of the ODE system, I would say you already ARE using an implicit method.