你不会找到这些方程的解析解,因为它们是超越的,包含三角函数的内部和外部.
我认为你在数值解决方案中遇到的麻烦是a的可接受值范围受到arcsin的约束.由于arcsin仅为-1和1之间的参数定义(假设你想要一个真实的),你的alpha和beta公式要求a> s / 2和a> (S-C)/ 2.
在Python中,您可以使用brentq函数找到第三个等式的零(以f(a)= 0形式重写):
import numpy as np
from scipy.optimize import brentq
s = 10014.6
c = 6339.06
mu = 398600.0
dt = 780.0
def f(a):
alpha = 2*np.arcsin(np.sqrt(s/(2*a)))
beta = 2*np.arcsin(np.sqrt((s-c)/(2*a)))
return alpha - beta - (np.sin(alpha)-np.sin(beta)) - np.sqrt(mu/a**3)*dt
a0 = max(s/2, (s-c)/2)
a = brentq(f, a0, 10*a0)
编辑:
澄清brentq(f,a,b)的工作方式是它在区间[a,b]上搜索f的零.在这里,我们知道a至少是max(s / 2,(s-c)/ 2).我只是猜测这是一个看似合理的上限的10倍,并且对于给定的参数起作用.更一般地说,您需要确保f在a和b之间改变符号.您可以在SciPy docs中详细了解该功能的工作原理.