Editing SSC/Dynamics/IsothermalSimilaritySolution (section)

=Numerical Integration of Coupled ODEs=

Let's develop a finite-difference expression that allows us to straightforwardly integrate the [[#CoupledODEs|above pair of coupled ODEs]], namely,
<div align="center">
<table border="0" cellpadding="5" align="center">

<tr>
  <td align="right">
<math>~\frac{dU}{d\zeta}</math>
  </td>
  <td align="center">
<math>~=</math>
  </td>
  <td align="left">
<math>~
\frac{(\zeta U +1) [\Rho (\zeta U +1) -2 ]}{[ (\zeta U +1)^2 - \zeta^2]} \, ,
</math>
  </td>
</tr>

<tr>
  <td align="right">
<math>~\frac{dP}{d\zeta}</math>
  </td>
  <td align="center">
<math>~=</math>
  </td>
  <td align="left">
<math>~\frac{\zeta \Rho [2-\Rho (\zeta U +1)]}{[ (\zeta U +1)^2 - \zeta^2]} \, .</math>
  </td>
</tr>
</table>
</div>

Assume that the starting values of <math>~U</math> and <math>~\Rho</math> have been provided by, for example, one of the above detailed series expansions.  Let the subscript "1" denote these known values at coordinate-location, <math>~\zeta_1</math>, and let the subscript "2" denote the unknown values of <math>~U</math> and <math>~\Rho</math> at <math>~\zeta_2 = \zeta_1 + \Delta\zeta</math>.  We should be able to construct a 2<sup>nd</sup>-order accurate integration scheme by treating <math>~U</math> and <math>~\Rho</math> as ''average values'' everywhere they occur on the right-hand sides of the pair of ODEs.  That is,
<div align="center">
<table border="0" cellpadding="5" align="center">

<tr>
  <td align="right">
<math>~\bar{U}</math>
  </td>
  <td align="center">
<math>~\rightarrow</math>
  </td>
  <td align="left">
<math>~\frac{1}{2}\biggl(U_2 + U_1 \biggr) \, ,</math>
  </td>
  <td align="center">&nbsp; &nbsp; &nbsp; and &nbsp; &nbsp; &nbsp; 
  <td align="right">
<math>~\bar{\Rho}</math>
  </td>
  <td align="center">
<math>~\rightarrow</math>
  </td>
  <td align="left">
<math>~\frac{1}{2}\biggl(\Rho_2 + \Rho_1 \biggr) \, .</math>
  </td>
</tr>
</table>
</div>

==First ODE==
We'll begin by using the first ODE to provide one expression for <math>~\Rho_2</math> in terms of <math>~U</math> and <math>~\Rho_1</math>.  We have,
<div align="center">
<table border="0" cellpadding="5" align="center">

<tr>
  <td align="right">
<math>~
\bar{\Rho} (\zeta \bar{U} +1) 
</math>
  </td>
  <td align="center">
<math>~=</math>
  </td>
  <td align="left">
<math>~
2+\biggl[ \frac{ (\zeta \bar{U} +1)^2 - \zeta^2  }{ (\zeta \bar{U} +1)  } \biggr]\frac{\Delta U}{\Delta\zeta}
</math>
  </td>
</tr>

<tr>
  <td align="right">
<math>~\Rightarrow~~~
\Rho_2  
</math>
  </td>
  <td align="center">
<math>~=</math>
  </td>
  <td align="left">
<math>~
\frac{2}{(\zeta \bar{U} +1)} \biggl\{ 2+\biggl[ \frac{ (\zeta \bar{U} +1)^2 - \zeta^2  }{ (\zeta \bar{U} +1)  } \biggr]\frac{\Delta U}{\Delta\zeta} \biggr\} - P_1 \, .
</math>
  </td>
</tr>

<tr>
  <td align="right">
&nbsp;
  </td>
  <td align="center">
<math>~=</math>
  </td>
  <td align="left">
<math>~
\frac{2}{\ell} \biggl\{ 2+\biggl[ \frac{ \ell^2 - \zeta^2  }{ \ell  } \biggr]\frac{\Delta U}{\Delta\zeta} \biggr\} - P_1 \, ,
</math>
  </td>
</tr>
</table>
</div>
where,
<div align="center">
<math>~\ell \equiv (\zeta \bar{U} + 1)</math> &nbsp;  &nbsp; &nbsp; <math>~\Rightarrow</math> &nbsp; &nbsp; &nbsp; <math>~\bar{U} =\biggl[\frac{\ell - 1}{\zeta}\biggr] \, .</math>
</div>

==Second ODE==
The second ODE can presumably provide a second, independent expression for <math>~\Rho_2</math> in terms of <math>~U</math> and <math>~\Rho_1</math>.  In this case we have, 
<div align="center">
<table border="0" cellpadding="5" align="center">

<tr>
  <td align="right">
<math>~\frac{(P_2 - P_1)}{\Delta\zeta}</math>
  </td>
  <td align="center">
<math>~=</math>
  </td>
  <td align="left">
<math>~\frac{\zeta \bar{\Rho} [2-\bar{\Rho} (\zeta \bar{U} +1)]}{[ (\zeta \bar{U} +1)^2 - \zeta^2]} </math>
  </td>
</tr>

<tr>
  <td align="right">
&nbsp;
  </td>
  <td align="center">
<math>~=</math>
  </td>
  <td align="left">
<math>~\frac{\zeta \tfrac{1}{2}(\Rho_1 + \Rho_2) [2-\tfrac{1}{2}(\Rho_1 + \Rho_2) (\zeta \bar{U} +1)]}{[ (\zeta \bar{U} +1)^2 - \zeta^2]} </math>
  </td>
</tr>

<tr>
  <td align="right">
<math>~\Rightarrow ~~~ \frac{ 4[ (\zeta \bar{U} +1)^2 - \zeta^2]}{ \zeta \Delta\zeta}  \biggl(P_2 - P_1 \biggr)</math>
  </td>
  <td align="center">
<math>~=</math>
  </td>
  <td align="left">
<math>~(\Rho_1 + \Rho_2) [4 - (\Rho_1 + \Rho_2) (\zeta \bar{U} +1)] </math>
  </td>
</tr>

<tr>
  <td align="right">
&nbsp;
  </td>
  <td align="center">
<math>~=</math>
  </td>
  <td align="left">
<math>~
4(\Rho_1 + \Rho_2)
- (\Rho_1^2 + 2\Rho_1\Rho_2 +\Rho_2^2) (\zeta \bar{U} +1)
</math>
  </td>
</tr>

<tr>
  <td align="right">
&nbsp;
  </td>
  <td align="center">
<math>~=</math>
  </td>
  <td align="left">
<math>~
\Rho_2 [4
- 2\Rho_1 (\zeta \bar{U} +1) ]
- \Rho_2^2 (\zeta \bar{U} +1) 
+ [ 4\Rho_1 - \Rho_1^2  (\zeta \bar{U} +1) ]
</math>
  </td>
</tr>

<tr>
  <td align="right">
<math>~\Rightarrow ~~~ 
4[ (\zeta \bar{U} +1)^2 - \zeta^2] P_2 
</math>
  </td>
  <td align="center">
<math>~=</math>
  </td>
  <td align="left">
<math>~
\Rho_2 [4 - 2\Rho_1 (\zeta \bar{U} +1) ]\zeta \Delta\zeta
- \Rho_2^2 (\zeta \bar{U} +1)\zeta \Delta\zeta 
+ \biggl\{ [ 4\Rho_1 - \Rho_1^2  (\zeta \bar{U} +1) ]\zeta \Delta\zeta
+ 4[ (\zeta \bar{U} +1)^2 - \zeta^2] P_1 \biggr\} \, .
</math>
  </td>
</tr>
</table>
</div>
This is a quadratic  equation of the form,
<div align="center">
<table border="0" cellpadding="5" align="center">

<tr>
  <td align="right">
<math>~aP_2^2 + bP_2 + c </math>
  </td>
  <td align="center">
<math>~=</math>
  </td>
  <td align="left">
<math>~0 \, ,</math>
  </td>
</tr>
</table>
</div>
where,
<div align="center">
<table border="0" cellpadding="5" align="center">

<tr>
  <td align="right">
<math>~a</math>
  </td>
  <td align="center">
<math>~\equiv</math>
  </td>
  <td align="left">
<math>~ 
(\zeta \bar{U} +1)\zeta \Delta\zeta 
= \ell \zeta \Delta\zeta 
\, ,</math>
  </td>
</tr>

<tr>
  <td align="right">
<math>~b</math>
  </td>
  <td align="center">
<math>~\equiv</math>
  </td>
  <td align="left">
<math>~
- \{ [4 - 2\Rho_1 (\zeta \bar{U} +1) ]\zeta \Delta\zeta
- 4[ (\zeta \bar{U} +1)^2 - \zeta^2] \}
=
- [ (4 - 2\Rho_1 \ell )\zeta \Delta\zeta - 4( \ell^2 - \zeta^2) ]
\, ,</math>
  </td>
</tr>

<tr>
  <td align="right">
<math>~c</math>
  </td>
  <td align="center">
<math>~\equiv</math>
  </td>
  <td align="left">
<math>~
- \{ [ 4\Rho_1 - \Rho_1^2  (\zeta \bar{U} +1) ]\zeta \Delta\zeta
+ 4[ (\zeta \bar{U} +1)^2 - \zeta^2] P_1 \} 
=
- [ ( 4\Rho_1 - \Rho_1^2 \ell )\zeta \Delta\zeta + 4( \ell^2 - \zeta^2) P_1 ] 
\, .
</math>
  </td>
</tr>
</table>
</div>

The pair of roots of this equation are, then,

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

<tr>
  <td align="right">
<math>~2 \ell \zeta \Delta\zeta P_2</math>
  </td>
  <td align="center">
<math>~=</math>
  </td>
  <td align="left">
<math>~
[ (4 - 2\Rho_1 \ell )\zeta \Delta\zeta - 4( \ell^2 - \zeta^2) ] 
\pm \biggl\{ 
[ (4 - 2\Rho_1 \ell )\zeta \Delta\zeta - 4( \ell^2 - \zeta^2) ]^2 
+ 4\ell \zeta \Delta\zeta[ ( 4\Rho_1 - \Rho_1^2 \ell )\zeta \Delta\zeta + 4( \ell^2 - \zeta^2) P_1 ] 
\biggr\}^{1 / 2} \, .
</math>
  </td>
</tr>
</table>
</div>
Or,

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

<tr>
  <td align="right">
<math>~\biggl\{ \frac{2 \ell \zeta \Delta\zeta}{ [ (4 - 2\Rho_1 \ell )\zeta \Delta\zeta - 4( \ell^2 - \zeta^2) ] } \biggr\}P_2</math>
  </td>
  <td align="center">
<math>~=</math>
  </td>
  <td align="left">
<math>~
1 
\pm \biggl\{ 1
+ \frac{4\ell \zeta \Delta\zeta[ ( 4\Rho_1 - \Rho_1^2 \ell )\zeta \Delta\zeta + 4( \ell^2 - \zeta^2) P_1 ]}{[ (4 - 2\Rho_1 \ell )\zeta \Delta\zeta - 4( \ell^2 - \zeta^2) ]^2 } 
\biggr\}^{1 / 2} \, .
</math>
  </td>
</tr>
</table>
</div>

==Combined==
Now, before returning to the first ODE, let's write <math>~\Delta U</math> in terms of <math>~\ell</math> and, hereafter, use <math>~\ell</math> as the unknown instead of <math>~U_2</math>.
<div align="center">
<table border="0" cellpadding="5" align="center">

<tr>
  <td align="right">
<math>~\Delta U \equiv U_2 - U_1</math>
  </td>
  <td align="center">
<math>~=</math>
  </td>
  <td align="left">
<math>~(2\bar{U} - U_1) - U_1</math>
  </td>
</tr>

<tr>
  <td align="right">
&nbsp;
  </td>
  <td align="center">
<math>~=</math>
  </td>
  <td align="left">
<math>~2\biggl[\bar{U} - U_1 \biggr]</math>
  </td>
</tr>

<tr>
  <td align="right">
&nbsp;
  </td>
  <td align="center">
<math>~=</math>
  </td>
  <td align="left">
<math>~2\biggl[\frac{(\ell - 1)}{\zeta} - U_1 \biggr]</math>
  </td>
</tr>
</table>
</div>

Hence, the first ODE gives,
<div align="center">
<table border="0" cellpadding="5" align="center">

<tr>
  <td align="right">
<math>~
\ell^2 \Rho_2  
</math>
  </td>
  <td align="center">
<math>~=</math>
  </td>
  <td align="left">
<math>~
2\biggl\{ 2\ell +\biggl[ \ell^2 - \zeta^2 \biggr]\frac{\Delta U}{\Delta\zeta} \biggr\} - \ell^2 P_1 
</math>
  </td>
</tr>

<tr>
  <td align="right">
<math>~\Rightarrow~~~
\Delta\zeta \ell^2 \Rho_2  
</math>
  </td>
  <td align="center">
<math>~=</math>
  </td>
  <td align="left">
<math>~
(4\ell - \ell^2 P_1)\Delta\zeta  + 2( \ell^2 - \zeta^2 ) \Delta U  
</math>
  </td>
</tr>

<tr>
  <td align="right">
&nbsp;
  </td>
  <td align="center">
<math>~=</math>
  </td>
  <td align="left">
<math>~
(4\ell - \ell^2 P_1)\Delta\zeta  + 4( \ell^2 - \zeta^2 ) \biggl[\frac{(\ell - 1)}{\zeta} - U_1 \biggr] 
</math>
  </td>
</tr>

<tr>
  <td align="right">
<math>~\Rightarrow~~~
(2 \zeta \ell^2 \Delta\zeta )\Rho_2  
</math>
  </td>
  <td align="center">
<math>~=</math>
  </td>
  <td align="left">
<math>~
2(4\ell - \ell^2 P_1) \zeta\Delta\zeta  + 8( \ell^2 - \zeta^2 ) (\ell - 1 - \zeta U_1 ) \, . 
</math>
  </td>
</tr>
</table>
</div>

Finally, using this to replace <math>~\Rho_2</math> in the second ODE expression gives,
<div align="center">
<table border="0" cellpadding="5" align="center">

<tr>
  <td align="right">
<math>~
2(4\ell - \ell^2 P_1) \zeta\Delta\zeta  + 8( \ell^2 - \zeta^2 ) (\ell - 1 - \zeta U_1 ) 
</math>
  </td>
  <td align="center">
<math>~=</math>
  </td>
  <td align="left">
<math>~
\ell [ (4 - 2\Rho_1 \ell )\zeta \Delta\zeta - 4( \ell^2 - \zeta^2) ] 
</math>
  </td>
</tr>

<tr>
  <td align="right">
&nbsp;
  </td>
  <td align="center">
&nbsp;
  </td>
  <td align="left">
<math>~
\pm \ell \biggl\{ 
[ (4 - 2\Rho_1 \ell )\zeta \Delta\zeta - 4( \ell^2 - \zeta^2) ]^2 
+ 4\ell \zeta \Delta\zeta[ ( 4\Rho_1 - \Rho_1^2 \ell )\zeta \Delta\zeta + 4( \ell^2 - \zeta^2) P_1 ] 
\biggr\}^{1 / 2} 
</math>
  </td>
</tr>

<tr>
  <td align="right">
<math>~\Rightarrow~~~
\frac{8}{\ell} ( \ell^2 - \zeta^2 ) (\ell - 1 - \zeta U_1 ) 
</math>
  </td>
  <td align="center">
<math>~=</math>
  </td>
  <td align="left">
<math>~
-4 \zeta \Delta\zeta   - 4( \ell^2 - \zeta^2)
</math>
  </td>
</tr>

<tr>
  <td align="right">
&nbsp;
  </td>
  <td align="center">
&nbsp;
  </td>
  <td align="left">
<math>~
\pm \biggl\{ 
\biggl[ (4 - 2\Rho_1 \ell )\zeta \Delta\zeta - 4( \ell^2 - \zeta^2) \biggr]^2 
+ 4\ell \zeta \Delta\zeta \biggl[ ( 4\Rho_1 - \Rho_1^2 \ell )\zeta \Delta\zeta + 4( \ell^2 - \zeta^2) P_1 \biggr] 
\biggr\}^{1 / 2} 
</math>
  </td>
</tr>

<tr>
  <td align="right">
<math>~\Rightarrow~~~
\frac{1}{\ell}  ( \ell^2 - \zeta^2 )\biggl[ 2  (\ell - 1 - \zeta U_1 ) + \ell \biggr] + \zeta \Delta\zeta
</math>
  </td>
  <td align="center">
<math>~=</math>
  </td>
  <td align="left">
<math>~
\pm \frac{1}{4} \biggl\{ 
\biggl[ (4 - 2\Rho_1 \ell )\zeta \Delta\zeta - 4( \ell^2 - \zeta^2) \biggr]^2 
+ 4\ell \zeta \Delta\zeta \biggl[ ( 4\Rho_1 - \Rho_1^2 \ell )\zeta \Delta\zeta + 4( \ell^2 - \zeta^2) P_1 \biggr] 
\biggr\}^{1 / 2} 
</math>
  </td>
</tr>

<tr>
  <td align="right">
<math>~\Rightarrow~~~
( \ell^2 - \zeta^2 )\biggl[ 3\ell - 2 - 2\zeta U_1 \biggr] + \ell \zeta \Delta\zeta
</math>
  </td>
  <td align="center">
<math>~=</math>
  </td>
  <td align="left">
<math>~
\pm \frac{\ell}{4} \biggl\{ 
\biggl[ (4 - 2\Rho_1 \ell )\zeta \Delta\zeta - 4( \ell^2 - \zeta^2) \biggr]^2 
+ 4\ell \zeta \Delta\zeta \biggl[ ( 4\Rho_1 - \Rho_1^2 \ell )\zeta \Delta\zeta + 4( \ell^2 - \zeta^2) P_1 \biggr] 
\biggr\}^{1 / 2} \, .
</math>
  </td>
</tr>
</table>
</div>
This doesn't look particularly useful because, after squaring both sides, it is a sixth-order polynomial in <math>~\ell</math>, which generally has no analytic solution.