Numerical integration of density hydrostatic balance equation

I am a phd student trying to reproduce similar results to those found in Mizuno, Nakazawa & Hayashi, 1978. The idea is integrating the following set of first order differential equations

$\frac{d \rho}{dr} = - \left(\frac{G\mu}{k_{b} T}\right)\frac{M_{r} \rho}{r^{2}}$

$\frac{d M_{r}}{dr} = - 4\pi r^{2}\rho$

Given initial conditions $M_{r}\left(r=r_{p}\right) = M_{0}$, $\rho\left(r=_{p}\right) = \rho_{0}$. My domain is the Hill Sphere of a planet immersed in a protoplanetary disk, from the planet radius ($r_{p}$) up to the Hill Sphere boundary. $G$ is the gravitation constant, $k_{b}$ is the boltzmann constant, $T$ is the temperature and $\mu$ is the mass of the average particle. In particular, $M_{0}$ equals the mass of a Jupiter type planet. I explore the other parameter $\rho_{0}$, the density at the planetary top boundary.

I have explored several numerical schemes which are kind of straightforward to me and probably well known to all of you; Explicit Euler, Implicit Euler and 4º order Runge Kutta.

Setting $M_{0}$ to its proper value seems to trigger an instability making the coefficients in the previous schemes to diverge at the few first steps. Increasing the domain discretization looked like a good solution to circle this issue but the level of discretization has to be higher than my computational capacities are able to put up with. Implicit Euler on the other side seems to be stable but yields a non physical result (oscillatory). If $M_{0}= 0$ the solutions are ok but do not represent what I want to model

I feel I am missing something when dealing with this system, and I do not know which other approach I may take to tackle this problem. The previous paper do not mention anything related to numerical schemes used. Do you have any suggestions?

Thank you very much.

• How are you discretizing the first equation? For example, using an explicit Euler, I get something like: $\rho^{i+1} = \rho^i - \Delta r_i^{i+1} \left(G\mu/k_bT\right) M_r^i \rho^i / (r^i)^2$ if I do it in a simple way. That won't work for updating $\rho^{i+1}$ because there is a division by $r^i = 0$ for the first point... If $M_r(0) = 0$ then your division by zero may just not matter in the floating point math (depending on your programming language). – tpg2114 May 31 '17 at 9:27
• Yes I exactly do that. But you need to take into account that I am starting at $r = r_{p}$ where $r_{p}$ is the planet radius. So that it is not a problem. By the way I use a uniform spacing, hence $\Delta r$ is constant. Thanks!. – Juan Luis Gómez González May 31 '17 at 9:47
• How are you taking into account the boundary conditions at $r=0$ if you don't start at/near $r=0$? And using a uniform grid is potentially okay, but your equation is very stiff for small values of $r$ so a globally uniform spacing that is fine enough to resolve the function may use too many points for your computer, while a reasonable number of points may be too coarse to resolve the gradients. You might want to look at non-uniform grid spacing to make sure you keep costs down while putting points where they are needed. Or, you can just use a pre-packaged ODE solver and not roll your own. – tpg2114 May 31 '17 at 9:58
• About taking initial conditions at $r=0$ it is a mistake I have made explaining my problem. I think now I have written it more clearly. I haven't gone through the idea of making a non uniform grid. Maybe I will try to include a spacing that somehow increases the density according to the gradient of $1/r^{2}$. Until you have mention the word "stiff" I was thinking that problem do not arise in this case. I am just working in the analytical solution arising from making $M_{r}$ constant, but it seems it does not work either, the system requires more gravity. Next I will implement your idea!. – Juan Luis Gómez González May 31 '17 at 10:13

$dr \leqslant \frac{r^{2}}{\frac{G \mu}{k_{b}T}}$
$\frac{1}{\rho}\frac{d \rho}{dr} = \frac{G \mu}{k_{b}T} \frac{1}{r^{2}}$
• Glad it worked! That looks something like a CFL number used when solving PDE's to restrict the step sizes based on the coefficients in the equation. You may find that you can take bigger $\Delta r$ with the RK4 scheme because that generally has an extended stability range relative to the Euler scheme (CFL = $2 \sqrt{2}$ vs. CFL = 1). The implicit Euler is usually A-stable for all step sizes, but it might have poor convergence and accuracy if the size is too big. – tpg2114 May 31 '17 at 16:29