solving the Thomas-Fermi-Dirac equation for planetary interior structure modeling

Question

I am currently working on implementing high pressure equations of state to model interiors of planets and exoplanets. To this end, what I need to do is to numerically solve the Thomas-Fermi-Dirac equation for a given material. The problem I have is purely mathematical, i.e. I am having trouble solving the differential equation. For the simpler case of the Thomas-Fermi equation I successfully proceeded in analogy to S. Esposito 2001 (link 1). After lengthy research I found in F. Ren 2014 (link 2) that they transformed the TFD equation and applied a 4th order Runge-Kutta scheme to solve it. I am, however, struggling to practically apply this idea my self. The equation I seek to solve is:

$$ {d^2 \phi}/{d x^2} = x(\epsilon + (\phi/x)^{1/2} )^3$$

With the initial conditions:

$$ d^2\phi/dx^2|_{x_0} = \phi(x_0)/x_0, \ \ \ \phi(0) = 1 $$

The first problem I have is, that the rhs depends on x and $\phi$ itself rather than on x and the first derivative, which would probably be easier. The second thing is that the initial conditions are given for the second derivative and the function itself at two different points. Of course, transforming the equation into a first order one would considerably simplify things. The simplest substitution $ f (\phi,x) = d \phi /dx$ would yield an integral term on the rhs so another approach seems to be required. Since I have no experience with 2nd order nonlinear differential equations I would appreciate it if some one could guide me through this one.

1: https://core.ac.uk/download/pdf/25325537.pdf

2: https://pdfs.semanticscholar.org/f0b8/0210ae2c2da807cd874c5e8c4c37091a4a46.pdf

Answer 1

Plug the differential equation into the first boundary condition, and set $\phi(x_0)^{1/2}=\psi$:

$$x_0^2(\epsilon+x_0^{-1/2} \psi)^3-\psi^2=0.$$

This is a cubic equation for $\psi$, in more detail

$$x_0^{1/2} \psi^3 + (3\epsilon x_0 - 1) \psi^2 + (3\epsilon^2 x_0^{3/2}) \psi + \epsilon^3=0.$$

If you knew which solution for $\psi$ you would have, then this would be a standard second order boundary value problem, with the value of $\phi$ given at both endpoints. Such a thing is usually solved by a collocation method; for example you can look at Matlab's bvp4c or bvp5c. It can also be solved by a shooting method, which can be achieved by using an IVP solver in combination with a scalar root finder.

If you don't know which $\psi$ you should have, then the solution to the problem is non-unique anyway, so you can try computing the solution to the problem for each of the admissible values of $\psi$ using one of the above methods. Hopefully the physics selects a particular one of them.

To aid in doing that selection, you can do some perturbation analysis.

Consider $\psi=\epsilon^\alpha \psi_0$ (where $\alpha$ is not necessarily positive), then you have

$$x_0^{1/2} \epsilon^{3\alpha} \psi_0^3 + (3\epsilon x_0-1) \epsilon^{2\alpha} \psi_0^2 + 3\epsilon^{2+\alpha} x_0^{3/2} \psi_0 + \epsilon^3=0.$$

So now there are a bunch of powers: $3\alpha,2\alpha,2+\alpha,2\alpha+1$ and $3$. There are now ${5 \choose 2}=10$ balances to consider:

1. $3\alpha = 2 \alpha \Rightarrow \alpha=0$, so you neglect $\epsilon^2,\epsilon^1$ and $\epsilon^3$ and cancel $\epsilon^0$. This is good, because the term being cancelled is the biggest.
2. $3\alpha=2+\alpha \Rightarrow \alpha=1$, so you neglect $\epsilon^2$ and cancel $\epsilon^3$. This is bad, because a term being ignored is bigger than the terms being canceled.
3. $3\alpha=2\alpha+1 \Rightarrow \alpha=1$, this is case 2.
4. $3\alpha=3 \Rightarrow \alpha=1$, this is case 2.
5. $2\alpha=2+\alpha \Rightarrow \alpha=2$, so you neglect $\epsilon^6,\epsilon^5$ and $\epsilon^3$ and cancel $\epsilon^4$. This is bad again.
6. $2\alpha=2\alpha+1$ is impossible.
7. $2\alpha=3 \Rightarrow \alpha=3/2$, so you neglect $\epsilon^{9/2},\epsilon^{7/2}$ and $\epsilon^4$ and cancel $\epsilon^3$. This is good.
8. $2+\alpha=2\alpha+1 \Rightarrow \alpha=1$, this is case 2.
9. $2+\alpha=3 \Rightarrow \alpha=1$, this is case 2.
10. $2\alpha+1=3 \Rightarrow \alpha=1$, this is case 2.

So the only self-consistent scalings are $O(1)$ and $O(\epsilon^{3/2})$. The leading order equations in each case are

$$x_0^{1/2} \psi_0^3 - \psi_0^2 = 0 \\ -\psi_0^2 + 1 = 0.$$

The consistent solution to the first equation is $\psi_0=x_0^{-1/2}$ which corresponds to $\psi=x_0^{-1/2}$. There are two solutions $\psi_0=\pm 1$ to the second equation which correspond to $\psi=\pm \epsilon^{3/2}$.

So you have solutions $x_0^{-1/2}$ (with forward and backward error both $O(\epsilon)$) and $\pm \epsilon^{3/2}$ (with forward error $O(\epsilon^{7/2})$ from which you can read off that the backward error is $O(\epsilon^2)$). Presumably you discard $-\epsilon^{3/2}$ as nonphysical, and then you're left to wonder which of the other two is physical. Having these in hand, I wouldn't continue to do perturbation analysis even if I wanted more accuracy, I would just use these as seeds for an iterative root finder.

From comments it appears the physical solution is $\psi=\epsilon^{3/2}+\dots$, so having solved this polynomial equation to the desired level of accuracy you can now pass off the problem to your favorite ODE BVP solver.