代码之家 › 专栏 › 技术社区 › Nickpick

Sympy与天体物理学相结合,用单位进行数值求解?

astropy solver sympy python

Nickpick · 技术社区 · 4 年前

有没有一种方法可以将Sympy的求解器与Astropy结合起来?例如,如果我想在所有其他变量的给定值下求解p的开普勒第三定律:

from sympy import solveset, Poly, Eq, Function, exp
from sympy import Symbol, pi
Ms = Symbol('Ms') # mass of sun
Mp = Symbol('Mp') # mass of planet
G = Symbol('G') # Gravitational consatant
a = Symbol('a') # semi jamor axes
P = Symbol('P') # period
solveset(G * (Ms + Mp) / 4 * pi**2 - a**3/P**2, a)

然后,我如何使用占星术将具有相应单位的值添加到所有变量中,并获得数值输出?

import astropy.units as u
from astropy.constants import G
Ms = 1 * u.M_sun
Mp = 1 * u.M_jup
s = 5 * u.pc
P = 1 * u.yr

这可以与Sympy收到的结果相结合,以求解具有预期单位和数值结果的a中的值吗?

任何建议都将不胜感激。

0 回复 | 直到 4 年前

Matt Pitkin 3 年前

我一直在玩用 lambdify 能够在Sympy方程中使用天体测量单位来实现。以下是我为这个问题提出的解决方案:

from sympy import pi, Symbol, lambdify, sympify, symbols, solve

import astropy.units as u
from astropy.constants import G as astroG

Ms = Symbol('Ms') # mass of sun
Mp = Symbol('Mp') # mass of planet
G = Symbol('G') # Gravitational consatant
a = Symbol('a') # semi major axes
P = Symbol('P') # period
solutions = solve(
    (G * (Ms + Mp) / (4 * pi**2)) - (a**3/P**2),
    a
)

# set values at which to evaluate the solutions
values = {
    "P": 2.1 * u.yr,
    "Ms": 1.4 * u.M_sun,
    "Mp": 0.5 * u.M_jup, # (0.5 * u.M_jup.to("M_sun")),
    "G": astroG,
}

evalsolutions = []

reqsymbols = [Ms, Mp, P, G]

for sol in solutions:
    # lambdify solutions
    f = lambdify(
        reqsymbols,
        sol,
        modules=["numpy"],
    )
    
    res = f(**values)
    if res.dtype != "complex":
        evalsolutions.append(res.si.decompose())

print(evalsolutions)

这给出了:

[<Quantity 2.74467861e+11 m>]

旧解决方案

这是我最初发布的另一个解决方案。它更复杂,而且被证明是矫枉过正的,但确实展示了如何利用 modules lambdify-ed函数的参数。

from sympy import pi, Symbol, lambdify, sympify, symbols, solve

import astropy.units as u
from astropy.constants import G as astroG

Ms = Symbol('Ms') # mass of sun
Mp = Symbol('Mp') # mass of planet
Munit = sympify('unit(msol)')
G = sympify('const(G)') # Gravitational consatant
a = Symbol('a') # semi major axes
P = Symbol('P') # period
Punit = Symbol('unit(yr)')
solutions = solve(
    (G * ((Ms + Mp) * Munit) / (4 * pi**2)) - (a**3/(P * Punit)**2),
    a
)

UNITS = {
    "msol": u.M_sun,
    "au": u.au,
    "yr": u.yr,
}

CONSTANTS = {
    "G": astroG,
}

def unitfunc(key):
    return 1.0 * UNITS.get(key, 1)

def constfunc(key):
    return CONSTANTS.get(key, 1)

# set values at which to evaluate the solutions
values = {
    "P": 2.1,
    "Ms": 1.4,
    "Mp": (0.5 * u.M_jup.to("M_sun")),
}

values.update({key: key for key in UNITS})
values.update({key: key for key in CONSTANTS})

evalsolutions = []

reqsymbols = [Ms, Mp, P]
reqsymbols.extend([symbols(key) for key in list(UNITS.keys()) + list(CONSTANTS.keys())])

for sol in solutions:
    # lambdify solutions
    f = lambdify(
        reqsymbols,
        sol,
        modules=[{"unit": unitfunc, "const": constfunc}, "numpy"],
    )
    
    res = f(**values)
    if res.dtype != "complex":
        evalsolutions.append(res.si.decompose())

print(evalsolutions)