Набор тестов и проверок, что все правильно считается в неустойчивостях



In [1]:

    
from IPython.display import HTML
from IPython.display import Image
from PIL import Image as ImagePIL

%pylab
%matplotlib inline









    



Using matplotlib backend: Qt4Agg
Populating the interactive namespace from numpy and matplotlib



In [2]:

    
%run ../../utils/load_notebook.py # импортируем функцию загрузки ноутбука



In [3]:

    
from instabilities import * # загружаем функции









    



importing Jupyter notebook from instabilities.ipynb
Using matplotlib backend: Qt4Agg
Populating the interactive namespace from numpy and matplotlib



In [4]:

    
%run ../../utils/show_notebook.py



In [5]:

    
show_notebook("instabilities.ipynb") # показываем функции из ноутбука, чтобы можно было подглядеть в определение









    




markdown cell
# Неустойчивости

Тут будут все функции, имеющие отношение к неустойчивостям, чтобы не копировать их каждый раз.
Тесты на эти функции в соседнем ноутбуке.
code cell
from IPython.display import HTML
from IPython.display import Image
from PIL import Image as ImagePIL

%pylab
%matplotlib inline


markdown cell
# Одножидкостный критерий
Устойчиво, когда > 1:
$$Q_g = \frac{\Sigma_g^{cr}}{\Sigma_g}=\frac{\kappa c_g}{\pi G \Sigma_g}$$
$$Q_s = \frac{\Sigma_s^{cr}}{\Sigma_s}=\frac{\sigma_R}{\sigma_R^{min}}=\frac{\kappa \sigma_R}{3.36 G \Sigma_s}$$
code cell
G = 4.32 #гравитационная постоянная в нужных единицах

def Qs(epicycl=None, sigma=None, star_density=None):
    '''Вычисление безразмерного параметра Тумре для звездного диска. 
    Зависит от плотности звезд, дисперсии скоростей и эпициклической частоты.'''
    return epicycl * sigma / (3.36 * G * star_density)


def Qg(epicycl=None, sound_vel=None, gas_density=None):
    '''Вычисление безразмерного параметра Тумре для газового диска. 
    Зависит от плотности газа и эпициклической частоты, скорости звука в газе.'''
    return epicycl * sound_vel / (math.pi * G * gas_density)


markdown cell
# Двухжидкостный критерий

Кинетическое приближение:
$$\frac{1}{Q_{\mathrm{eff}}}=\frac{2}{Q_{\mathrm{s}}}\frac{1}{\bar{k}}\left[1-e^{-\bar{k}^{2}}I_{0}(\bar{k}^{2})\right]+\frac{2}{Q_{\mathrm{g}}}s\frac{\bar{k}}{1+\bar{k}^{2}s^{2}}>1\,$$

Гидродинамическое приближение:
$$\frac{2\,\pi\, G\, k\,\Sigma_{\mathrm{s}}}{\kappa+k^{2}\sigma_{\mathrm{s}}}+\frac{2\,\pi\, G\, k\,\Sigma_{\mathrm{g}}}{\kappa+k^{2}c_{\mathrm{g}}}>1$$ или $$\frac{1}{Q_{\mathrm{eff}}}=\frac{2}{Q_{\mathrm{s}}}\frac{\bar{k}}{1+\bar{k}^{2}}+\frac{2}{Q_{\mathrm{g}}}s\frac{\bar{k}}{1+\bar{k}^{2}s^{2}}>1$$ для безразмерного волнового числа ${\displaystyle \bar{k}\equiv\frac{k\,\sigma_{\mathrm{s}}}{\kappa}},\, s=c/\sigma$
code cell
from scipy.special import i0e, i1e

def inverse_hydro_Qeff_from_k(dimlK, Qg=None, Qs=None, s=None):
    return 2.*dimlK / Qs / (1 + dimlK**2) + 2*s*dimlK / Qg / (1 + dimlK**2 * s**2)

def inverse_kinem_Qeff_from_k(dimlK, Qg=None, Qs=None, s=None):
    return 2. / dimlK / Qs * (1 - i0e(dimlK ** 2)) + 2*s*dimlK / Qg / (1 + dimlK**2 * s**2)


markdown cell
### Нахождение максимума: 
markdown cell
Найти максимум функции в гидродинамическом приближении вообще просто - это многочлен
$$\frac{2}{Q_{\mathrm{s}}}\frac{\bar{k}}{1+\bar{k}^{2}}+\frac{2}{Q_{\mathrm{g}}}s\frac{\bar{k}}{1+\bar{k}^{2}s^{2}}>1$$
и у него можно найти максимум методами Sympy, взяв производную:
code cell
from sympy import *

def findInvHydroQeffSympy(Qs, Qg, s):
    '''Решаем уравнение deriv()=0 чтобы найти максимум функции в гидродинамическом приближении.'''
    k = Symbol('k') #solve for complex because it may returns roots as 1.03957287978471 + 0.e-20*I 
    foo = 2./Qs*k/(1+k**2) + 2/Qg*s*k/(1+k**2 * s**2)
    foo2 = 2./Qs * (1-k)*(1+k * s**2)**2 + 2/Qg*s*(1-k*s**2)*(1+k)**2
    roots = solve(foo2.simplify(), k)
    roots = [np.sqrt(float(abs(re(r)))) for r in roots]
    _tmp = [foo.evalf(subs={k:r}) for r in roots]
    max_val = max(_tmp)
    return (roots[_tmp.index(max_val)], max_val)


code cell
def findInvHydroQeffBrute(Qs, Qg, s, krange):
    '''Находим максимум функции в гидродинамическом приближении перебором по сетке.'''
    _tmp = [inverse_hydro_Qeff_from_k(l, Qg=Qg, Qs=Qs, s=s) for l in krange]
    max_val = max(_tmp)
    root_for_max = krange[_tmp.index(max_val)]
    if abs(root_for_max-krange[-1]) < 0.5:
        print 'WARNING! For Qs={} Qg={} s={} root of max near the max of k-range'.format(Qs, Qg, s)
    return (root_for_max, max_val)


code cell
from scipy.optimize import brentq

def findInvHydroQeffBrentq(Qs, Qg, s, krange):
    '''Решение уравнения deriv(9) = 0 для нахождения максимума исходной функции. Запускается brentq на исходной сетке,
    в случае если на концах сетки разные знаки функции (промежуток содержит корень),
    затем выбираются лучшие корни, после чего ищется, какой их них дает максимум. Возвращается только этот корень.'''
    grid = krange
    args = [Qs, Qg, s]
    signs = [derivTwoFluidHydroQeff(x, *args) for x in grid]
    signs = map(lambda x: x / abs(x), signs)
    roots = []
    for i in range(0, signs.__len__() - 1):
        if signs[i] * signs[i + 1] < 0:
            roots.append(brentq(lambda x: derivTwoFluidHydroQeff(x, *args), grid[i], grid[i + 1]))
    original = [inverse_hydro_Qeff_from_k(l, Qg=Qg, Qs=Qs, s=s) for l in roots]
    root_for_max = roots[original.index(max(original))]
    if abs(root_for_max-krange[-1]) < 0.5:
        print 'WARNING! For Qs={} Qg={} s={} root of max near the max of k-range'.format(Qs, Qg, s)
    return (root_for_max, max(original))

def derivTwoFluidHydroQeff(dimlK, Qs, Qg, s):
    '''Производная по \bar{k} от левой части (9) для того, чтобы найти максимум.'''
    part1 = (1 - dimlK ** 2) / (1 + dimlK ** 2) ** 2
    part3 = (1 - (dimlK * s) ** 2) / (1 + (dimlK * s) ** 2) ** 2
    return (2 * part1 / Qs) + (2 * s * part3 / Qg)


markdown cell
Теперь кинематическое приближение:
$$\frac{2}{Q_{\mathrm{s}}}\frac{1}{\bar{k}}\left[1-e^{-\bar{k}^{2}}I_{0}(\bar{k}^{2})\right]+\frac{2}{Q_{\mathrm{g}}}s\frac{\bar{k}}{1+\bar{k}^{2}s^{2}}>1\,$$
Тут сложнее, честно уже не решить. остается два способа - брутфорсом и brentq, производная известна.
code cell
def findInvKinemQeffBrute(Qs, Qg, s, krange):
    '''Находим максимум функции в кинематическом приближении перебором по сетке.'''
    _tmp = [inverse_kinem_Qeff_from_k(l, Qg=Qg, Qs=Qs, s=s) for l in krange]
    max_val = max(_tmp)
    root_for_max = krange[_tmp.index(max_val)]
    if abs(root_for_max-krange[-1]) < 0.5:
        print 'WARNING! For Qs={} Qg={} s={} root of max near the max of k-range'.format(Qs, Qg, s)
    return (root_for_max, max_val)


code cell
def findInvKinemQeffBrentq(Qs, Qg, s, krange):
    '''Решение уравнения deriv(13) = 0 для нахождения максимума исходной функции. Запускается brentq на исходной сетке,
    в случае если на концах сетки разные знаки функции (промежуток содержит корень),
    затем выбираются лучшие корни, после чего ищется, какой их них дает максимум. Возвращается только этот корень.'''
    grid = krange
    args = [Qs, Qg, s]
    signs = [derivTwoFluidKinemQeff(x, *args) for x in grid]
    signs = map(lambda x: x / abs(x), signs)
    roots = []
    for i in range(0, signs.__len__() - 1):
        if signs[i] * signs[i + 1] < 0:
            roots.append(brentq(lambda x: derivTwoFluidKinemQeff(x, *args), grid[i], grid[i + 1]))
    original = [inverse_kinem_Qeff_from_k(l, Qg=Qg, Qs=Qs, s=s) for l in roots]
    root_for_max = roots[original.index(max(original))]
    if abs(root_for_max-krange[-1]) < 0.5:
        print 'WARNING! For Qs={} Qg={} s={} root of max near the max of k-range'.format(Qs, Qg, s)
    return (root_for_max, max(original))


def derivTwoFluidKinemQeff(dimlK, Qs, Qg, s):
    '''Производная по \bar{k} от левой части (13) для того, чтобы найти максимум. Коррекция за ассимптотику производится
    с помощью встроенных функций бесселя, нормированных на exp.'''
    part1 = (1 - i0e(dimlK ** 2)) / (-dimlK ** 2)
    part2 = (2 * dimlK * i0e(dimlK ** 2) - 2 * dimlK * i1e(dimlK ** 2)) / dimlK
    part3 = (1 - (dimlK * s) ** 2) / (1 + (dimlK * s) ** 2) ** 2
    return 2 * (part1 + part2) / Qs + 2 * s * part3 / Qg


code cell

Для случая бесконечного тонкого диска: $$\kappa=\frac{3}{R}\frac{d\Phi}{dR}+\frac{d^2\Phi}{dR^2}$$ где $\Phi$ - гравпотенциал, однако его знать не надо, т.к. есть проще формула: $$\kappa=\sqrt{2}\frac{\vartheta_c}{R}\sqrt{1+\frac{R}{\vartheta_c}\frac{d\vartheta_c}{dR}}$$



In [6]:

    
def model_kappa(R):
    return 100.*np.exp(-R/90.)

test_points = np.arange(20., 250., 10.)

fig = plt.figure(figsize=[8, 6])
plt.plot(test_points, model_kappa(test_points), '-')
plt.xlabel('$R, arcsec$')
plt.ylabel('$\kappa,\, km/s/kpc$', fontsize=15)
plt.show()



In [7]:

    
def model_sigmaR(R):
    return 200.*np.exp(-R/50.)

fig = plt.figure(figsize=[8, 6])
plt.plot(test_points, model_sigmaR(test_points), '-')
plt.xlabel('$R, arcsec$')
plt.ylabel('$\sigma_R,\, km/s$', fontsize=15)
plt.show()



In [8]:

    
def model_star(R):
    return 500.*np.exp(-R/50.)

def model_gas(R):
    return 20.*np.exp(-R/50.) + 15.*np.exp(-(R-50.)**2/50.)

fig = plt.figure(figsize=[8, 6])
plt.plot(test_points, model_gas(test_points), '-')
plt.plot(test_points, model_star(test_points), '-')
plt.xlabel('$R, arcsec$')
plt.ylabel('$\Sigma,\, M_{o}/pc^2$', fontsize=15)
plt.show()

Одножидкостный

Устойчиво, когда > 1: $$Q_g = \frac{\Sigma_g^{cr}}{\Sigma_g}=\frac{\kappa c_g}{\pi G \Sigma_g}$$ $$Q_s = \frac{\Sigma_s^{cr}}{\Sigma_s}=\frac{\sigma_R}{\sigma_R^{min}}=\frac{\kappa \sigma_R}{3.36 G \Sigma_s}$$



In [9]:

    
sound_vel = 6  #скорость звука в газе, км/с



In [10]:

    
fig = plt.figure(figsize=[8, 6])
plt.plot(test_points, [Qs(epicycl=epicycl, sigma=sigma, star_density=star_density) for (epicycl, sigma, star_density) in zip(
        model_kappa(test_points), model_sigmaR(test_points), model_star(test_points))], label='Qs')
plt.plot(test_points, [Qg(epicycl=epicycl, sound_vel=sound_vel_, gas_density=gas_density) for 
                       (epicycl, sound_vel_, gas_density) in zip(model_kappa(test_points), 
                                                                 [sound_vel]*len(test_points), 
                                                                 model_gas(test_points))], label='Qg')
plt.legend()
plt.show()

Двухжидкостная:

Кинетическое приближение: $$\frac{1}{Q_{\mathrm{eff}}}=\frac{2}{Q_{\mathrm{s}}}\frac{1}{\bar{k}}\left[1-e^{-\bar{k}^{2}}I_{0}(\bar{k}^{2})\right]+\frac{2}{Q_{\mathrm{g}}}s\frac{\bar{k}}{1+\bar{k}^{2}s^{2}}>1\,$$

Гидродинамическое приближение: $$\frac{2\,\pi\, G\, k\,\Sigma_{\mathrm{s}}}{\kappa+k^{2}\sigma_{\mathrm{s}}}+\frac{2\,\pi\, G\, k\,\Sigma_{\mathrm{g}}}{\kappa+k^{2}c_{\mathrm{g}}}>1$$ или $$\frac{1}{Q_{\mathrm{eff}}}=\frac{2}{Q_{\mathrm{s}}}\frac{\bar{k}}{1+\bar{k}^{2}}+\frac{2}{Q_{\mathrm{g}}}s\frac{\bar{k}}{1+\bar{k}^{2}s^{2}}>1$$ для безразмерного волнового числа ${\displaystyle \bar{k}\equiv\frac{k\,\sigma_{\mathrm{s}}}{\kappa}},\, s=c/\sigma$

✓Test:

проверим, что на каком-то участке значения inverse_kinem_Qeff_from_k и inverse_hydro_Qeff_from_k совпадают



In [11]:

    
fig = plt.figure(figsize=[16, 6])
pp_ = np.linspace(0.1, 100., 1000)
plt.plot(pp_, [inverse_hydro_Qeff_from_k(l, Qg=1., Qs=1., s=0.2) for l in pp_], '-')
plt.plot(pp_, [inverse_kinem_Qeff_from_k(l, Qg=1., Qs=1., s=0.2) for l in pp_], '-')
plt.ylabel(r'$\frac{1}{Q_{eff}}$', fontsize=25)
plt.show()

Видно, что различия есть (как и показывал Rafikov), но с хорошей точностью значения совпадают.

Зависимость от величины параметров:



In [12]:

    
inv_qss = np.arange(0.01, 1.2, 0.1)
inv_qgs = np.arange(0.01, 1.2, 0.1)

inv_qeffs = []

for qg in inv_qgs:
    for qs in inv_qss:
        inv_qeff = max([inverse_kinem_Qeff_from_k(l, Qg=1./qg, Qs=1./qs, s=0.2) for l in np.arange(0.1, 10000., 1.0)])
        inv_qeffs.append(inv_qeff)



In [13]:

    
fig = plt.figure(figsize=[8, 8])

ax1 = fig.add_subplot(111)
qeffs = np.array(inv_qeffs).reshape(len(inv_qgs), len(inv_qgs))
levels = np.linspace(qeffs.min(), qeffs.max(), 20)

cset=ax1.contour(qeffs, levels, hold='on', colors = 'k', origin='lower', extent=[inv_qgs[0],inv_qgs[-1],inv_qss[0],inv_qss[-1]])
ax1.clabel(cset, inline=1, fontsize=10, fmt='%1.2f')
ax1.grid()
ax1.set_xlabel('1/Qg')
ax1.set_ylabel('1/Qs')
ax1.set_title('1/Qeff, s = 0.2')
ax1.set_xticks(np.arange(0, 1.1, 0.1))
ax1.set_yticks(np.arange(0, 1.1, 0.1))

levels = np.array([1.0])
cset=ax1.contour(qeffs, levels, hold='on', colors = 'r', origin='lower', extent=[inv_qgs[0],inv_qgs[-1],inv_qss[0],inv_qss[-1]])

plt.show()

✓Test:

воспроизвести следующую иллюстрацию из Rafikov 2001 (Fig.3)



In [14]:

    
Image('rafikov_fig3.png')









    Out[14]:



In [15]:

    
inv_qss = np.arange(0.01, 1.2, 0.1)
inv_qgs = np.arange(0.01, 1.2, 0.1)

s_params = [1.0, 0.5, 0.3, 0.2, 0.1]
inv_qeffs_k = [[],[],[],[],[]]
inv_qeffs_h = [[],[],[],[],[]]

for ind, s in enumerate(s_params):
    for qg in inv_qgs:
        for qs in inv_qss:
            inv_qeff = max([inverse_kinem_Qeff_from_k(l, Qg=1./qg, Qs=1./qs, s=s) for l in np.arange(0.1, 10000., 1.0)])
            inv_qeffs_k[ind].append(inv_qeff)
            inv_qeff = max([inverse_hydro_Qeff_from_k(l, Qg=1./qg, Qs=1./qs, s=s) for l in np.arange(0.1, 10000., 1.0)])
            inv_qeffs_h[ind].append(inv_qeff)



In [16]:

    
fig = plt.figure(figsize=[8, 8])

ax1 = fig.add_subplot(111)
levels = np.array([1.])

for ind, s in enumerate(s_params):
    qeffs = np.array(inv_qeffs_h[ind]).reshape(len(inv_qgs), len(inv_qgs))
    cset=ax1.contour(qeffs, levels, hold='on', colors = 'k', origin='lower', 
                     extent=[inv_qgs[0],inv_qgs[-1],inv_qss[0],inv_qss[-1]], linestyles='dashed')
    qeffs = np.array(inv_qeffs_k[ind]).reshape(len(inv_qgs), len(inv_qgs))
    cset=ax1.contour(qeffs, levels, hold='on', colors = 'k', origin='lower', 
                     extent=[inv_qgs[0],inv_qgs[-1],inv_qss[0],inv_qss[-1]])
ax1.grid()
ax1.set_xlabel('1/Qg')
ax1.set_ylabel('1/Qs')
ax1.set_title('1/Qeff, s = 0.2')
ax1.set_xticks(np.arange(0, 1.1, 0.1))
ax1.set_yticks(np.arange(0, 1.1, 0.1))

plt.show()

Успешно!

Нахождение максимума функции:

В гидродинамическом приближении вообще просто - это многочлен $$\frac{2}{Q_{\mathrm{s}}}\frac{\bar{k}}{1+\bar{k}^{2}}+\frac{2}{Q_{\mathrm{g}}}s\frac{\bar{k}}{1+\bar{k}^{2}s^{2}}>1$$ и у него можно найти максимум методами Sympy, взяв производную:



In [17]:

    
from sympy import *
x = Symbol('x', real=True)
a = 1.5
f = 2./a*x/(1+x**2) + 2/0.9*0.4*x/(1+x**2 * 0.4**2)
roots = solve(f.diff(), x)
print roots
print max([f.evalf(subs={x:r}) for r in roots])









    



[-1.85855358154608, 1.85855358154608]
1.62034285003827



In [18]:

    
from sympy.plotting import plot
plot(f, (x, 0., 10.))









    












    Out[18]:





<sympy.plotting.plot.Plot at 0x9e4c358>



In [19]:

    
findInvHydroQeffSympy(1.32, 1.22, 0.52)









    Out[19]:





(1.408446000988677, 1.49668341702704)

Второй метод - тупо перебор по сетке:



In [20]:

    
findInvHydroQeffBrute(1.32, 1.22, 0.52, np.arange(0.01, 300., 0.01))









    Out[20]:





(1.4100000000000001, 1.4966826880822235)

Последний метод - используя производную, запускать brentq на нужном промежутке:



In [21]:

    
findInvHydroQeffBrentq(1.32, 1.22, 0.52, np.arange(0.01, 300., 0.01))









    Out[21]:





(1.4084460009886812, 1.496683417027044)

✓Test:

все три метода должны работать и давать почти одинаковые результаты на случайных данных.



In [22]:

    
cases_n = 1000
test_cases = zip(np.random.uniform(0.01, 3., cases_n), 
                 np.random.uniform(0.01, 3., cases_n), 
                 np.abs(np.random.normal(0.01, 0.5, cases_n)))



In [23]:

    
test_cases[0:10]









    Out[23]:





[(0.46590634709598588, 2.195401062730074, 0.4577795604792968),
 (1.6856797233429963, 2.9857988221438978, 0.03361509402127872),
 (2.4378383945209969, 1.9381950746086971, 0.71032872826613636),
 (1.439191818170283, 1.570165255204824, 0.81454629915080223),
 (2.0097893940272669, 0.46676098323732246, 0.54957897491578811),
 (2.7021777615670568, 2.0005488702185859, 0.048587842247888501),
 (0.55202016032025203, 2.2711912705100601, 0.17031920753056703),
 (0.94317892462366748, 1.8827633727852395, 0.25490235087831092),
 (2.3081298858531469, 1.9526260405642728, 0.46565157296700804),
 (1.4572992434684517, 0.63900453781634048, 0.91694709294906118)]



In [24]:

    
%time bren = [findInvHydroQeffBrentq(case[0], case[1], case[2], np.arange(0.01, 60000., 1.)) for case in test_cases]









    



Wall time: 4min 16s



In [25]:

    
%time sym = [findInvHydroQeffSympy(case[0], case[1], case[2]) for case in test_cases]









    



Wall time: 7min 22s



In [26]:

    
%time bru = [findInvHydroQeffBrute(case[0], case[1], case[2], np.arange(0.01, 3000., 0.01)) for case in test_cases]









    



WARNING! For Qs=0.842660154861 Qg=0.724569189699 s=0.000234916421128 root of max near the max of k-range
Wall time: 12min 26s



In [27]:

    
%time bru = [findInvHydroQeffBrute(c[0], c[1], c[2], np.arange(b[0]-100., b[0]+100., 0.001)) for c, b in zip(test_cases, bren)]









    



Wall time: 8min 2s



In [28]:

    
fig = plt.figure(figsize=[16, 8])
plt.plot(range(cases_n), abs(np.array(zip(*bren)[1])-np.array(zip(*bru)[1])), 's-')
plt.plot(range(cases_n), abs(np.array(zip(*bren)[1])-np.array(zip(*sym)[1])), 'd-')
plt.plot(range(cases_n), abs(np.array(zip(*sym)[1])-np.array(zip(*bru)[1])), 'v-')
plt.show()



In [29]:

    
fig, axes = plt.subplots(2, 3, figsize=[16, 10])
axes[0][0].plot(zip(*sym)[1], zip(*bru)[1], 'v')
axes[0][1].plot(zip(*sym)[1], zip(*bren)[1], '>')
axes[0][2].plot(zip(*bren)[1], zip(*bru)[1], '<')

axes[1][0].plot(zip(*sym)[0], zip(*bru)[0], '8')
axes[1][1].plot(zip(*sym)[0], zip(*bren)[0], 'p')
axes[1][2].plot(zip(*bren)[0], zip(*bru)[0], 'h')

for ax in axes:
    for a in ax:
        a.set_xlim(0., 10.)
        a.set_ylim(0., 10.)
        a.grid()
        a.plot([0,10], [0, 10], '--')
plt.show()

Внешнее согласие отличное.

Теперь кинематическое приближение: $$\frac{2}{Q_{\mathrm{s}}}\frac{1}{\bar{k}}\left[1-e^{-\bar{k}^{2}}I_{0}(\bar{k}^{2})\right]+\frac{2}{Q_{\mathrm{g}}}s\frac{\bar{k}}{1+\bar{k}^{2}s^{2}}>1\,$$ Тут сложнее, честно уже не решить. остается два способа - брутфорсом и brentq, производная известна.



In [30]:

    
findInvKinemQeffBrute(1.32, 1.22, 0.52, np.arange(0.01, 300., 0.01))









    Out[30]:





(1.3, 1.5278157625496003)



In [31]:

    
findInvKinemQeffBrentq(1.32, 1.22, 0.52, np.arange(0.01, 300., 0.01))









    Out[31]:





(1.3012816962470588, 1.5278164121945141)

✓Test:

оба метода должны давать схожие величины и работать в принципе



In [32]:

    
cases_n = 1000
test_cases = zip(np.random.uniform(0.01, 3., cases_n), 
                 np.random.uniform(0.01, 3., cases_n), 
                 np.abs(np.random.normal(0.01, 0.5, cases_n)))



In [33]:

    
test_cases[0:10]









    Out[33]:





[(2.9816301250366029, 0.49583480891652859, 0.50679046379497839),
 (2.4546699444881042, 2.2081201057179007, 0.034452973748545251),
 (2.5373262448757834, 0.56041404738584999, 0.21219789505497341),
 (2.7595149207194245, 0.54472571428973338, 0.37995760059045336),
 (2.3631763592386861, 2.0123314345868781, 0.24750093410065591),
 (2.9186558079909353, 2.3848753894318402, 0.15991386706971089),
 (1.4969430951231046, 2.5105768114111662, 0.49298576757316248),
 (0.29186782067183092, 1.9526461798016765, 0.71353884108189758),
 (0.14696390709855933, 2.0559147869484455, 0.89347274746895133),
 (2.8028372130470709, 0.24344467998663694, 0.23645961488818681)]



In [34]:

    
%time bren = [findInvKinemQeffBrentq(case[0], case[1], case[2], np.arange(0.01, 60000., 1.)) for case in test_cases]









    



Wall time: 8min 34s



In [35]:

    
%time bru = [findInvKinemQeffBrute(case[0], case[1], case[2], np.arange(0.01, 3000., 0.01)) for case in test_cases]









    



Wall time: 18min 49s



In [36]:

    
%time bru = [findInvKinemQeffBrute(c[0], c[1], c[2], np.arange(b[0]-100., b[0]+100., 0.001)) for c, b in zip(test_cases, bren)]









    



Wall time: 12min 40s



In [37]:

    
fig = plt.figure(figsize=[16, 5])
plt.plot(range(cases_n), abs(np.array(zip(*bren)[1])-np.array(zip(*bru)[1])), 's-')
plt.show()



In [38]:

    
fig, axes = plt.subplots(1, 2, figsize=[12, 6])
axes[0].plot(zip(*bren)[1], zip(*bru)[1], '<')
axes[1].plot(zip(*bren)[0], zip(*bru)[0], 'h')

for ax in axes:
    ax.set_xlim(0., 10.)
    ax.set_ylim(0., 10.)
    ax.grid()
    ax.plot([0,10], [0, 10], '--')
plt.show()

Работает

✓Test:

Кинематические решения должны быть не очень далеко от гидродинамических - Kinem ~ Hydro



In [39]:

    
cases_n = 200
test_cases = zip(np.random.uniform(0.01, 3., cases_n), 
                 np.random.uniform(0.01, 3., cases_n), 
                 np.abs(np.random.normal(0.01, 0.5, cases_n)))



In [40]:

    
test_cases[0:10]









    Out[40]:





[(2.9974977463186541, 0.67174053464992034, 0.00028979580141474254),
 (0.27586990963569941, 1.9891871268476615, 0.17832811773198312),
 (1.2431541173938563, 1.2658035904694929, 0.56687704808003869),
 (0.12360972017731459, 2.9365423924979299, 0.30774064769198683),
 (1.6237534243173859, 1.0358788334940487, 0.21076809551684791),
 (1.2612691708436041, 0.46286099591573066, 0.49896674494201582),
 (0.47145867798130348, 2.8478533353855702, 0.34794639546405076),
 (2.7616221056932218, 2.0884660215131254, 0.50281815084474069),
 (2.4560504424927392, 0.42193433374259398, 0.2913354237251522),
 (1.8025734919806378, 1.1782129525195597, 0.21234063995262512)]



In [41]:

    
%time kin = [findInvKinemQeffBrentq(case[0], case[1], case[2], np.arange(0.01, 60000., 1.)) for case in test_cases]









    



Wall time: 1min 46s



In [42]:

    
%time hyd = [findInvHydroQeffBrentq(case[0], case[1], case[2], np.arange(0.01, 60000., 1.)) for case in test_cases]









    



Wall time: 53.9 s



In [43]:

    
fig = plt.figure(figsize=[16, 5])
plt.plot(range(cases_n), abs(np.array(zip(*kin)[1])-np.array(zip(*hyd)[1])), 's-', color='red')
plt.plot(range(cases_n), abs(np.array(zip(*kin)[0])-np.array(zip(*hyd)[0])), 'v-')
plt.ylim(0, 5.)
plt.show()



In [44]:

    
fig, axes = plt.subplots(1, 2, figsize=[16, 7])
axes[0].plot(zip(*kin)[1], zip(*hyd)[1], '<')
axes[1].plot(zip(*kin)[0], zip(*hyd)[0], 'h')

for ax in axes:
    ax.set_xlim(0., 10.)
    ax.set_ylim(0., 10.)
    ax.grid()
    ax.plot([0,10], [0, 10], '--')
plt.show()

Видно, что максимум иногда съезжает и Q иногда различается, но в целом все достаточно близко.

✓Test:

$$Q_{2F} \le Q_{1F}$$



In [45]:

    
# используем сосчитанные выше кейсы
%time qgs = [1./case[1] for case in test_cases]









    



Wall time: 0 ns



In [46]:

    
fig, axes = plt.subplots(1, 2, figsize=[16, 7])
axes[0].plot(qgs, zip(*hyd)[1], '<')
axes[1].plot(qgs, zip(*kin)[1], 'h')

for ax in axes:
    ax.set_xlim(0., 10.)
    ax.set_ylim(0., 10.)
    ax.grid()
    ax.plot([0,10], [0, 10], '--')
    ax.set_xlabel('1/Qg')
    ax.set_ylabel('1/Qeff')
plt.show()

Верно

✓Test:

воспроизвести график из следующей иллюстрации Rafikov 2001



In [47]:

    
Image('rafikov_fig1.png')









    Out[47]: