6

編集：方程式を取得した参照には、いくつかのエラーが含まれていました。ここで修正しました。解決策は、今では実際に理にかなっているかもしれません!

2 層の流体が地形上を流れる場合、流体内の流速と波の速度の相対的なサイズに応じて、さまざまな解が存在します。

クリティカルフロー

これらは「超臨界」、「亜臨界」、「臨界」と呼ばれます (ここでは最初の 2 つを「超臨界」と呼びます)。

次の方程式は、(h, U0) パラメーター空間での臨界動作と超臨界動作の間の境界線を定義します。

式1

eq2

を消去してd_1c(つまり、それが何であるかは気にしません)、これらの方程式の解を で見つけたいと思い(h, U_0)ます。

単純化要因:

  • 与えられた答えだけが必要です d_0
  • 正確な解は必要ありません。解曲線の概要だけが必要なので、これは解析的または数値的に解決できます。
  • 領域 (h, U0) = (0,0) から (0.5, 1) だけをプロットしたいだけです。

Enthoughtディストリビューション（numpy、scipy、sympy）で利用可能なモジュールを使用してこれを解決したいのですが、どこから始めればよいのか本当にわかりません。私を本当に混乱させるのは、変数 d1c の削除です。

Python の方程式は次のとおりです。

def eq1(h, U0, d1c, d0=0.1):
    f = (U0) ** 2 * ((d0 ** 2 / d1c ** 3) + (1 - d0) ** 2 / (1 - d1c - d0) ** 3) - 1
    return f

def eq2(h, U0, d1c, d0=0.1):
    f = 0.5 * (U0) ** 2 * ((d0 ** 2 / d1c ** 2) - (1 - d0) ** 2 / (1 - d1c - d0) ** 2) + d1c + (h - d_0)
    return f

私は、多くのソリューション ブランチ (必ずしも物理的ではありませんが、心配する必要はありません) を持ち、おおよそ次のようなソリューションを期待しています。

重要体制図

これを実装するにはどうすればよいですか？

4

3 に答える 3

4

準形式的には、あなたが解決しようとしている問題は次のとおりです: 与えられた d0 に対して、論理式 "eq1(h, U0, d1c, d0) = eq2(h, U0, d1c, d0) = h と U0 の場合は 0"。

式を多項式「P(h, U0) = 0」に還元するアルゴリズムが存在します。これは「量指定子の除去」と呼ばれ、通常は別のアルゴリズム「円筒代数分解」に依存しています。残念ながら、これは sympy では (まだ) 実装されていません。

ただし、U0 は簡単に削除できるため、答えを見つけるために sympy でできることがあります。皮切りに

h, U0, d1c, d0 = symbols('h, U0, d1c, d0')
f1 = (U0) ** 2 * ((d0 ** 2 / d1c ** 3) + (1 - d0) ** 2 / (1 - d1c - d0 * h) ** 3) - 1
f2 = U0**2 / 2 * ((d0 ** 2 / d1c ** 2) + (1 - d0) ** 2 / (1 - d1c - d0 * h)) + d1c + d0 * (h - 1)

ここで、f1 から U0 を削除し、値を f2 に挿入します (より適切な式を得るために、solve() を使用するのではなく「手動で」実行しています)。

U2_val = ((f1 + 1)/U0**2)**-1
f3 = f2.subs(U0**2, U2_val)

f3 は h と d1c のみに依存します。また、これは有理分数であるため、その分子が 0 になるときだけを気にするので、2 つの変数で単一の多項式が得られます。

p3 = fraction(cancel(f3))

ここで、特定の d0 に対して、p3.subs(d0, .1) を数値的に反転して h(d1c) を取得し、それを U0 に戻し、(h, U0) のパラメトリック プロットを次の関数として作成できるはずです。 d1c。

于 2012-12-11T21:51:17.000 に答える
3

最初の除去を扱いましょうd1c。最初の方程式を の形にすることができたと想像してくださいd1c = f(U, h, d0)。次に、これを 2 番目の方程式に代入し、との間Uに特定の関係を持たせます。固定の場合、これは 2 つの変数との単一の方程式を定義し、原則として、そこから任意の を見つけることができます。最後のスケッチに基づいて、これはあなたがソリューションと呼んでいるもののようです。悪いニュースは、どちらの方程式からも簡単に得られないことです。良いニュースは、あなたがする必要がないということです。hd0d0UhUhd1c

fsolve連立方程式を取ることができます。たとえば、2 つの変数に依存し、解を与える 2 つの方程式です。この場合: fix h(d0は既に修正されています) であり、システムにフィードし、変数およびfsolveとして扱います。の値を記録し、 の次の値について繰り返します。U0d1cU0h

@duffymo の提案に反して、私は を使用することをお勧めしますfsolve。または、少なくともそれから始めて、蒸気がなくなった場合にのみ他のソルバーを探すことをお勧めします。

U0考えられる注意点は、与えられた に対して複数の解があることを期待していることですh:fsolve開始時の推測が必要であり、解の分岐の 1 つに収束するように指示する簡単な方法はありません。それが問題であることが判明した場合は、brentqソルバーを見てください。

別の方法は、システムから簡単に排除できることを確認することですU0。hこのようにして、との単一の方程式を取得し、の各値についてd1cそれを解き、元の方程式のいずれかを使用して与えられたとを計算します。d1chU0d1ch

使用例fsolve:

>>> from scipy.optimize import fsolve
>>> def f(x, p):
...   return x**2 -p
... 
>>> fsolve(f, 0.5, args=(2,))
array([ 1.41421356])
>>>

の場合に実際に解きたいことargs=(2,)を伝える構文を次に示します。は の値の開始推測です。fsolvef(x,2)=00.5x

于 2012-12-11T19:21:00.687 に答える
0

Newton Raphson や BFGS などの非線形ソルバーを使用して、連立非線形方程式を解きます。それらは、マトリックスの開始条件と調整に敏感であるため、注意が必要です。

于 2012-12-11T15:59:22.110 に答える