Editing SSC/Structure/Polytropes/Numerical (section)

===Techniques===
====HSCF Technique====
{| class="PGEclass" style="float:left; margin-right: 20px; border-style: solid; border-width: 3px border-color: black"
|- 
! style="height: 125px; width: 125px; background-color:white;" |
<font size="-1">[[H_BookTiledMenu#Equilibrium_Structures|<b>via<br />Self-Consistent<br />Field (SCF)<br />Technique</b>]]</font>
|}
On the first spreadsheet within the workbook, we establish the following columns of number:
<ul>
  <li>'''Column A:''' &nbsp; Labeled <math>~r_i/R</math> (for ''i'' between 1 and N), that represents a discrete radial grid of spacing, <math>~~\Delta = (N-1)^{-1}</math>; each row gives the radial coordinate location of the ''i''<sup>th</sup> zone, starting from <math>~r_1/R = 0</math> and ending at <math>~r_N/R = 1</math>.</li>
  <li>'''Column B:''' &nbsp; Labeled <math>~\mathrm{rhf}_i</math> (for ''i'' between 1 and N-1); each row gives the radial coordinate of the mid-point of a grid zone.</li>
<div align="center">
<math>~\mathrm{rhf}_i \equiv \frac{1}{2}\biggr[\frac{r_i}{R} + \frac{r_{i+1}}{R}  \biggr] \, .</math>
</div>
  <li>'''Column C:''' &nbsp; Labeled <math>~\rho_i</math> (for ''i'' between 1 and N-1); each row provides an initial ''guess'' for the mass-density of the grid zone.  Usually it is sufficient to guess, <math>~\rho_i = 1</math> throughout.  For an <math>~n=0</math> polytrope, this proves also to be the correct ''final'' density profile.</li>
  <li>'''Column D:''' &nbsp; Labeled <math>~M_i</math> (for ''i'' between 1 and N); the ''i''<sup>th</sup> row gives the integrated mass enclosed interior to the radial grid coordinate, <math>~r_i/R</math>.  Specifically, <math>~M_1 = 0</math>, and thereafter, beginning with zone, <math>~i = 2</math>,</li>
<div align="center">
<math>~M_i = M_{i-1} + \frac{4\pi \rho_{i-1}}{3}\biggl[ \biggl(\frac{r_{i}}{R} \biggr)^3 - \biggl(\frac{r_{i-1}}{R} \biggr)^3  \biggr] \, .</math>
</div>
  <li>Note that, <math>~M_\mathrm{tot} = M_N \, .</math></li>
  <li>'''Column E:''' &nbsp; Labeled <math>~g_i</math> (for ''i'' between 2 and N); each row tabulates the inwardly directed gravitational acceleration that is felt at the outer edge of each grid zone.  Specifically,</li>
<div align="center">
<math>g_i = \frac{GM_i}{ (r_i/R)^2} \, .</math>
</div>
  <li>'''Column F:''' &nbsp; Labeled <math>~\Phi_i</math> (for ''i'' between 1 and N); each row gives the value of the gravitational potential at the mid-point of a grid zone.  Here, we start by specifying the ''value'' of the potential just (specifically, half a radial grid-zone) outside the surface of the configuration, where it should be, <math>~\Phi_N = -GM_\mathrm{tot}/(1+\Delta/2)</math>.  Then, working from the surface, inward and, given that, <math>~g = d\Phi/dr</math>, we use the corresponding finite-difference representation of the radial derivative and set,</li>
<div align="center">
<table border="0" cellpadding="5" align="center">

<tr>
  <td align="right">
<math>~\frac{\Phi_i - \Phi_{i-1}}{\Delta} </math>
  </td>
  <td align="center">
<math>~=</math>
  </td>
  <td align="left">
<math>~g_i </math>
  </td>
</tr>

<tr>
  <td align="right">
<math>~\Rightarrow ~~~ \Phi_{i-1} </math>
  </td>
  <td align="center">
<math>~=</math>
  </td>
  <td align="left">
<math>~ \Phi_i - g_i \Delta \, .</math>
  </td>
</tr>
</table>
</div>
<li>Note that the value of the gravitational potential ''at'' the surface is not <math>~\Phi_N</math> but, rather, must be <math>~\Phi_\mathrm{surf} = -GM_\mathrm{tot}/R</math>.</li>
<li>Furthermore, note that a lop-sided Taylor-series expansion about the center of the configuration provides the following good approximation to the gravitational potential ''at'' the center: &nbsp; <math>~\Phi_c \approx (9\Phi_1 - \Phi_2)/8</math>.</li>
<li>Note as well that all of these numerically determined values of the gravitational potential can be checked against the [[SSC/Structure/UniformDensity#UniformSpherePotential|known analytic expression]] for the radial profile of the potential in a uniform-density sphere.</li>
<li>'''Column G:''' &nbsp; Labeled <math>~H_i</math> (for ''i'' between 1 and N-1); each row provides the value of the fluid enthalpy at the center of a grid cell.  Adopting the convention that the enthalpy is zero at the surface of the configuration, and given that [[SSCpt2/SolutionStrategies#Technique_3|the enthalpy and the gravitational potential must sum to zero]] throughout the configuration, we have,</li>
<div align="center">
<math>H_i = \Phi_\mathrm{surf} - \Phi_i \, .</math>
</div>
<li>At the center of the configuration, we have, <math>~H_c = \Phi_\mathrm{surf} - \Phi_c</math>. </li>
<li>'''Column H:''' &nbsp; Labeled <math>~H_\mathrm{norm}</math> (for ''i'' between 1 and N-1); each row provides the value of the fluid enthalpy, renormalized to the central value, specifically,</li>
<div align="center">
<math>~[H_\mathrm{norm}]_i = \frac{H_i}{H_c} \, .</math>
</div>
</ul>


The second spreadsheet within the workbook should be initially created by generating a copy of the first spreadsheet.  Then:
<ul>
  <li>'''Column C:''' &nbsp; Labeled <math>~\rho_i</math> (for ''i'' between 1 and N-1); generate a new, improved ''guess'' for the normalized mass-density at each grid zone based on the corresponding value of the normalized enthalpy from the previous spreadsheet/iteration. Specifically, given that the [[SR#Barotropic_Structure|relationship between the density and enthalpy in a polytrope]] of index, <math>~n</math>, is, <math>~\rho \propto H^n</math>, we should set,</li>
<div align="center">
<math>~\biggl\{ \frac{\rho_i}{\rho_c} \biggr\}_\mathrm{sheet2}= \biggr\{[H_\mathrm{norm}]_i^n \biggr\}_\mathrm{sheet1} \, .</math>
</div>
</ul>

====Straight Numerical Integration====
{| class="PGEclass" style="float:left; margin-right: 20px; border-style: solid; border-width: 3px border-color: black"
|- 
! style="height: 125px; width: 125px; background-color:white;" |
<font size="-1">[[H_BookTiledMenu#Equilibrium_Structures|<b>via<br />Direct<br />Numerical<br />Integration</b>]]</font>
|}
The [[#LaneEmdenEquation|above governing relation]] may be rewritten as,
<div align="center">
<table border="0" cellpadding="5" align="center">

<tr>
  <td align="right">
<math>~ 
\xi \theta^{''} + 2 \theta^' 
</math>
  </td>
  <td align="center">
<math>~=</math>
  </td>
  <td align="left">
<math>~- \xi \theta^n \, .</math>
  </td>
</tr>
</table>
</div>

We'll adopt the following finite-difference approximations for the first and second derivatives on a grid of radial spacing, <math>~\Delta_\xi</math>:
<div align="center">
<table border="0" cellpadding="5" align="center">

<tr>
  <td align="right">
<math>~\theta_i'</math>
  </td>
  <td align="center">
<math>~\approx</math>
  </td>
  <td align="left">
<math>~\frac{\theta_+ - \theta_-}{2\Delta_\xi}</math>
  </td>
</tr>
</table>
</div>
and,
<div align="center">
<table border="0" cellpadding="5" align="center">

<tr>
  <td align="right">
<math>~\theta_i''</math>
  </td>
  <td align="center">
<math>~\approx</math>
  </td>
  <td align="left">
<math>~\frac{\theta_+ - 2\theta_i +\theta_-}{\Delta_\xi^2} \, .</math>
  </td>
</tr>
</table>
</div>
Our finite-difference approximation of the governing equation is, then,
<div align="center">
<table border="0" cellpadding="5" align="center">

<tr>
  <td align="right">
<math>~\xi_i \biggl[ \frac{\theta_+ - 2\theta_i +\theta_-}{\Delta_\xi^2} \biggr] + 2\biggl[  \frac{\theta_+ - \theta_-}{2\Delta_\xi}\biggr]
</math>
  </td>
  <td align="center">
<math>~=</math>
  </td>
  <td align="left">
<math>~- \xi_i \theta_i^n </math>
  </td>
</tr>

<tr>
  <td align="right">
<math>~\Rightarrow ~~~ \xi_i [ \theta_+ - 2\theta_i +\theta_-] + \Delta_\xi [  \theta_+ - \theta_- ]
</math>
  </td>
  <td align="center">
<math>~=</math>
  </td>
  <td align="left">
<math>~- \Delta_\xi^2\xi_i \theta_i^n </math>
  </td>
</tr>

<tr>
  <td align="right">
<math>~\Rightarrow ~~~ 
\theta_+
</math>
  </td>
  <td align="center">
<math>~=</math>
  </td>
  <td align="left">
<math>~\frac{2\xi_i \theta_i  + \theta_-(\Delta_\xi - \xi_i)  - \Delta_\xi^2\xi_i \theta_i^n }{\Delta_\xi + \xi_i} \, .</math>
  </td>
</tr>
</table>
</div>

Now, for the first two steps away from the center &#8212; where, <math>~\theta_i = \theta_0 = 1</math> and <math>~\xi_i = \xi_0 = 0</math> &#8212;  we will use the following [[Appendix/Ramblings/PowerSeriesExpressions#PolytropicLaneEmden|power-series expansion]] (see, for example, eq. 62 from &sect;5 in Chapter IV of [<b>[[Appendix/References#C67|<font color="red">C67</font>]]</b>]) to determine the value of <math>~\theta_i</math>:
<div align="center">
<table border="0" cellpadding="5" align="center">

<tr>
  <td align="right">
<math>~\theta_1 
</math>
  </td>
  <td align="center">
<math>~=</math>
  </td>
  <td align="left">
<math>1 - \frac{\Delta_\xi^2}{6} + \frac{n \Delta_\xi^4}{120}  - \frac{n}{378} \biggl( \frac{n}{5} - \frac{1}{8} \biggr) \Delta_\xi^6 \, ,</math> 
  </td>
</tr>
</table>
</div>
and,
<div align="center">
<table border="0" cellpadding="5" align="center">

<tr>
  <td align="right">
<math>~\theta_2 
</math>
  </td>
  <td align="center">
<math>~=</math>
  </td>
  <td align="left">
<math>1 - \frac{(2\Delta_\xi)^2}{6} + \frac{n (2\Delta_\xi)^4}{120}  - \frac{n}{378} \biggl( \frac{n}{5} - \frac{1}{8} \biggr) (2\Delta_\xi)^6 \, .</math> 
  </td>
</tr>
</table>
</div>