Python(sympy)で、2変数を伴った1本の式（陰関数）から両者の関係を表す数式を導きたい

Question

前提・実現したいこと

2つの変数を伴った陰関数（1本の式）において、片方の変数をもう片方の変数の関数として表したい（関係を表す式を導きたい）
と考えております。

Python 初心者のため、ヒント・解決のための方法等を教えていただけると有り難いです。

よろしくお願いいたします。

発生している問題・エラーメッセージ

以下の数式（JCxbcul）において、θ を x_b の関数として導きたいのですが、実行結果が得られません。
与えた数式に対して sympy.solve を実行すると、そのまま動かなくなってしまいます。
単に式が複雑だから時間がかかっているのか、それとも根本的に対応が間違っているのか、また他に効果的な方法があるのか、悩んでいます。

なお以下の JCxbcul = で表されている式は、私が現在取り組んでいる経済学の理論モデルから導いたものです。

該当のソースコード

import sympy as sm
import numpy as np
import matplotlib.pyplot as plt
sm.init_printing()
sm.var('x_b, θ')

JCxbcul = 12.4144736842105*θ**0.5 - 0.5*((-0.0192084388185654*θ**0.5/(0.04598*θ**0.5 + 0.043659) + 0.252742616033755*(0.076*θ**0.5 + 0.1468)/((x_b + 0.605)*(0.076*θ**0.5 + 0.1078) + 0.0078) - 1.26371308016878*sm.log(13.1578947368421*θ**(-0.5)*(0.04598*θ**0.5 + 0.043659)) + 1.26371308016878*sm.log(((x_b + 0.605)*(0.076*θ**0.5 + 0.1078) + 0.0078)/(0.076*θ**0.5 + 0.1468)))/(0.076*θ**0.5 + 0.1078) + (-0.0192084388185654*θ**0.5/((x_b + 0.605)*(0.076*θ**0.5 + 0.0078) - 0.00156) - 1.26371308016878*sm.log(13.1578947368421*θ**(-0.5)*((x_b + 0.605)*(0.076*θ**0.5 + 0.0078) - 0.00156)) + 2.03291512368581)/(0.076*θ**0.5 + 0.0078) + (0.020042194092827*θ**0.5*(0.4792*x_b + 0.33476912)/((x_b + 0.605)*(0.076*θ**0.5 + 0.0078) - 0.00156) - 0.263713080168776*(0.4792*x_b + 0.33476912)*(0.076*θ**0.5 + 0.1468)/((x_b + 0.605)*(0.076*θ**0.5 + 0.1078) + 0.0078) + 0.185513080168776*sm.log(13.1578947368421*θ**(-0.5)*((x_b + 0.605)*(0.076*θ**0.5 + 0.0078) - 0.00156)) - 0.185513080168776*sm.log(((x_b + 0.605)*(0.076*θ**0.5 + 0.1078) + 0.0078)/(0.076*θ**0.5 + 0.1468)))/(0.0035568*θ**0.5 + 0.00114504))/(0.096042194092827*θ**0.5/(0.04598*θ**0.5 + 0.043659) - 0.263713080168776)

sm.solve(JCxbcul, θ)

####################################################

#図による確認 
xrange = np.linspace(0.1, 2, 20)
yrange = np.linspace(0.1, 2, 20)
x_b, θ = np.meshgrid(xrange,yrange)
plt.figure(figsize=(5,5))

plt.axis([0.1, 2, 0.1, 2])

Z=12.4144736842105*θ**0.5 - 0.5*((-0.0192084388185654*θ**0.5/(0.04598*θ**0.5 + 0.043659) + 0.252742616033755*(0.076*θ**0.5 + 0.1468)/((x_b + 0.605)*(0.076*θ**0.5 + 0.1078) + 0.0078) - 1.26371308016878*np.log(13.1578947368421*θ**(-0.5)*(0.04598*θ**0.5 + 0.043659)) + 1.26371308016878*np.log(((x_b + 0.605)*(0.076*θ**0.5 + 0.1078) + 0.0078)/(0.076*θ**0.5 + 0.1468)))/(0.076*θ**0.5 + 0.1078) + (-0.0192084388185654*θ**0.5/((x_b + 0.605)*(0.076*θ**0.5 + 0.0078) - 0.00156) - 1.26371308016878*np.log(13.1578947368421*θ**(-0.5)*((x_b + 0.605)*(0.076*θ**0.5 + 0.0078) - 0.00156)) + 2.03291512368581)/(0.076*θ**0.5 + 0.0078) + (0.020042194092827*θ**0.5*(0.4792*x_b + 0.33476912)/((x_b + 0.605)*(0.076*θ**0.5 + 0.0078) - 0.00156) - 0.263713080168776*(0.4792*x_b + 0.33476912)*(0.076*θ**0.5 + 0.1468)/((x_b + 0.605)*(0.076*θ**0.5 + 0.1078) + 0.0078) + 0.185513080168776*np.log(13.1578947368421*θ**(-0.5)*((x_b + 0.605)*(0.076*θ**0.5 + 0.0078) - 0.00156)) - 0.185513080168776*np.log(((x_b + 0.605)*(0.076*θ**0.5 + 0.1078) + 0.0078)/(0.076*θ**0.5 + 0.1468)))/(0.0035568*θ**0.5 + 0.00114504))/(0.096042194092827*θ**0.5/(0.04598*θ**0.5 + 0.043659) - 0.263713080168776)
plt.grid(True)
#Z = 10*θ**2 + μ_g - 1
plt.contour(x_b,θ, Z, [0])
plt.show()

試したこと

matplotlib で図を作成した所、右下がりの曲線が描かれました（そちらのコードも上に記しておきました）。

補足情報（FW/ツールのバージョンなど）

Jupyter Notebook (Python 3) や Spyder(Python 3.7)を使用しています。

Answer 1

私の環境でも、1時間以上経っても終わりませんね。
θが全て0.5乗されていたので別変数に置き換えてみたりもしましたが、それでも終わりません。
そのまま解析的に解くには式が複雑すぎるようです。

・プログラム走らせる前に、数学的処理でもう少し簡易な式に変形する。
・速いPCを使う。
・ひたすら放置して回す。
というあたりになるかと思います。後ろ2つはオススメできませんが…

それと、プロットはすぐ描けるので、数値解は簡単に出るようです。
離散的に解くとか、近似解を出すみたいなやり方は（もとの研究方針が許すなら）ありかもしれません。

Comments

[marimoswordsman]

ぺん様：
ご回答くださり、有り難うございます。
また返信が送れてしまい申し訳ございません。

他の方の環境でも解析的に解くことが難しいことが分かったのは大きな収穫です。
近似解を出す方法を模索してみたり、解析的に解かなくて済むような分析の方向性がないか検討してみようと思います。

お忙しい中ご対応いただき、本当に有り難うございました。