ちょっと頑張ってL/D効果も含めた侵徹深さの評価をしたい話で
\[ \left(1+\frac{\rho_tu^2}{Y_t}\right)\sqrt{K_t-\rho_tu^2\alpha^2}= \left(1+\frac{\rho_tu^2\alpha^2}{2G_t}\right)\sqrt{K_t-\rho_tu^2} \]
に対してニュートンラフソン法、あるいは二分法を用いれば$\alpha$を求められると書いていました。これは、この式が平方根を含むために直接解析解を求めることができないためで、AWモデルを用いた計算の大部分を$\alpha$を求めるために消費していました。
一方で、$\alpha$の値域は$\alpha \gt 1$の実数というのは明らかですので、
\[ \left(1 + \frac{\rho u^{2}}{Y}\right)^{2} \left(K_{t} - \alpha^{2} \rho u^{2}\right) = \left(1 + \frac{\alpha^{2} \rho u^{2}}{2 G}\right)^{2} \left(K_{t} - \rho u^{2}\right) \]
とした解を求めても基本的には等価です。ここで$\alpha^2$については明らかに解析解が存在していて、sympyで適当に求めると、
…ここに示しておきたかったんですが長すぎて数式に表示されませんでした。
import sympy
rho,u2, Y, K_t, alpha, G = sympy.symbols("rho, u^2, Y, K_t, alpha^2, G")
lh = (1 + rho*u2 /Y) * (1 + rho*u2 / Y)* (K_t - rho*u2 * alpha)
rh = (1 + rho*u2 * alpha * 0.5 / G) * (1 + rho*u2 * alpha /2 / G) * (K_t - rho*u2)
eq = sympy.Eq(lh,rh)
s1,s2 = sympy.solve(eq,alpha)
s2
と書くことが出来ます、とだけ示しておきます。これの平方根を取れば、$\alpha$になります。数値解との対応は以下のとおりです。
個人的には計算時間が4msから1msまで減らせたので満足です。
