python - Pythonを使用して非線形方程式を解く

Question

私はPythonを使ったことがありませんが、Mathematicaは私が解こうとしている方程式を処理できません。次の方程式の変数「a」を解こうとしています。ここで、s、c、mu、およびdeltatは既知のパラメーターです。

MathematicaでNSolveやSolveなどを試してみましたが、1時間実行されて運がありませんでした。私はPythonに精通していないので、Pythonを使用してこの方程式を解く方法はありますか？

Answer 1

aこれらの方程式は超越的であり、三角関数の内側と外側の両方を含んでいるため、これらの方程式の解析解を見つけることはできません。

数値解法で問題になっているのは、の許容値の範囲がaによって制約されていることだと思いますarcsin。は-1と1の間の引数に対してのみ定義されているため（実際にarcsinなりたいと仮定）、の式はそれとを必要とします。aalphabetaa > s/2a > (s-c)/2

Pythonでは、次の関数f(a) = 0を使用して、3番目の方程式（形式で書き直された）のゼロを見つけることができます。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)でゼロを検索することです。ここでは、それが少なくともであることがわかります。私はちょうどそれがもっともらしい上限であり、与えられたパラメータに対して機能したと推測しました。より一般的には、変更がとの間で署名されることを確認する必要があります。関数がどのように機能するかについて詳しくは、SciPyのドキュメントをご覧ください。f[a,b]amax(s/2, (s-c)/2)fab

Answer 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) )

# <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での以前の解は、アークサインを動作させるために必要なクリッピングによるものです。実際、この方程式はa=5000以下では意味がありません。（正確な値は、レイズソリューションで指定されたa0です）。これにより、Raysソリューションの手法で使用できる優れた範囲が得られます。