# Solving differential equation for density in a gas sphere

I'm trying to derive the equation for density as a function of height in a gas sphere due to gravitational force, and I have derived the following equation: $$\frac{\mathrm{d}\rho}{\mathrm{d}z}=\frac{4 \pi G}{z^2RT}\rho(z)\int_{0}^{z}z'^2\rho(z')\mathrm{d}z'$$ Is there a way to solve it, even if numerically, for $$\rho(z)$$?

• Is $T$ constant or a function of $z$? Working on an idea – Stephen Donovan Apr 12 at 14:33
• $T$ is the the temperature of the gas. We can consider it to be constant, though it would be interesting to consider it a function of $z$. – ordptt Apr 12 at 14:36
• Either works, I guess I'll give the variable version and if you want the constant it shouldn't be difficult to obtain. I assumed $\rho$ can never be zero, is that okay? Might be a dumb question but I figured I should ask – Stephen Donovan Apr 12 at 14:38
• That is a good assumption. If there is gas, there must a non-zero density. – ordptt Apr 12 at 14:39
• Good, just wanted to make sure. Working on an answer now – Stephen Donovan Apr 12 at 14:42

My thought process for this is that I really wanted to get rid of that integral, so I started by isolating it on the right-hand side:

$$\dfrac{z^2RT}{4\pi G} \dfrac{\rho'}{\rho} = \int_0^z z'^2 \rho(z') dz'$$

Now we can differentiate both sides with respect to $$z$$, using the first part of the Fundamental Theorem of Calculus for the right-hand side.

$$\dfrac{2zRT + z^2RT'}{4\pi G} \dfrac{\rho'}{\rho} + \dfrac{z^2RT}{4\pi G} \dfrac{\rho\rho''-(\rho')^2}{\rho^2} = z^2\rho$$

$$(2zRT + z^2RT')\rho\rho' + z^2RT(\rho\rho''-(\rho')^2) = 4\pi Gz^2\rho^3$$

I definitely don't know how to solve that analytically, but with a given temperature profile you should be able to solve numerically using the form $$\rho'' = \dfrac{4\pi G}{RT}\rho^2 + \dfrac{(\rho')^2}{\rho}-\dfrac{2T+zT'}{zT}\rho'$$, from here I'd probably split it into a system and do Runge-Kutta, but the exact best way to do that is beyond my personal knowledge.

• Trying to find how to solve that numerically in Python haha – ordptt Apr 12 at 15:57
• I'm not an expert on numerical methods but I could bang something simple out if it's giving you trouble. It wouldn't be perfect but it should at least be usable. – Stephen Donovan Apr 12 at 16:08
• I know that I most certainly need a initial condition $\rho(0)$, but do I also need a initial condition for $\rho'(0)$? Or is there a way to calculate it? – ordptt Apr 13 at 18:43
• Do you have that $\rho’$ must be continuous? If so I think there might be a trick with taking the limit of your original equation as $z$ goes to zero. That said I’m not really sure if that works considering your original equation doesn’t work at $z = 0$. – Stephen Donovan Apr 14 at 2:44

Here's a different rewrite of your integro-differential equation. I'm swallowing all of the constants up into one, called $$C$$.

$$\frac{1}{\rho(z)}\frac{d\rho}{dz} = \frac{C}{z^2}\int_0^z \zeta^2 \rho(\zeta)\,d\zeta$$

The left hand side can be seen as $$\displaystyle\frac{d}{dz}\log(\rho(z))$$. Differentiating both sides with respect to $$z$$ gives

$$\frac{d^2}{dz^2}\log(\rho(z)) = -\frac{2C}{z^3}\int_0^z \zeta^2 \rho(\zeta)\,d\zeta + \frac{C}{z^2} z^2 \rho(z).$$

The middle term can be recognized as $$\displaystyle-\frac{2}{z}\frac{d}{dz}\log(\rho)$$, giving

$$\frac{d^2}{dz^2} \log(\rho(z)) + \frac{2}{z} \frac{d}{dz}\log(\rho(z)) = C\rho(z).$$

This can be recast as

$$\frac{1}{z^2}\frac{d}{dz} \bigg(z^2\frac{d}{dz} \log(\rho(z))\bigg) = C\rho(z).$$

The operator on the left hand side is closely related to the radial Laplacian, which makes me think this might have been where you started? At any rate, this is horribly non-linear, so a numerical method is your best bet for sure.

Not a full answer, but here's what I have:

I isolated the integral, and then performed integration by parts twice. I ended up defining a function $$Q(z)$$, such that $$Q'''(z) = \rho(z)$$. This gives a nonlinear ODE:

$$\frac{z^2RT}{4\pi G}\frac{Q''''}{Q'''} = z^2Q'' - 2zQ' + 2Q(z)$$ where I used the arbitrariness of constants to define $$Q(0) = 0$$. I can see that $$\frac{d}{dz}\ln Q''' = \frac{Q''''}{Q'''}$$, so if we went through the slog of integrating one more time, repeating the integration by parts trick, we would finally have an (I think) 5th order linear ODE with non constant coefficients. I don't know how useful that is to you with the techniques you know, but my impression is that linear ODEs are far more tractable.