In-Class Assignment 1

In-Class Assignment 1#

Credit: These exercises follow from work by Mike Zingale’s computational lectures. Used here with explicit permission.

I’m using SymPy – the symbolic math library for python – to do the algebra here. The next 2 cells load SymPy for interactive work and define the variables we are going to use (and treat as math symbols)

Let’s start by importing the necessary components:

from sympy import init_session
init_session(use_latex="mathjax")
%matplotlib inline

IPython console for SymPy 1.13.3 (Python 3.11.8-64-bit) (ground types: python)

These commands were executed:
>>> from sympy import *
>>> x, y, z, t = symbols('x y z t')
>>> k, m, n = symbols('k m n', integer=True)
>>> f, g, h = symbols('f g h', cls=Function)
>>> init_printing()

Documentation can be found at https://docs.sympy.org/1.13.3/

1. stellar model#

Below are some examples of defining variables that we will use in the follow exercises.

rho = symbols('rho', cls=Function)
rhoc = symbols('rho_c')
qc = symbols('q_c')
Pc = symbols('P_c')
G = symbols('G')
Mstar, Rstar = symbols('M_\star R_\star')
r = symbols('r')
xi = symbols('xi')
mu = symbols('mu')
k = symbols('k')
m_u = symbols('m_u')

Imagine a star has a density profile of the form:

\[\rho(r) = \rho_c \left [ 1 - \left (\frac{r}{R_\star} \right )^3 \right ]\]

where \(\rho_c\) is the central density.

a.#

first, plot a new quantity, lets call it \(z = \rho (\xi * R_{\star})/\rho_c \), where we have taken the equation above, normalized by the still unknown central density and substituted \(r\).

https://docs.sympy.org/latest/modules/core.html#sympy.core.basic.Basic.subs

# plot z as a function of xi from 0 to 1 below

b.#

Find an expression for the central density in terms of \(R_\star\) and \(M_\star\).

You need to use the mass equation from Chapter 1.1 in HKT.

# define rho here as a function of r, not normalized

# Next, we have to compute M_star here

https://docs.sympy.org/latest/modules/core.html#sympy.core.expr.Expr.integrate

# Now, use M_star to finally find the central density

Now we will symbolically replace \(\rho_c\) in the density expression with this result for the future calculations

# Do substitution here for new rho(r) expression

c.#

Use HSE and the fact that the pressure vanishes at the surface to find the pressure as a function of radius. Your answer will be in the form of \(P(r) = P_c \times\) (polynomial in \(r/R_\star\)).

What is \(P_c\) in terms of \(M_\star\) and \(R_\star\)?

Express \(P_c\) numerically with \(M_\star\) and \(R_\star\) in solar units.

We need to find a general expression for \(M(r)\)

# plot the resulting expression here for Pc here

Now we can integrate HSE. We will do this as

\[P = P_c + \int_0^{R_\star} \frac{-GM(r)}{r^2} \rho dr\]

We’ll substitute in our expressions for \(M(r)\) and \(\rho(r)\) in the integrand.

# Solve for Pc here

However, recall that \(P(R_\star) = 0\), so we can enforce that here to find \(P_c\)

# Solve here given this boundary condition, did this change the result of the integrand?

For our general pressure expression. we integrate to some shell \(r\)

# Compute our general expression here for P

Or factoring \(P_c\) out (i.e. \(P/P_c\)), we have

# Simplify the expression by factoring out P_c

The resulting equation should be

\begin{equation} P(r) = P_c \left [1 - \frac{40}{21} \left (\frac{r}{R_\star}\right )^2 + \frac{8}{7} \left (\frac{r}{R_\star}\right )^3 - \frac{5}{21} \left ( \frac{r}{R_\star}\right )^8 \right ] \end{equation}

# Finally, plot a rescaled version of the result for P, as we did in a.
Contents