我从未使用过 python,但 Mathematica 无法处理我试图求解的方程。我正在尝试求解以下方程的变量“a”,其中 s、c、mu 和 delta t 是已知参数。
我尝试在 Mathematica 中进行 NSolve、Solve 等,但它已经运行了一个小时,但没有运气。由于我不熟悉 Python,有没有办法可以使用 Python 来解决这个方程?
我从未使用过 python,但 Mathematica 无法处理我试图求解的方程。我正在尝试求解以下方程的变量“a”,其中 s、c、mu 和 delta t 是已知参数。
我尝试在 Mathematica 中进行 NSolve、Solve 等,但它已经运行了一个小时,但没有运气。由于我不熟悉 Python,有没有办法可以使用 Python 来解决这个方程?
你不会找到这些方程的解析解,因为它们是超越的,包含a
三角函数的内部和外部。
我认为您在使用数值解决方案时遇到的问题是可接受值的范围a
受arcsin
. 由于arcsin
仅为 -1 和 1 之间的参数定义(假设您想要a
是真实的),因此您的公式alpha
和beta
需要a > s/2
and a > (s-c)/2
。
在 Python 中,您可以f(a) = 0
使用以下brentq
函数找到第三个方程的零(以 形式重写):
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)
搜索零。在这里,我们知道至少是. 我只是猜测 10 次是一个合理的上限,并且适用于给定的参数。更一般地说,您需要确保和之间的更改符号。您可以在SciPy 文档中阅读有关该函数如何工作的更多信息。f
[a,b]
a
max(s/2, (s-c)/2)
f
a
b
我认为在尝试解决它之前检查函数的行为是值得的。如果不这样做,您将不知道是否有唯一的解决方案、许多解决方案或没有解决方案。(最大的问题是许多解决方案,其中数值方法可能无法为您提供您需要/期望的解决方案 - 如果您盲目使用它,“坏事”可能会发生)。您使用 scipy 和 ipython 很好地检查了该行为。这是一个示例笔记本
# -*- coding: utf-8 -*-
# <nbformat>3.0</nbformat>
# <codecell>
s = 10014.6
c = 6339.06
mu = 398600.0
dt = 780.0
# <codecell>
def sin_alpha_2(x):
return numpy.sqrt(s/(2*x))
def sin_beta_2(x):
return numpy.sqrt((s-c)/(2*x))
def alpha(x):
return 2*numpy.arcsin( numpy.clip(sin_alpha_2(x),-0.99,0.99) )
def beta(x):
return 2*numpy.arcsin( numpy.clip(sin_beta_2(x),-0.99,0.99) )
# <codecell>
def fn(x):
return alpha(x)-beta(x)-numpy.sin(alpha(x))+numpy.sin(beta(x)) - dt * numpy.sqrt( mu / numpy.power(x,3) )
# <codecell>
xx = numpy.arange(1,20000)
pylab.plot(xx, numpy.clip(fn(xx),-2,2) )
# <codecell>
xx=numpy.arange(4000,10000)
pylab.plot(xx,fn(xx))
# <codecell>
xx=numpy.arange(8000,9000)
pylab.plot(xx,fn(xx))
这表明我们希望找到一个介于 8000 和 9000 之间的解。曲线中大约 5000 处的奇怪扭结和大约 4000 处的较早解是由于使 arcsin 行为所需的剪裁。实际上,方程在 a=5000 以下是没有意义的。(精确值是 Rays 解决方案中给出的 a0)。然后,这提供了一个很好的范围,可以与 Rays 解决方案中的技术一起使用。