Self-consistent temperatures with Newton's Method

[mabruzzo]

I just wanted to jot down my notes for how we could use Newton's method to self-consistently solve for the gas temperature (rather than fixed point iteration).

While the math is a little ugly, it turns out that using Newton's method with analytic derivatives is relatively cheap (from a computational perspective).

In that scenario we would be solving for the temperature for which $f(T) = 0$, where $f(T) = e(T, \mu, \eta_{H2}) - e_{\rm cur}$. To further elaborate:

$$\begin{align} e(T, \mu, \eta_{\rm H2}) &= \frac{k_B T}{\mu m_H} (\gamma - 1)^{-1} \\\ \eta_{H2} &= n_{H2} / (n_{H2} + n_{\rm other})\\\ (\gamma - 1)^{-1} &= (1 - \eta_{H2}) (\gamma_{\rm other} - 1)^{-1} + \eta_{H2} (\gamma_{H2}-1)^{-1} \\\ &= (1 - \eta_{H2}) (\gamma_{\rm other} - 1)^{-1} + 2.5 \eta_{H2} + \eta_{H2} h(x) \\\ h(x) &= exp(x) \left(\frac{x}{\rm exp(x) - 1}\right)^2 \\\ x &= (6100\ {\rm K}) / T \end{align}$$

Note that $f(x)$ is 1 for $x=0$ (i.e. infinite temperature) and asymptotes to 0 as $x$ approaches infinity (i.e. T goes to zero).

Newton's method requires us to know $d f(T) / dT = d e / d T$. It turns out that

$$\begin{align} \frac{d e}{dT} &= \frac{k_B }{\mu m_H} \left( (\gamma - 1)^{-1} + \eta_{H2} T \frac{d h}{dT} \right)\\\ \frac{d x}{dT} &= x^2 / (-6100\ {\rm K})\\\ \frac{d h}{dT} &= \frac{d h}{dx}\, \frac{d x}{dT}\\\ &= \frac{x + 2 + (x - 2) \exp(x)}{-6100} \left(\frac{x}{\rm exp(x) - 1}\right) h(x) \end{align}$$

putting this together you end up with

$$\frac{d e}{dT} =\frac{k_B }{\mu m_H} \left( (\gamma - 1)^{-1} - \eta_{H2} \frac{x + 2 + (x - 2) \exp(x)}{\exp(x) - 1} h(x) \right)$$

It's worth highlighting that if you're computing $e(T, \mu, \eta_{\rm H2})$, then you just need a few extra multiplications to also get $\frac{d e}{dT}$ since you can reuse the values of $k_B/(\mu mH)$, $ \eta_{H2}$, $(\gamma - 1)^{-1}$, $x$, $\exp(x)$ $(\exp(x)-1)^{-1}$, and $h(x)$.

Since a newton-raphson solver requires an extra division to get the next guessed root compared to fixed point iteration, I would guess that a single newton-raphson iteration takes no more than ~4/3 the time of a single fixed point iteration.¹ Thus, you start to save time compared to fixed-point iteration if the newton raphson solver converges in under 3/4 of the steps required for fixed point iteration (considering that newton-raphson is 2nd order I bet this is a common occurrence).

Being able to robustly compute $de/dT$ would also be useful in some other parts of the code (currently, when we need that quantity, we make large approximations)

---

If we are worried about numerical errors for extreme values of $x$, we can rewrite $h(x)$ as $(x/2)^2 {\rm csch}^2 (x/2)$ and, according to mathematica, $dh/dx = -x ((x/2) {\rm coth}(x/2) - 1) {\rm csch}^2(x/2)$ (the derivative should be double-checked).

Then, we could make use of the various hyperbolic trigonometric identities to cheaply compute $dh/dT$. From a quick glance at the wikipedia page, there seems to be a few options (there may be something even better):

• I think we just need to compute ${\rm coth}(x/2)$ since ${\rm coth}^2 (x) - 1 = {\rm csch}^2(x)$
• alternatively, via the tangent half argument formulas, we could rewrite all hyperbolic functions in terms of ${\rm tanh}(x/4)$ (which we could get with std::tanh)

Footnotes

1. This is a conservative guess. I bet that a single newton-raphson iteration takes more like 5/4 or 6/5 of the time of a fixed point iteration... ↩