<jupyter_start><jupyter_text>Let's generate a random 2x2 matrix and look at the dynamics of the resulting linear system$$\begin{bmatrix}\dot{x}_1 \\ \dot{x}_2\end{bmatrix}=\begin{bmatrix}a & b \\c & d \end{bmatrix}\begin{bmatrix}x_1 \\ x_2\end{bmatrix}$$$$\dot{x} = Ax$$<jupyter_code>A = 20*rand(2, 2) - 10
print('A=')
print(A)
L, V = eig(A) ## the function `eig` computes numerical approximation for eigenvalues and eigenvectors
print('eigenvalues=')
print(L)
V = matrix(V) # the eigenvectors v1 and v2 are the columns of the matrix V
def xsol(t, x0): # the solution to the linear system, given a time t and an initial condition x0
    c = V.I*x0 # the constants are determined by the initial condition
    return V*diag(exp(L*t))*c # this is equivalent to c1*v1*exp(lambda1*t) + c2*v2*exp(lambda2*t)<jupyter_output>A=
[[ 4.00414813  9.10666051]
 [-7.93856479  2.30465515]]
eigenvalues=
[3.15440164+8.46000859j 3.15440164-8.46000859j]<jupyter_text>We will pick 10 random inital conditions and plot the trajectories<jupyter_code>t = linspace(0, 2, 200)
figure(1, [10, 6])
for n in arange(10):
    x0 = matrix(2*rand(2) -1).reshape(2, 1) # pick a random initial condition
    x = array([xsol(tn, x0) for tn in t]) # make the solution
    plot(x0[0], x0[1], 'ok') # plot the inital condition
    plot(x[:, 0], x[:, 1]) # plot the trajectory
xlim(-2, 2)
ylim(-2, 2);
xlabel('$x_1$', fontsize=24)
ylabel('$x_2$', fontsize=24);<jupyter_output>/anaconda3/lib/python3.6/site-packages/numpy/core/_asarray.py:85: ComplexWarning: Casting complex values to real discards the imaginary part
  return array(a, dtype, copy=False, order=order)<jupyter_text>Fishing model 3: Allee effectAn 'allee' effect occurs when the model has three fixed points: a stable fixed point at $x=0$, an unstable fixed point at $\hat{x}_1>0$ and a stable carrying capacaty at $\hat{x}_2 > \hat{x}_1$. When the population drops below the unstable fixed point $\hat{x}_1$, the population converges to the zero fixed point which corresponds to extinction. This models the effect of a minimum viable population necessary to thrive with a stable long term population level.For certain parameter values, the fishing model can display this behavior. Here, I've modified the birth rate to be $\epsilon + rN$ so that extinction is impossible. What we see is an example of stochastic switching between the stable fixed points of the ODE model.<jupyter_code>r = 2.
K = 10.
H = 2.5
A = 2.

x = linspace(-0.5, 12, 200)

y0 = r*x*(1 - x/K)
y = r*x*(1 - x/K) - H*x/(A + x)
y1 = r*x*(1 - x/K) - 5.5*x/(A + x)

figure(1, [10, 5])
plot(x, y)
plot(x, y0)
plot(x, y1)
plot(x, 0*x, 'k')

ylim(-0.5, None)
xlim(-0.1, 12,5)

xlabel('x', fontsize=20)
ylabel('f(x)', fontsize=20);
epsilon = 0.2 # TRY CHANGING THIS VALUE
r = 2.
K = 10.
H = 5.5
A = 2.
# K = 1.
n0 = int(5/epsilon)
Nsteps = 50000
N = zeros(Nsteps)
N[0] = n0
T = zeros(Nsteps)
T[0] = 0
for j in arange(1, Nsteps):
    n = N[j-1]
    combined_rate = epsilon + n*r + n**2*epsilon*r/K + H*n/(A + epsilon*n)
    u = rand(1)
    tau = -log(u)/combined_rate
    T[j] = T[j-1] + tau # update time array with new time value
    ## now we make a random variable to decide if birth or death occurs
    probability_of_birth = (epsilon + n*r)/combined_rate
    probability_of_death = 1 - probability_of_birth
    u2 = rand(1)
    if u2 < probability_of_birth:
        N[j] = N[j-1] + 1 # a birth means N -> N+1
    else:
        N[j] = N[j-1] - 1
        
#     if N[j] == 0:
#         T = T[:j+1]
#         N = N[:j+1]
#         break
        
def ode(x0, Tmax): # change `Nsteps` to alter
    ### Euler's method parameters
    Delta_t = 0.1 # The size of the time step
    Nsteps = int(Tmax/Delta_t)
    x = zeros(Nsteps) # an array of zeros to store the solution in
    x[0] = x0 # store the initial condition
    def f_ode(x): # the ligistic growth RHS function
        return r*x*(1-x/K) - H*x/(A + x)
    ####### The simulation loop #########
    for n in arange(1, Nsteps): # loop from n=1 up to n=Nsteps-1
        x[n] = x[n-1] + Delta_t*f_ode(x[n-1])
    t = Delta_t*arange(Nsteps)
    return t, x 
    
Tode1, Xode1 = ode(.5, T[-1])
Tode2, Xode2 = ode(5, T[-1])
######## plot the result
figure(1)
plot(T, N, '-o') # the option '-o' makes a line with 'o' symbols at the jump times
plot(Tode1, Xode1/epsilon, 'k')
plot(Tode2, Xode2/epsilon, 'k')
xlabel('t', fontsize=24)
ylabel('N(t)', fontsize=24);

figure(2)
P, b = histogram(N[1:], bins=50, weights=diff(T));
# semilogy(b[:-1], P)
# bar(b[:-1], P)
hist(N[1:],  bins=arange(70), weights=diff(T));<jupyter_output><empty_output>
{'fraction_non_alphanumeric': 0.1564001747487986, 'fraction_numerical': 0.055701179554390565, 'mean_word_length': 3.6581892166836214, 'pattern_counts': {'":': 0, '<': 12, '<?xml version=': 0, '>': 14, 'https://': 0, 'lorem ipsum': 0, 'www.': 0, 'xml': 0}, 'pii_count': 0, 'substrings_counts': 0, 'word_list_counts': {'cursed_substrings.json': 0, 'profanity_word_list.json': 0, 'sexual_word_list.json': 0, 'zh_pornsignals.json': 0}}
{'chain_length': '3', 'dir': 'jupyter-structured-clean-dedup', 'id': '104067', 'license': 'Apache-2.0', 'n_tokens_mistral': 1914, 'n_tokens_neox': 1773, 'n_words': 558, 'path': 'Week 3 - ODE linear systems.ipynb', 'repo_name': 'newby-jay/MATH371-Winter2020-JupyterNotebooks', 'size': '4578'}