Homework 4 solutions

We want to invert and equation of state (EOS), i.e., given a pressure, \(p_\star\), and specific internal energy, \(e_\star\), and an equation of state that is expressed in terms of density and temperature:

\[\begin{align*} p &= p(\rho, T) \\ e &= e(\rho, T) \end{align*}\]

we want to find the \(\rho\) and \(T\) that give our desired pressure and energy.

First we’ll write our equation of state. We’ll create a class, EOSState, that simply serves as a type to hold all of the thermodynamic information for our state.

class EOSState:
    def __init__(self, rho=None, T=None,
                 p=None, e=None,
                 dpdrho=None, dpdT=None,
                 dedrho=None, dedT=None, mu=4):
        self.rho = rho
        self.T = T
        self.mu = mu
        
        self.p = p
        self.e = e
        self.dpdrho = dpdrho
        self.dpdT = dpdT
        self.dedrho = dedrho
        self.dedT = dedT

and now our EOS—it simply takes an EOSState and will fill the components based on the density and temperature in the state.

As stated in the assignment, this is simply an ideal gas + radiation. We can compute all of the thermodynamic derivatives analytically.

def eos(state):
    # constants
    k = 1.38e-16  # erg/K
    a = 7.56e-15  # erg/cm^3/K^4
    m_u = 1.66e-24  # g

    state.p = (state.rho * k * state.T / (state.mu * m_u) +
               (1./3.) * a * state.T**4)
    state.e = 1.5 * k * state.T / (state.mu * m_u) + a * state.T**4 / state.rho

    state.dpdrho = k * state.T / (state.mu * m_u)
    state.dpdT = state.rho * k / (state.mu * m_u) + (4./3.) * a * state.T**3

    state.dedrho = -a * state.T**4 / state.rho**2
    state.dedT = 1.5 * k / (state.mu * m_u) + 4 * a * state.T**3 / state.rho

We use a Newton’s method for our system, writing it as:

\[\begin{split}\Psi(\rho_0 + \delta_\rho, T_0 + \delta T) = 0 = \Psi(\rho_0, T_0) + {\bf J} \left ( \begin{array}{c} \delta \rho \\ \delta T \end{array} \right ) + \ldots\end{split}\]

where \((\rho_0, T_0)\) is the initial guess for the density and temperature that satisfy our thermodynamics.

We find the corrections, \((\delta\rho, \delta T)\) by solving the linear system:

\[\begin{split}{\bf J} \left (\begin{array}{c} \delta \rho \\ \delta T \end{array} \right ) = - \Psi(\rho_0, T_0)\end{split}\]

and iterating until the corrections are small. Here, \({\bf J}\) is the Jacobian, which has the form:

\[\begin{split}{\bf J} = \left ( \begin{array}{cc} \partial p/\partial \rho |_T & \partial p/\partial T |_\rho \\ \partial e/\partial \rho |_T & \partial e/\partial T |_\rho \end{array} \right )\end{split}\]

import numpy as np

Here’s our implementation—we pass in the pressure and energy we want, and optionally a guess for \(\rho_0\) and \(T_0\).

def rhoT_from_pe(p_in, e_in, *, rho0=1, T0=1, tol=1.e-5):
    state = EOSState(rho=rho0, T=T0)

    err = 1.e30
    while (err > tol):
        # get the current thermodynamics
        eos(state)

        # construct the Jacobian and Psi
        J = np.array([[state.dpdrho, state.dpdT],
                      [state.dedrho, state.dedT]])
        psi = np.array([state.p - p_in, state.e - e_in])

        # solve for the corrections
        delta = np.linalg.solve(J, -psi)

        # update our guesses
        state.rho += delta[0]
        state.T += delta[1]

        # compute the error
        err = max(abs(delta[0]/state.rho), abs(delta[1]/state.T))

    return state.rho, state.T

Now we can test it out.

p_star = 2.3e10
e_star = 3.87e13

rho, T = rhoT_from_pe(p_star, e_star, rho0=1.e2, T0=1.e5)
(rho, T)

(np.float64(0.000988366164960824), np.float64(997994.6760719108))

Now we can check how well we did by calling the EOS with the \(\rho\) and \(T\) we found and comparing to the pressure and energy we wanted.

s_new = EOSState(rho=rho, T=T)
eos(s_new)

print(f"pressure: we found {s_new.p} and wanted {p_star}")
print(f"energy: we found {s_new.e} and wanted {e_star}")

pressure: we found 23000000000.01443 and wanted 23000000000.0
energy: we found 38700000000007.32 and wanted 38700000000000.0

Clearly we did quite well.