Examples of Numerical Solutions

In our numerical approach here, we assume that the external radius $x_{2}=\Upsilon_{0}r_{2}$ is given. In order to complete the calculation we have to determine the interior radius $x_{1}$. This can be done recalling that the dimensionless pressure $p(x)$ is zero for $x=x_{2}$ and $x=x_{1}$. Since according to Equation (53) $p(x)=1/z(x)-1$, this is equivalent to the determination of the values of $x$ for which $z(x)=1$. By the determination of $x_{1}$ we would have solved the problem in the entire matter region. Note that since $x=\Upsilon_{0}r=\sqrt{\kappa\rho_{0}}\;r$ we have obtained a family of solutions parametrized by two parameters, the external radius $r_{2}$ and the parameter $\eta$.

If the discriminant $\Delta\ne 0$ the integral in Equation (84) is expressed in terms of elliptic integrals and the result is not very transparent. It is more convenient to integrate the differential Equation (79) using the fourth-order Runge-Kutta algorithm (RK4) [#!NumericalRecipes!#]. We start by choosing a value of $x=x_{2}$ for which the cubic polynomial is positive and we put $z(x_{2})=1$. This determines the outer radius of the matter shell. We then iterate the differential equation given in Equation (79) in the decreasing $x$ direction until we reach the first point for which the value of $z$ returns to $1$. This point is chosen as $x_{1}$. If a value for $x_{1}$ cannot be found, we conclude that there is no solution to the problem with the given values of $x_{2}$ and $x_{M}$. A good test for the efficiency of the algorithm is to compare the exact analytic result given in Equation (88) with the result from the numerical integration in that same case. These results are shown in Figure 1. On any current $64$-bit desktop computer one can easily reach a high degree of precision with little numerical effort. After iterating the RK4 algorithm from $x_{2}$ to $x_{1}$ the difference between the exact and the numerical results for $z(x)$ stays below $1.03536\times 10^{-29}$ for an iteration step of $\delta x\approx 10^{-7}$.

In the comments that follow $x_{\mu}\equiv\Upsilon_{0}r_\mu$, where $r_\mu$ is the integration constant that results from the solution of the Einstein equations in the inner vacuum region, given in Equation (31). In the matter region the input parameters are $\eta$ and $x_{2}$. The parameter $x_{1}$ is obtained from the iteration of Equation (79). The value of $x_{M}$ that is necessary for plotting the curves is given in Equation (83). The expressions for $\lambda (x)$ and $\nu (x)$ are given in Table 1. Figure 2 shows the plots of the functions $\nu (x)$ and $\lambda (x)$ for $\eta =2.0$ and $x_{2}=5^{1/3}$. The curves were obtained analytically using Equation (88) and the expressions in Table 1, but using the numerically calculated parameters $x_{1}=0.594881$ and $x_{\mu }=0.596494$.

In Figure 3 we plot the dimensionless pressure $p(x)$ as a function of $x$, in a case in which there is no analytic expression in terms of elementary functions and the calculation is performed numerically. The parameters are $x_{1}=1.24050$ and $x_{\mu}=1.03035$. Comparing Figures 1 and 3, that depict the dimensionless pressure $p(x)$ as a function of $x$ for $\eta =2.0$ and $\eta=5.0$, one notes that the two graphs are similar but for larger values of $\eta$ the graph becomes less symmetric.

Figure 4 shows the plots of the functions $\nu (x)$ and $\lambda (x)$, for $\eta=5.0$ and $x_{2}=2.0$. In this case there are no analytical solutions in terms of elementary functions available in the matter region and the values of $\nu (x)$ and $\lambda (x)$ were obtained numerically. In the vacuum regions we used the analytical expressions given in Table 1 with the parameters $x_{1}=1.24050$ and $x_{\mu}=1.03035$.