Using a Polytrope To Mathematically Model Stellar Structure

Version 1.0 of 12/9/2004-1:00 a.m.

Let's see if we can model the structure of the sun, using some interesting mathematics. We start by using the polytrope model [1], so we use Maple to get some interesting results. Let's load some parameters:

> restart;
> n:=3; #polytrope index for sun
> dc:=160000; #density at the core (kg/m^{3})
> Pc:=10^(16.53); #pressure at the core (Pascal)
> G:=6.67428*10^(-11); #gravitational constant (m^{3}/kg/s^{2})
> r:=6.955*10^8; #radius of the sun (m)
> zeta:=r*(4*Pi*G*dc^2/(n+1)/Pc)^(1/2); #run extent for solution
> evalf(%);

Let's now build the fundamental differential equation for the model, the Lane-Emden equation [2/3]:

Now we calculate the density and pressure as a function of radius:

> d:=x->dc*y(x)^n; #density as a function of radius
> K:=Pc/dc^(1+1/n); #calculate constant K
> P:=x->K*d(x)^(1+1/n); #pressure as a function of radius
> rext:=fsolve(d(x)=0,x=0..zeta); #find where density drops to 0, in [0,zeta]

We now normalize the radius with respect to rext and what do we get?

> with(plots):
> plot(d(R*rext),R=0..1);

Density of the sun as a function of radius

> plot(P(R*rext),R=0..1);

Pressure of the sun as a function of radius

The first graph agrees fairly well with the data on the graph of Wikipedia's polytrope page [1], and with [4], [5] and [6], which give the density of the sun as a function of its radius.