2

我从未使用过 python,但 Mathematica 无法处理我试图求解的方程。我正在尝试求解以下方程的变量“a”,其中 s、c、mu 和 delta t 是已知参数。

在此处输入图像描述

我尝试在 Mathematica 中进行 NSolve、Solve 等,但它已经运行了一个小时,但没有运气。由于我不熟悉 Python,有没有办法可以使用 Python 来解决这个方程?

4

2 回答 2

3

你不会找到这些方程的解析解,因为它们是超越的,包含a三角函数的内部和外部。

我认为您在使用数值解决方案时遇到的问题是可接受值的范围aarcsin. 由于arcsin仅为 -1 和 1 之间的参数定义(假设您想要a是真实的),因此您的公式alphabeta需要a > s/2and 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]amax(s/2, (s-c)/2)fab

于 2013-03-27T01:37:54.403 回答
2

我认为在尝试解决它之前检查函数的行为是值得的。如果不这样做,您将不知道是否有唯一的解决方案、许多解决方案或没有解决方案。(最大的问题是许多解决方案,其中数值方法可能无法为您提供您需要/期望的解决方案 - 如果您盲目使用它,“坏事”可能会发生)。您使用 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) )

fn 图

# <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 解决方案中的技术一起使用。

于 2013-03-27T04:59:57.153 回答