Notebook

これは日々の作業を通して学んだことや毎日の生活で気づいたことをを記録しておく備忘録である。

HTML ファイル生成日時: 2024/05/18 17:09:42.458 (台灣標準時)

Lane-Emden 方程式を数値的に解く

Scipy の scipy.integrate.solve_ivp を使って Lane-Emden 方程式を数値的 に解いてみたでござる。コードは以下の通りにござる。同じコードを GitHub にも置いておいたでござる。


#!/usr/pkg/bin/python3.9

#
# Time-stamp: <2022/05/02 20:15:11 (CST) daisuke>
#

# importing argparse module
import argparse

# importing sys module
import sys

# importing pathlib module
import pathlib

# importing numpy module
import numpy

# importing scipy module
import scipy.integrate

# initialising a parser
desc   = 'solving Lane-Emden equation numerically'
parser = argparse.ArgumentParser (description=desc)

# adding arguments
parser.add_argument ('-n', '--n', type=float, default=0.0,
                     help='index n of polytrope (default: 0.0)')
parser.add_argument ('-s', '--step', type=float, default=0.00001,
                     help='step size of radius (default: 0.00001)')
parser.add_argument ('-o', '--output', default='', \
                     help='output file name')

# parsing arguments
args = parser.parse_args ()

# parameters
n           = args.n
step        = args.step
file_output = args.output

# existence check of output file
path_output = pathlib.Path (file_output)
if (path_output.exists ()):
    # printing message
    print ("ERROR: output file '%s' exists, exiting...")
    # exit
    sys.exit ()

# check of polytropic index n
if (n > 4.95):
    # printing message
    print ("ERROR: choose polytropic index between 0 and 4.95, exiting...")
    # exit
    sys.exit ()

#
# Lane-Emden equation
#
def laneemden (t, y):
    # making a list with a dimension same as y
    dy = numpy.zeros_like (y)
    # Lane-Emden equation
    dy[0] = y[1]
    dy[1] = -y[0]**n - 2.0 * y[1] / t
    # returning values
    return dy

#
# event trigger
#
def cross_zero (t, y):
    return y[0]
cross_zero.terminal  = True
cross_zero.direction = -1

# initial values
y_init  = [ 1.0, 0.0 ]

# output dimensionless radius values
xi_max = 350.0
n_step = int (xi_max / step)
t_eval = numpy.linspace (step, xi_max, n_step)

# solving equation
sol = scipy.integrate.solve_ivp (laneemden, [step, xi_max], y_init, \
                                 t_eval=t_eval, dense_output=True, \
                                 events=cross_zero, \
                                 rtol=10**-12, atol=10**-15)

# quick-and-dirty job for finding the surface
for i in range ( len (sol.y[0]) ):
    if (sol.y[0][i] > 0.0):
        xi_surface = (i+1)*step
    else:
        break
        
with open (file_output, 'w') as fh:
    # writing a header
    header = "#\n"
    header += "# xi, theta, r_over_R, rho_over_rho_c, P_over_P_c\n"
    header += "#\n"
    header += "# polytropic index n = %f\n" % n
    header += "# xi_surface = %f\n" % xi_surface
    header += "#\n"
    fh.write (header)
    
    # printing results
    for i in range ( len (sol.y[0]) ):
        # xi
        xi = (i + 1) * step
        # theta
        theta = sol.y[0][i]
        # stop writing data if theta is negative
        if (theta < 0.0):
            break
        # normalised radius
        r_over_R = xi / xi_surface
        # normalised density
        rho_over_rho_c = theta**n
        # normalised pressure
        P_over_P_c = theta**(n+1)
        # output data
        data = "%15.13f %15.13f %15.13f %15.13f %15.13f\n" \
            % (xi, theta, r_over_R, rho_over_rho_c, P_over_P_c)
        # writing data to output file
        fh.write (data)

例えば、 n=3 の場合を計算させるには、以下のようにすればよいでござる。


% ./laneemden.py -n 3.0 -o laneemden_n30.data

幾つかの polytropic index について計算して、その結果を図にしてみたでご ざる。

fig_202205/polytrope_temperature.png
fig_202205/polytrope_density.png
fig_202205/polytrope_pressure.png


Frequently accessed files

  1. Computer___Python/20220518_0.html
  2. Computer___FreeBSD/20220621_0.html
  3. Computer___Network/20230516_00.html
  4. Food___Taiwan/20220429_0.html
  5. Computer___Python/20220715_0.html
  6. Computer___NetBSD/20220817_3.html
  7. Computer___Python/20220410_0.html
  8. Computer___Debian/20210223_1.html
  9. Food___Taiwan/20230526_00.html
  10. Food___Taiwan/20210205_5.html
  11. Computer___Network/20230508_00.html
  12. Computer___Python/20210124_0.html
  13. Computer___NetBSD/20220818_1.html
  14. Computer___NetBSD/20220428_0.html
  15. Food___Taiwan/20210205_1.html
  16. Computer___Network/20230726_00.html
  17. Science___Astronomy/20220503_0.html
  18. Computer___TeX/20230503_00.html
  19. Travel___Taiwan/20220809_2.html
  20. Science___Math/20220420_0.html
  21. Computer___NetBSD/20210204_0.html
  22. Computer___NetBSD/20230515_00.html
  23. Food___Taiwan/20230618_00.html
  24. Computer___NetBSD/20220808_0.html
  25. Science___Astronomy/20220420_1.html
  26. Computer___NetBSD/20230119_00.html
  27. Computer___Python/20220816_1.html
  28. Food___Taiwan/20220424_8.html
  29. Computer___Python/20221013_0.html
  30. Food___Taiwan/20220424_1.html


HTML file generated by Kinoshita Daisuke.