<jupyter_start><jupyter_text>Allee Effect $$\frac{d N}{dt} = r\, N\, \left( \frac{N}{A}-1\right)\, \left( 1-\frac{N}{K}\right)$$ with population size $N$, growth rate $r$, carrying capacity $K$, a critical. point $A$.The population has a negative growth rate for $0<N<A$ and a positive. growth rate for $A<N<K$ (assuming that $0<A<K$).<jupyter_code>def allee_growth(t, y, r=1, A=1, K=1):
    yp = r * y * (y/A - 1) * (1 - y/K)
    #yp = r * y *  (1 - y/K)
    return yp
# various plotting functions

import matplotlib.gridspec as gridspec

def find_nearest_value_index(array, value):
    array = np.asarray(array)
    idx = (np.abs(array - value)).argmin()
    return idx


def add_arrow(line, value=None, position=None, direction='right', size=15, color=None):
    """
    add an arrow to a line.

    line:       Line2D object
    position:   x-position of the arrow. If None, mean of xdata is taken
    direction:  'left' or 'right'
    size:       size of the arrow in fontsize points
    color:      if None, line color is taken.
    """
    if color is None:
        color = line.get_color()

    xdata = line.get_xdata()
    ydata = line.get_ydata()

    if value is not None:
        if value<0:
            direction = 'left'
        elif value>0:
            direction = 'right'
    if position is None:
        position = xdata.mean()
    # find closest index
    start_ind = np.argmin(np.absolute(xdata - position))
    if direction == 'right':
        end_ind = start_ind + 1
    else:
        end_ind = start_ind - 1

    line.axes.annotate('',
        xytext=(xdata[start_ind], ydata[start_ind]),
        xy=(xdata[end_ind], ydata[end_ind]),
        arrowprops=dict(arrowstyle="->", color=color, linewidth=2),
        size=size
    )

    
def solve_and_plot(y0, r, A, K=1, t_max=100):
    # x' vs x
    x_max = max([A, K, y0]) +0.1
    x_min = -0.1
    x = np.arange(x_min,x_max,0.01)
    
    xp = allee_growth(0, x, r=r, A=A, K=K)
    # y vs t
    t = np.arange(0,t_max,1)
    sol = solve_ivp(fun=lambda t, y: allee_growth(t, y, r=r, A=A, K=K), 
                    t_span=[0, t_max], y0=[y0], t_eval=t)

    # plot
    color_A = 'darkorange'
    color_K = 'green'
    color_y0= 'red'
    color_line = 'cornflowerblue'
    
    fig = plt.figure(figsize=(15, 10))
    gs = fig.add_gridspec(nrows=2, ncols=1, height_ratios=[1, 2])
    ax1 = fig.add_subplot(gs[0])
    ax2 = fig.add_subplot(gs[1])

    # x' vs x
    ax1.plot(x, xp, color=color_line)
    #ax1.axhline(y=0, color='k', linestyle=':')
    zero_line = ax1.plot(x, np.zeros(x.shape), color='k', linestyle=':')
    # arrows: 
    x_arr1 = x_min/2
    if A<K:
        x_arr2 = A/2
        x_arr3 = A+(K-A)/2
        x_arr4 = K+(x_max-K)/2
    else:
        x_arr2 = K/2
        x_arr3 = K+(A-K)/2
        x_arr4 = A+(x_max-A)/2
    for x_arr in [x_arr1, x_arr2, x_arr3, x_arr4]:
        val = xp[find_nearest_value_index(x, x_arr)]
        add_arrow(zero_line[0], position=x_arr, value=val, size=30, color='k')
    
    ax1.axvline(x=y0, color=color_y0, linestyle=':', label="$y_0=$%.2f"%y0)
    ax1.axvline(x=A, color=color_A, linestyle='--', label="$A=$%.2f"%A)
    ax1.axvline(x=K, color=color_K, linestyle='--', label="$K=$%.2f"%K)
    ax1.plot(0,0, color=color_line, marker='o')
    ax1.plot(A,0, color=color_line, marker='o')
    ax1.plot(K,0, color=color_line, marker='o')
    ax1.set_xlabel('$y$')
    ax1.set_ylabel("$y'$")
    #axes[0].legend(loc='lower right')
    ax1.grid(which='minor', alpha=0.5)
    ax1.grid(which='major', alpha=0.5)
    ax1.minorticks_on()
    ax1.set_xlim(x[0],x[-1])
    # y vs t
    ax2.plot(sol.t, sol.y[0],color=color_line)
    ax2.axhline(y=0, color='k', linestyle=':')
    ax2.axhline(y=y0, color=color_y0, linestyle=':', label="$y_0=$%.2f"%y0)
    ax2.axhline(y=A, color=color_A, linestyle='--', label="$A=$%.2f"%A)
    ax2.axhline(y=K, color=color_K, linestyle='--', label="$K=$%.2f"%K)
    ax2.plot(0, y0, color=color_y0, marker='o')
    ax2.set_xlabel("time")
    ax2.set_ylabel("population size")
    y_range_max = max([A, K])
    if y0<y_range_max:
        ax2.set_ylim(-0.05,y_range_max+0.05)
    # Put a legend to the right of the current axis
    ax2.legend(loc='center', bbox_to_anchor=(0.5, -0.35), ncol=3)
    ax2.grid(which='minor', alpha=0.5)
    ax2.grid(which='major', alpha=0.5)
    ax2.minorticks_on()
    plt.tight_layout(pad=0.4, w_pad=0.5, h_pad=1.0)
    plt.show()
    
#solve_and_plot(y0=1.3, r=0.1, A=0.5, K=1, t_max=200)
# interactive ipython widgets
import ipywidgets as widgets
from IPython.display import display

#-- initial values
w_y0 = widgets.FloatSlider(value=0.1, min=0.0, max=2.0, step=0.01, continuous_update=False, 
                                   description=r'\(y_0\)', orientation='horizontal', readout=True)
w_r  = widgets.FloatSlider(value=0.1, min=0.0, max=1.0, step=0.01, continuous_update=False, 
                           description=r'\(r\)', orientation='horizontal', readout=True)
w_A  = widgets.FloatSlider(value=0.2, min=0, max=2, step=0.01, continuous_update=False, 
                         description=r'\(A\)', orientation='horizontal', readout=True)
w_K  = widgets.FloatSlider(value=1, min=0, max=2, step=0.01, continuous_update=False, 
                         description=r'\(K\)', orientation='horizontal', readout=True)
w_t  = widgets.IntSlider(value=200, min=10, max=1000, step=1, continuous_update=False, 
                         description=r'\(t_{max}\)', orientation='horizontal', readout=True)

box1  = widgets.VBox( [w_y0, w_t] )
box2  = widgets.VBox( [w_r] )
box3 = widgets.VBox([w_A, w_K])
ui = widgets.HBox([box1, box2, box3])

out = widgets.interactive_output(solve_and_plot, { 
                                                'y0': w_y0,                                            
                                                'r' : w_r,
                                                'A' : w_A,
                                                'K' : w_K,
                                                't_max' : w_t})

# # display widget
display(ui, out)<jupyter_output><empty_output><jupyter_text>Bifurcation Diagram $$\frac{d N}{dt} = r\, N\, \left( \frac{N}{A}-1\right)\, \left( 1-\frac{N}{K}\right)$$**Fix Points** in function of $A$:- $N_1^*=0$- $N_{2, 3}^*(A)=\frac{1}{2}\left (A+K \pm \sqrt{(A+K)^2-4AK}\right)$, for $N\neq 0,\; r\neq 0$<jupyter_code># bifurcation diagram in function of A

def x_2(A, K=1):
    return  ((A+K)+np.sqrt(np.power(A+K,2)-4*A*K))/2.

def x_3(A, K=1):
    return  ((A+K)-np.sqrt(np.power(A+K,2)-4*A*K))/2.

A=np.arange(0,2,0.1)
K=1
x1 = np.zeros(A.shape)
x2 = x_2(A, K=K)
x3 = x_3(A, K=K)


color_A = 'darkorange'
color_K = 'green'

fig, ax = plt.subplots(1,1, figsize=(15, 10))

ax.axhline(y=K, linestyle=':', color=color_K)
ax.axvline(x=K, linestyle=':', color=color_K)

ax.plot(A, x1, color='darkmagenta')
ax.plot(A, x2, color='darkblue')
ax.plot(A, x3, '--', color='darkred')
ax.set_xlabel('A')
ax.set_ylabel('$N^*$')



plt.plot([0, K],[0, K], 'ko')<jupyter_output><empty_output>
{'fraction_non_alphanumeric': 0.18267312624679397, 'fraction_numerical': 0.036905101168424055, 'mean_word_length': 2.7059134107708553, 'pattern_counts': {'":': 0, '<': 18, '<?xml version=': 0, '>': 11, 'https://': 0, 'lorem ipsum': 0, 'www.': 0, 'xml': 0}, 'pii_count': 0, 'substrings_counts': 0, 'word_list_counts': {'cursed_substrings.json': 12, 'profanity_word_list.json': 0, 'sexual_word_list.json': 0, 'zh_pornsignals.json': 0}}
{'chain_length': '2', 'dir': 'jupyter-structured-clean-dedup', 'id': '155583', 'license': 'CC-BY-4.0', 'n_tokens_mistral': 2945, 'n_tokens_neox': 2853, 'n_words': 550, 'path': 'GrowthModels/AlleeEffect.ipynb', 'repo_name': 'sbranciamore/biosci670', 'size': '7018'}