分叉与逻辑斯蒂映射

分叉与逻辑斯蒂映射

上节回顾：

logistic方程：$$\{\begin{array}{ll}微分式：& N^{\prime }=r\cdot N(1-\frac{N}{K})\\ 积分式：& N_t=\frac{K\cdot b}{e^{-r\cdot t}+b}\end{array}$$

一、系统的概念

系统（System）：多元素有序排列、相互关联的功能单位（个人理解）。e.g.生态学系统 vs. 生态系统
系统一般具有三大要素：
$$\{\begin{array}{l}1、元素（Elements）\\ 2、关系（Relationship）\\ 3、功能（Function）\end{array}$$
此时的logistic方程 $N^{\prime }=r\cdot N(1-\frac{N}{K})$ 就是一个微分动力系统：
$$\{\begin{array}{l}左边=N^{\prime }=\frac{dN}{dt}：包含有参数t，是时变系统\\ 右边=r\cdot N(1-\frac{N}{K})：不包含t参数，为自治系统\end{array}$$
对于常微分方程（未知函数是一元函数）的求解，一般分为以下步骤：

1、Find a steady state

即寻求“静止”状态下的函数关系，e.g. logistic方程：$N^{\prime }=r\cdot N(1-\frac{N}{K})$ 在 $N^{\prime }=0$ 时，系统静止，即种群大小没有波动，即：
$$t时刻后种群大小的变化：N^{\prime }=\frac{dN}{dt}=r\cdot N(1-\frac{N}{K})=0$$

2、求解

即求解“静止”状态的函数方程，即：
$$\frac{dN}{dt}=r\cdot N(1-\frac{N}{K})=0$$
$$解：\{\begin{array}{l}{N}_1^{\ast }=0\\ N_2^{\ast }=K\end{array}$$
在数轴表示为：

3、寻求均衡解

$$N^{\prime }=r\cdot N(\frac{K-N}{K})$$
一共有以下几种情况：$$\{\begin{array}{l}N>K时，\frac{K-N}{K}<0，N\uparrow \\ N=K时，\frac{K-N}{K}=0，N不变\\ 0<N<K时，\frac{K-N}{K}>0，N\downarrow \\ N\le 0，无现实意义\end{array}$$
在数轴上表示为：

$\because N\to K$
$\therefore N_2^{\ast }=K是函数的均衡解$

二、logistic方程二次积分

$$\begin{array}{rl}{N}^{\prime \prime }& =(r\cdot N(1-\frac{N}{K}){)}^{\prime }\\ & =(r\cdot N-\frac{r\cdot N^2}{K}{)}^{\prime }\\ & =(r\cdot N{)}^{\prime }-(\frac{r\cdot N^2}{K}{)}^{\prime }\\ & =r-\frac{2\cdot N}{K}\end{array}$$
求解：
$令N^{\prime \prime }=0得：$
$$0=r-\frac{2r\cdot N}{K}$$
解得$N=\frac{K}{2}$，即收获理论

$N^{\prime \prime }=r-\frac{2r\cdot N}{K}，有几种情况：$
$$\{\begin{array}{l}{N}^{\prime \prime }>0时，N<\frac{K}{2}，N^{\prime }上升（上凹）\\ N^{\prime \prime }<0时，N>\frac{K}{2}，N^{\prime }下降（下凹）\\ N^{\prime \prime }=\frac{K}{2}时，N^{\prime }最大，种群增长率最大\end{array}$$

三、Bifurcation and Logistic Map

1、$N^{\prime }=4N-N^3+C$

构建函数：
$$N^{\prime }=4N-N^3$$
求解：
Steady state:$N^{\prime }=0$
方程 $4N-N^3=0$ 解得：
$$\{\begin{array}{l}{N}_1=-2\\ N_2=0\\ N_3=2\end{array}$$
函数关系可视化为：

$此时函数N^{\prime }=4N-N^3与N轴有三个交点N_1、N_2、N_3。$
$\because$ $N\in (0,2]时，N^{\prime }最大值为a$
$\therefore$ $当C=a时，函数与N轴有两个交点$

使用Plotly（python v.3.9）绘制散点图，以探究 $a$ 值：

from plotly.graph_objs import Scatter
from plotly import offline

#Three lists made visualization requested.
dNs, Ns, labels = [], [], []

#Make data lists
for N in range(-30,30):
    N = N * 0.1
    dN = 4 * N - N**3
    label = f"({N},{dN})"
    Ns.append(N)
    dNs.append(dN)
    labels.append(label)

#visulization
trace =Scatter(
    x = Ns, 
    y = dNs, 
    mode = "markers",  
    marker = dict(color = 'rgba(0, 255, 0, 1)'),
    text= labels
    )

data = [trace]
layout = dict(title = "N' = 4N - N^3",
              xaxis= dict(title= "N"),
              yaxis= dict(title= "N'")
             )
fig = dict(data = data, layout = layout)
offline.plot(fig,filename='N_dN.html')

从图中可以看出，当 $N\approx 1.2时，N^{\prime }\approx 3.072，即a\approx 3.072$
$为探究函数与N轴的交点数的变化，引入参数C得新函数：$
$$N^{\prime }=4N-N^3+C$$
$使C值逐渐增加，使函数曲线与N轴的交点逐渐减少：$

Python实现：

from plotly.graph_objs import Scatter
from plotly import offline

#Three lists made visualization requested.
Ns, data = [], []
C = 0

#Make Ns
for N in range(-250,250):
    N = N * 0.01
    Ns.append(N)

#Main loop
while C <= 5:
    dNs = []
    for N in Ns:
        dN = (4 * N - N**3) + C
        dNs.append(dN)

    #Mark function curve, if C > 3.072
    if C > 3.072:
        trace =Scatter(
            x = Ns, 
            y = dNs, 
            mode = "markers",  
            marker = dict(color = 'rgba(255, 0, 0, 1)'),
            )
    #Hight the point by green, if C = 0
    elif C == 0:
        trace =Scatter(
            x = Ns, 
            y = dNs, 
            mode = "markers",  
            marker = dict(color = 'rgba(0, 255, 0, 1)')
            )
    #Make points transparencied, if C isn't 3.072 or 0
    else:
        trace =Scatter(
            x = Ns, 
            y = dNs, 
            mode = "markers",  
            marker = dict(color = 'rgba(0, 0, 255, 0.1)')
            )
    data.append(trace)

    #iteration
    C += 0.1512

    
layout = dict(title = "N' = 4N - N^3 + C",
              xaxis= dict(title= "N"),
              yaxis= dict(title= "N'")
             )
fig = dict(data = data, layout = layout)
offline.plot(fig, filename = "N_dNs.html")

从可视化结果中可以发现，随着C的值增加，交点 $N_1$ 和 $N_2$ 越来越近，当 $C=a$ 时，$N_1=N_2$ 即两交点合并为一个交点。当 $C>a$ 时，函数与N轴只余一个交点 $N_3$ 。

2、Bifurcation

现在以 $C$ 为 $x$ 轴，$N$ 为 $y$ 轴：

此时 $C=a$ 时对应的 $N$ 值就是Critical point，即在此点发生了临界突变（Sudden Transition）。

3、微分原理

以logistic方程为例：
$$\begin{array}{rl}\frac{dN}{dt}& =r\cdot N(1-\frac{N}{K})\\ & =\underset{\Delta t\to 0}{\lim }\frac{N_{t+\Delta t}-N_t}{(t+\Delta t)-t}\end{array}$$
表示单位时间 ($\Delta t$,时间增量)内，种群大小 $N$ 的增量 ($\Delta N$ )的近似值（$\frac{dN}{dt}$）,但需要注意，$\Delta N\ne \frac{dN}{dt}$，后者只是前者的近似值，且随着$\Delta t\to 0$,二者间的差异也趋近于零。

4、Logistic映射

整理合并，得：
$$\frac{N_{t+\Delta t}-N_t}{\Delta t}=r\cdot N_t(1-\frac{N_t}{K})$$
令$\Delta t=1$（1表示1单位的时间增量，不是具体时间间隔）得：
$$N_{t+1}-N_t=r\cdot N_t(1-\frac{N_t}{K})$$
整理等式：
$$\begin{array}{rl}{N}_{t+1}& =N_t+r\cdot N_t(1-\frac{N_t}{K})\\ & =N_t[1+r(1-\frac{N_t}{K})]\\ & =N_t(r+1-\frac{r\cdot N_t}{K})\\ & =(r+1)N_t\cdot [1-\frac{r\cdot Nt}{(r+1)K}]\end{array}$$
设：$$\begin{array}{rl}& x=\frac{r\cdot N}{(r+1)K}\\ & x_t=\frac{r\cdot N_t}{(r+1)K}\\ & x_{t+1}=\frac{r\cdot N_{t+1}}{(r+1)K}\end{array}$$
将$x,x_t,xt+1$带入得：
$$\begin{array}{rl}{N}_{t+1}& =(r+1)N_t\cdot [1-\frac{r\cdot N_t}{(r+1)K}]\\ & =(r+1)N_t\cdot (1-x_t)\\ \frac{x_{t+1}\cdot (r+1)K}{r}& =\frac{x_t\cdot (r+1{)}^2\cdot K\cdot (1-x_t)}{r}\\ x_{t+1}& =x_t\cdot (r+1)\cdot (1-x_t)\end{array}$$
设 $\mu =r+1$ 得：
$$x_{t+1}=x_t\cdot \mu \cdot (1-x_t)$$
$\because$
$$N>0，r\in (0,1]，K>0，K>N$$
$\therefore$
$$x_t\in (0,1]，\mu \in (0,4)$$

5、映射可视化

据上述函数关系，利用python实现函数的可视化：

from tqdm import tqdm
import numpy as np
from plotly.graph_objs import Scatter
from plotly import offline


def LogisticMap():
    mu = np.arange(0, 4, 0.0001)
    x = 0.2  #initial value of x
    iters = 1000  # iterations that not be visualized
    last = 100  # iterations will be visualized

    #iteration
    for i in tqdm(range(iters+last)):
        x = mu * x * (1 - x)
        if i >= iters:
            #visualization
            trace =Scatter(
                x = mu, 
                y = x, 
                mode = "markers", 
                marker = dict(color = 'rgba(125, 255, 125, 1)',
                    size = 2)
                )
    data = [trace]
    layout = dict(title = "Logistic Map",
              xaxis= dict(title = "mu"),
              yaxis= dict(title = "x")
             )
    fig = dict(data = data, layout = layout)
    offline.plot(fig,filename = 'Logisticmap.html')


LogisticMap()

plotly可能存在兼容性的问题（更有可能是我使用不当），图中的某些点没有显示。

利用另一常用的模块matplotlib进行可视化：
模块一：

from tqdm import tqdm
import matplotlib.pyplot as plt
import numpy as np


def logisticmap(mu = 4):
    mu = np.arange(2, mu, 0.0001)
    x = 0.2  
    iters = 1000 
    last = 100 
    for i in tqdm(range(iters+last)):
        x = mu * x * (1 - x)
        if i >= iters:
            plt.plot(mu, x, ',k', alpha=0.25)
    plt.show()

可视化：

from image import logisticmap

logisticmap(mu = 4)