Editing Appendix/Ramblings/51BiPolytropeStability/BetterInterfacePt2 (section)

====Improve First Approximation====

After slogging through the preceding "Second thru Fifth" approximations, I have come to appreciate that the approach used way back in the "[[Appendix/Ramblings/51BiPolytropeStability/BetterInterface#First_Approximation|First Approximation]]" was a good one, but that the attempt to introduce an implicit dependence was misguided.  I now think I have discovered the preferable implicit treatment.  Let's repeat, while improving this approach.

=====2<sup>nd</sup>-Order Explicit Approach=====
As was done in our earlier [[Appendix/Ramblings/51BiPolytropeStability/BetterInterface#First_Approximation|First Approximation]], let's set up a grid associated with a uniformly spaced spherical radius, where the subscript <math>J</math> denotes the grid zone at which all terms in the finite-difference representation of the governing relations will be evaluated.  More specifically,

<table border="0" align="center" cellpadding="5">

<tr>
  <td align="right"><math>\tilde{r}_{J-1}</math></td>
  <td align="center"><math>=</math></td>
  <td align="left">
<math>
\tilde{r}_J - \Delta\tilde{r}
</math>
  </td>
<td align="center">&nbsp; &nbsp; &nbsp; and &nbsp; &nbsp; &nbsp;</td>
  <td align="right"><math>\tilde{r}_{J+1}</math></td>
  <td align="center"><math>=</math></td>
  <td align="left">
<math>
\tilde{r}_J + \Delta\tilde{r} \, ;
</math>
  </td>
</tr>
</table>
also,
<table border="0" align="center" cellpadding="5">

<tr>
  <td align="right"><math>\biggl(\frac{dx}{d\tilde{r}}\biggr)_{J}</math></td>
  <td align="center"><math>\approx</math></td>
  <td align="left">
<math>
\frac{(x_{J+1} - x_{J-1})}{2\Delta\tilde{r}}
</math>
  </td>
<td align="center">&nbsp; &nbsp; &nbsp; and &nbsp; &nbsp; &nbsp;</td>
  <td align="right"><math>\biggl(\frac{dp}{d\tilde{r}}\biggr)_{J}</math></td>
  <td align="center"><math>\approx</math></td>
  <td align="left">
<math>
\frac{(p_{J+1} - p_{J-1})}{2\Delta\tilde{r}} \, .
</math>
  </td>
</tr>
</table>
And at each grid location, the governing relations establish the local evaluation of the derivatives, that is,

<table border="0" align="center" cellpadding="5">

<tr>
  <td align="right"><math>\biggl(\frac{dx}{d \tilde{r}}\biggr)_J </math></td>
  <td align="center"><math>=</math></td>
  <td align="left">
<math>
- \frac{1}{\tilde{r}_J}\biggl[ 
3x + \frac{p}{\gamma_g}\biggr]_J \, ,
</math>
  </td>
<td align="center">&nbsp; &nbsp; &nbsp; and &nbsp; &nbsp; &nbsp;</td>
  <td align="right"><math>\biggl(\frac{dp}{d \tilde{r}}\biggr)_J </math></td>
  <td align="center"><math>=</math></td>
  <td align="left">
<math>
\frac{\tilde{\rho}_J}{\tilde{P}_J}\biggl[ 
(4x + p)\frac{\tilde{M}_r}{\tilde{r}^2}
+
\tau_c^2 \omega^2 \tilde{r} x \biggr]_J \, .
</math>
  </td>
</tr>
</table>
<span id="1stapprox">So, integrating</span> step-by-step from the center of the configuration, outward, once all the variable values are known at grid locations <math>J</math> and <math>(J-1)</math>, the values of <math>x</math> and <math>p</math> at <math>(J+1)</math> are given by the expressions,
<table border="0" align="center" cellpadding="5">

<tr>
  <td align="right"><math>x_{J+1}</math></td>
  <td align="center"><math>=</math></td>
  <td align="left">
<math>
x_{J-1} - 2\Delta\tilde{r} \biggl\{
\frac{1}{\tilde{r}}\biggl[ 
3x + \frac{p}{\gamma_g}\biggr]
\biggr\}_J \, ,
</math>
  </td>
<td align="center">&nbsp; &nbsp; &nbsp; and &nbsp; &nbsp; &nbsp;</td>
  <td align="right"><math>p_{J+1}</math></td>
  <td align="center"><math>=</math></td>
  <td align="left">
<math>
p_{J-1} + 2\Delta\tilde{r}
\biggl\{
\frac{\tilde{\rho}}{\tilde{P}} \cdot \frac{\tilde{M}_r}{\tilde{r}^2}
\biggl[ (4x + p)
+
\sigma_c^2 \biggl(\frac{2\pi}{3} \cdot \frac{\tilde{\rho}_c \tilde{r}^3 }{\tilde{M}_r}\biggr) x \biggr] 
\biggr\}_J\, .
</math>
  </td>
</tr>
</table>

Then we will obtain the "<math>x_J</math>" and "<math>p_J</math>" values via the ''average'' expressions,
<table border="0" align="center" cellpadding="5">

<tr>
  <td align="right"><math>x_{J}</math></td>
  <td align="center"><math>=</math></td>
  <td align="left">
<math>
\frac{1}{2}(x_{J-1} + x_{J+1})
\, ,
</math>
  </td>
<td align="center">&nbsp; &nbsp; &nbsp; and &nbsp; &nbsp; &nbsp;</td>
  <td align="right"><math>p_{J}</math></td>
  <td align="center"><math>=</math></td>
  <td align="left">
<math>
\frac{1}{2}(p_{J-1} + p_{J+1})
\, .
</math>
  </td>
</tr>
</table>

=====Convert to Implicit Approach=====

Consider implementing an ''implicit'' finite-difference analysis that improves on our earlier [[Appendix/Ramblings/51BiPolytropeStability/BetterInterface#First_Approximation|First Approximation]].  The general form of the source term expressions is,
<table border="0" align="center" cellpadding="5">

<tr>
  <td align="right"><math>x_{J+1}</math></td>
  <td align="center"><math>=</math></td>
  <td align="left">
<math>
x_{J-1} + 2\Delta\tilde{r} \biggl\{
\mathfrak{A}x_J + \mathfrak{B}p_J
\biggr\} 
</math>
  </td>
</tr>
</table>
where,
<table border="0" align="center" cellpadding="5">

<tr>
  <td align="right"><math>\mathfrak{A}</math></td>
  <td align="center"><math>\equiv</math></td>
  <td align="left">
<math>
- \biggl\{ \frac{3}{\tilde{r}} \biggr\}_J \, ,
</math>
  </td>
<td align="center">&nbsp; &nbsp; &nbsp; and &nbsp; &nbsp; &nbsp;</td>
  <td align="right"><math>\mathfrak{B}</math></td>
  <td align="center"><math>\equiv</math></td>
  <td align="left">
<math>
- \biggl\{ \frac{1}{\gamma_g \tilde{r}} \biggr\}_J \, ;
</math>
  </td>
</tr>
</table>
and,

<table border="0" align="center" cellpadding="5">

<tr>
  <td align="right"><math>p_{J+1}</math></td>
  <td align="center"><math>=</math></td>
  <td align="left">
<math>
p_{J-1} + 2\Delta\tilde{r} \biggl\{
\mathfrak{C}x_J + \mathfrak{D}p_J
\biggr\} 
</math>
  </td>
</tr>
</table>

where,
<table border="0" align="center" cellpadding="5">

<tr>
  <td align="right"><math>\mathfrak{C}</math></td>
  <td align="center"><math>\equiv</math></td>
  <td align="left">
<math>
\biggl\{
\mathfrak{D}
\biggl[ 4
+
\sigma_c^2 \biggl(\frac{2\pi}{3} \cdot \frac{\tilde{\rho}_c \tilde{r}^3 }{\tilde{M}_r}\biggr)  \biggr] 
\biggr\}_J\, ,
</math>
  </td>
<td align="center">&nbsp; &nbsp; &nbsp; and &nbsp; &nbsp; &nbsp;</td>
  <td align="right"><math>\mathfrak{D}</math></td>
  <td align="center"><math>\equiv</math></td>
  <td align="left">
<math>
\biggl\{
\frac{\tilde{\rho}}{\tilde{P}} \cdot \frac{\tilde{M}_r}{\tilde{r}^2}
\biggr\}_J\, .
</math>
  </td>
</tr>
</table>

Now, wherever a "<math>J+1</math>" index appears in the source term, replace it with the ''average expressions''; specifically, <math>x_{J+1} \rightarrow (2 x_J - x_{J-1})</math> and <math>p_{J+1} \rightarrow (2 p_J - p_{J-1})</math>.  For the fractional radial displacement, we have,
<table border="0" align="center" cellpadding="5">

<tr>
  <td align="right"><math>x_J </math></td>
  <td align="center"><math>=</math></td>
  <td align="left">
<math>
x_{J-1} + \Delta\tilde{r} \biggl\{
\mathfrak{A}x_J + \mathfrak{B}p_J
\biggr\} \, ;
</math>
  </td>
</tr>
</table>
and for the fractional pressure displacement,

<table border="0" align="center" cellpadding="5">

<tr>
  <td align="right"><math>p_{J} </math></td>
  <td align="center"><math>=</math></td>
  <td align="left">
<math>
p_{J-1} + \Delta\tilde{r} \biggl\{
\mathfrak{C}x_J + \mathfrak{D}p_J
\biggr\} \, .
</math>
  </td>
</tr>
</table>
Solving for <math>p_J</math> in this second expression, we obtain,

<table border="0" align="center" cellpadding="5">

<tr>
  <td align="right"><math>p_{J}\biggl[ 1 - \mathfrak{D}\Delta\tilde{r}\biggr] </math></td>
  <td align="center"><math>=</math></td>
  <td align="left">
<math>
p_{J-1} + (\mathfrak{C}\Delta\tilde{r}) x_J 
</math>
  </td>
</tr>
</table>
in which case the first expression gives,
<table border="0" align="center" cellpadding="5">

<tr>
  <td align="right"><math>x_J </math></td>
  <td align="center"><math>=</math></td>
  <td align="left">
<math>
x_{J-1} + (\mathfrak{A} \Delta\tilde{r})x_J 
+ (\mathfrak{B} \Delta\tilde{r}) \biggl[ p_{J-1} + (\mathfrak{C}\Delta\tilde{r}) x_J\biggr]\biggl[ 1 - \mathfrak{D}\Delta\tilde{r}\biggr]^{-1}
</math>
  </td>
</tr>

<tr>
  <td align="right"><math>\Rightarrow ~~~ x_J \biggl[ 1 - \mathfrak{D}\Delta\tilde{r}\biggr]</math></td>
  <td align="center"><math>=</math></td>
  <td align="left">
<math>
\biggl\{
x_{J-1} + (\mathfrak{A} \Delta\tilde{r})x_J 
\biggr\}\biggl[ 1 - \mathfrak{D}\Delta\tilde{r}\biggr]
+ (\mathfrak{B} \Delta\tilde{r}) \biggl[ p_{J-1} + (\mathfrak{C}\Delta\tilde{r}) x_J\biggr]
</math>
  </td>
</tr>

<tr>
  <td align="right">&nbsp;</td>
  <td align="center"><math>=</math></td>
  <td align="left">
<math>
x_{J-1} \biggl[ 1 - \mathfrak{D}\Delta\tilde{r}\biggr]
+ (\mathfrak{B} \Delta\tilde{r}) p_{J-1} 
+
(\mathfrak{A} \Delta\tilde{r})x_J \biggl[ 1 - \mathfrak{D}\Delta\tilde{r}\biggr]
+ (\mathfrak{B} \Delta\tilde{r}) (\mathfrak{C}\Delta\tilde{r}) x_J
</math>
  </td>
</tr>

<tr>
  <td align="right">&nbsp;</td>
  <td align="center"><math>=</math></td>
  <td align="left">
<math>
x_{J-1} ( 1 - \mathfrak{D}\Delta\tilde{r})
+ (\mathfrak{B} \Delta\tilde{r}) p_{J-1} 
+\biggl[
(\mathfrak{A} \Delta\tilde{r}) ( 1 - \mathfrak{D}\Delta\tilde{r})
+ (\mathfrak{B} \Delta\tilde{r}) (\mathfrak{C}\Delta\tilde{r}) 
\biggr] x_J
</math>
  </td>
</tr>

<tr>
  <td align="right"><math>
\Rightarrow ~~~ 
x_J \biggl\{ ( 1 - \mathfrak{D}\Delta\tilde{r} )- \biggl[
(\mathfrak{A} \Delta\tilde{r}) ( 1 - \mathfrak{D}\Delta\tilde{r})
+ (\mathfrak{B} \Delta\tilde{r}) (\mathfrak{C}\Delta\tilde{r}) 
\biggr] \biggr\}
</math></td>
  <td align="center"><math>=</math></td>
  <td align="left">
<math>
x_{J-1} ( 1 - \mathfrak{D}\Delta\tilde{r})
+ (\mathfrak{B} \Delta\tilde{r}) p_{J-1} 
</math>
  </td>
</tr>

<tr>
  <td align="right"><math>
\Rightarrow ~~~ 
x_J \biggl\{ 1 - \biggl[
(\mathfrak{A} \Delta\tilde{r}) 
+ \frac{(\mathfrak{B} \Delta\tilde{r}) (\mathfrak{C}\Delta\tilde{r})}{ ( 1 - \mathfrak{D}\Delta\tilde{r} ) } 
\biggr] \biggr\}
</math></td>
  <td align="center"><math>=</math></td>
  <td align="left">
<math>
x_{J-1} 
+ \biggl[\frac{(\mathfrak{B} \Delta\tilde{r})}{( 1 - \mathfrak{D}\Delta\tilde{r} )}\biggr] p_{J-1} 
</math>
  </td>
</tr>

<tr>
  <td align="right"><math>
\Rightarrow ~~~ 
x_J 
</math></td>
  <td align="center"><math>=</math></td>
  <td align="left">
<math>
\biggl\{ x_{J-1} 
+ \biggl[\frac{(\mathfrak{B} \Delta\tilde{r})}{( 1 - \mathfrak{D}\Delta\tilde{r} )}\biggr] p_{J-1} \biggr\}
\biggl\{ 1 - \biggl[
(\mathfrak{A} \Delta\tilde{r}) 
+ \frac{(\mathfrak{B} \Delta\tilde{r}) (\mathfrak{C}\Delta\tilde{r})}{ ( 1 - \mathfrak{D}\Delta\tilde{r} ) } 
\biggr] \biggr\}^{-1}
\, .
</math>
  </td>
</tr>
</table>
Then,

<table border="0" align="center" cellpadding="5">

<tr>
  <td align="right"><math>p_{J}( 1 - \mathfrak{D}\Delta\tilde{r} ) </math></td>
  <td align="center"><math>=</math></td>
  <td align="left">
<math>
p_{J-1} + (\mathfrak{C}\Delta\tilde{r}) x_J 
</math>
  </td>
</tr>

<tr>
  <td align="right">&nbsp;</td>
  <td align="center"><math>=</math></td>
  <td align="left">
<math>
p_{J-1} + (\mathfrak{C}\Delta\tilde{r}) 
\biggl\{ x_{J-1} 
+ \biggl[\frac{(\mathfrak{B} \Delta\tilde{r})}{( 1 - \mathfrak{D}\Delta\tilde{r} )}\biggr] p_{J-1} \biggr\}
\biggl\{ 1 - \biggl[
(\mathfrak{A} \Delta\tilde{r}) 
+ \frac{(\mathfrak{B} \Delta\tilde{r}) (\mathfrak{C}\Delta\tilde{r})}{ ( 1 - \mathfrak{D}\Delta\tilde{r} ) } 
\biggr] \biggr\}^{-1}

</math>
  </td>
</tr>
</table>

<table border=1 align="center" cellpadding="10" width="80%"><tr><td bgcolor="lightblue" align="left">
This is test of  our "implicit" scheme for the <math>(n_c, n_e) = (5, 1)</math> bipolytrope with <math>\mu_e/\mu_c = 0.31</math> and (Model A) <math>\xi_i = 9.12744</math>; here, we also assume <math>\sigma_c^2 = 0.000109</math> and <math>J = i+2</math>.  Here are the quantities that we assume are <b>known</b> &hellip;
<table border="1" align="center" cellpadding="5">

  <tr>
<td align="center" bgcolor="white"><math>\tilde{r}</math></td>
<td align="center" bgcolor="white"><math>\tilde\rho</math></td>
<td align="center" bgcolor="white"><math>\tilde{P}</math></td>
<td align="center" bgcolor="white"><math>\tilde{M}_r</math></td>
<td align="center" bgcolor="white"><math>\tilde{\rho}_c</math></td>

<td align="center" bgcolor="white"><math>\mathfrak{A}</math></td>
<td align="center" bgcolor="white"><math>\mathfrak{B}</math></td>
<td align="center" bgcolor="white"><math>\frac{\mathfrak{C}}{\mathfrak{D}}</math></td>
<td align="center" bgcolor="white"><math>\mathfrak{D}</math></td>
  </tr>

  <tr>
<td align="center" bgcolor="white">0.0193368</td>
<td align="center" bgcolor="white">192.21728</td>
<td align="center" bgcolor="white">1913.1421</td>
<td align="center" bgcolor="white">0.3403116</td>
<td align="center" bgcolor="white">3359266.406</td>

<td align="center" bgcolor="white">-155.14459</td>
<td align="center" bgcolor="white">-25.85743</td>
<td align="center" bgcolor="white">4.0162932</td>
<td align="center" bgcolor="white">91.443479</td>
  </tr>
</table>
<table border="1" align="center" cellpadding="5">

  <tr>
<td align="center" bgcolor="white"><math>\Delta\tilde{r}</math></td>
<td align="center" bgcolor="white"><math>x_{J-1}</math></td>
<td align="center" bgcolor="white"><math>p_{J-1}</math></td>
<td align="center" bgcolor="white"><font size="-1">determined</font><br /><math>x_J</math></td>
<td align="center" bgcolor="white"><font size="-1">determined</font><br /><math>p_J</math></td>
  </tr>

  <tr>
<td align="center" bgcolor="white">0.001936393</td>
<td align="center" bgcolor="white">-4.755073</td>
<td align="center" bgcolor="white">32.25497</td>
<td align="center" bgcolor="white">-4.999355</td>
<td align="center" bgcolor="white">34.874915</td>
  </tr>
</table>
</td></tr>
<tr><td bgcolor="white" align="center">

<table border="0" align="center" cellpadding="5">

<tr>
  <td align="right"><math>
x_J 
</math></td>
  <td align="center"><math>=</math></td>
  <td align="left">
<math>
\biggl\{ x_{J-1} 
+ \biggl[\frac{(\mathfrak{B} \Delta\tilde{r})}{( 1 - \mathfrak{D}\Delta\tilde{r} )}\biggr] p_{J-1} \biggr\}
\biggl\{ 1 - \biggl[
(\mathfrak{A} \Delta\tilde{r}) 
+ \frac{(\mathfrak{B} \Delta\tilde{r}) (\mathfrak{C}\Delta\tilde{r})}{ ( 1 - \mathfrak{D}\Delta\tilde{r} ) } 
\biggr] \biggr\}^{-1}
</math>
  </td>
</tr>

<tr>
  <td align="right">&nbsp;</td>
  <td align="center"><math>=</math></td>
  <td align="left">
<math>
\biggl\{ x_{J-1} 
+ \biggl[-0.06084379\biggr] p_{J-1} \biggr\}
\biggl\{ 1 - \biggl[-0.3436910 \biggr] \biggr\}^{-1}
=
-4.999354 \, ;
</math>
  </td>
</tr>

<tr>
  <td align="right"><math>p_{J}( 1 - \mathfrak{D}\Delta\tilde{r} ) </math></td>
  <td align="center"><math>=</math></td>
  <td align="left">
<math>
\biggl[ p_{J-1} + (\mathfrak{C}\Delta\tilde{r}) x_J \biggr]( 1 - \mathfrak{D}\Delta\tilde{r} )^{-1}
=
34.87491 \, .
</math>
  </td>
</tr>
</table>
<tr><td bgcolor="lightblue" align="left">
<b>Best values:</b><br />Nodes 0: &nbsp;n/a
<br />Nodes 1: &nbsp;<math>\sigma_c^2 = 8.958784\times 10^{-5}</math>
<br />Nodes 2: &nbsp;<math>\sigma_c^2 = 3.021\times 10^{-4}</math>
<br />Nodes 3: &nbsp;<math>\sigma_c^2 = 6.09\times 10^{-4}</math>
</td></tr>
</table>
</td></tr></table>