我从来没有使用蟒蛇,但数学不能处理我试图求解方程。 我试图解决对于其中S,C,μ,和delta吨是已知的参数下面的等式的变量“a”。
我试图做NSolve,解决等在数学,但它已经运行了,没有运气小时。 由于我不熟悉Python,是有办法,我可以使用Python来解决这个方程的?
我从来没有使用蟒蛇,但数学不能处理我试图求解方程。 我试图解决对于其中S,C,μ,和delta吨是已知的参数下面的等式的变量“a”。
我试图做NSolve,解决等在数学,但它已经运行了,没有运气小时。 由于我不熟悉Python,是有办法,我可以使用Python来解决这个方程的?
你不会找到一个解析解这些方程,因为他们是超然的,包含a
内部和三角函数之外。
我想你和数值解具有麻烦的是,可接受值的范围a
被约束arcsin
。 由于arcsin
为参数-1到1之间只定义了(假设你想a
是真实的),你的公式alpha
和beta
需要a > s/2
和a > (sc)/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)
的工作原理是,它搜索的零f
上的间隔[a,b]
在这里,我们知道, a
至少是max(s/2, (sc)/2)
我只是猜测,10倍,这是一个合理的上限,这工作给定参数。 更一般地,你需要确保f
变化之间签订a
和b
。 你可以阅读更多有关功能的工作原理SciPy的文档 。
我认为它值得研究atempting解决之前,该函数的行为。 如果不这样做,你不知道如果有一个独特的解决方案,许多解决方案,还是无解。 (最大的问题是许多解决方案,其中数值方法可能不给你你需要/期望的解决方案 - 如果你一味地使用它“不好的事情”可能发生)。 你检查的行为很好地使用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的行为所需的剪辑。 真正的方程式无厘头低于约= 5000。 (精确值是在射线溶液给出的A0)。 然后,这给出了可与射线溶液的技术一起使用一个很好的范围。