离散模型-北太天元学习18

离散模型-北太天元学习18

我们前面讲到一个例子，是描述鸡蛋的温度随着时间的变化的，我们让T_n表示鸡蛋在第n个时刻 的温度, 通过递推关系 T_{n+1} = f( T_n ) 来建立模型。这是一个离散模型。在很多情况下，我们定期跟踪发生的事情，例如每个小时，每天或每年。将时间视为连续变量可能没有额外的好处。离散模型主要由递归关系和一些模型参数组成。

我们今天先以人口模型为例来讲讲离散模型，其实这里的人口可以是其他生物的数量，如草原上的兔子，池塘里的鱼，某个培养皿中的细菌或者病毒。我们先统一的说成人口。  

例如，设x_n是某个时间n的人口数量。递归关系可以以如下形式的方程给出

x_n=f(x_{n−1}，x_{n−2}，…，x_{1}，x_{0}).

也就是说，在时间n，x_n的是所有先前所有时间的人口数量的函数。
形式为x_n=f(x_{n−1})的简单方程被认为是一阶递推关系，因为它只依赖于上一个时间的值。 形式为x_n=f(x_{n−1}，x_{n−2}) 的序列被认为是二阶递归关系。

如果一个模型使用一个或多个递推关系来描述某些现象，则称其为离散动力系统。

首先，我们研究了一个基本的人口增长模型。假设人口在下一个时间间隔的增长率与当前人口成正比，这是一阶关系。让x_n是时间n的人口， 按照上面的假设，我们就有
      x_{n+1}−x_{n}=r*x_{n}，
其中r是是一个常数，称为增长率。左侧的x_{n+1}−x_{n}表示在在时间n到时间 (n+1)这一个时间单位的人口的增加，我们看到它与当前人口x_{n}成比例。

对x_{n+1}与x_n的递推关系就是如下的模型

x_{n+1} = x_n + r* x_n=（1+r)* x_n.

这有时被称为指数增长模型，或马尔萨斯模型，以托马斯·马尔萨斯（1766-1834）的名字命名。

注意: 当|1+r|<1 时, 这个模型给出的人口指数是随着时间是指数衰减的;
当1+r>1时，这个模型给出的人口随着时间是指数增长的。

通常情况下，如果种群数量增加到某一个数以上，这个模型就不再适用了， 例如，一个国家的人口接近一定的最大容量，增长率就会放缓。逻辑斯蒂克模型描述了 增长率随人口接近最大容量。设K是环境支持人口的最大容量 （这有时被称为人口的承载能力）。我们假设当人口x与K相比相对较小时， 人口呈指数级增长，并且随着人口数量向承载能力的方向增加，增长放缓。 离散逻辑斯蒂克模型是

x_{n+1} = x_n + r* x_n * ( 1 -  x_n / K)

其中K是承载能力，x_n是时间n时的人口。注意，当x_n与K相比小时，x_n/K将接近于零， 这意味着x_{n+1}≈x_n+r*x_n。这意味着当人口较少时，逻辑斯蒂克模型接近指数增长模型。 另一方面，当x_n接近K时，x_{n}/K≈1，因此x_{n+1}≈x_{n}+r*x_{n}(1-1) = x_n。 换句话说，随着人口接近承载能力K，人口是相对恒定的。

逻辑斯蒂克增长方程的另一种写法

x_{n+1}=（1+r）x_{n}− r * K * ( x_{n} )^2 .

虽然仍然是一阶，但在这种形式下，更明显的是，逻辑斯蒂克增长是一种非线性递推关系，
因为它包含二次项− r * K * ( x_{n} )^2.


如何确定模型参数，如增长率r, 承载能力K?

对于指数增长模型,  根据数据 x_{0}, x_{1}, 可以确定r:
具体来说就是根据 x_{1} = (1+r)*x_{0}, 得到 r = x_{1} /x_{0} - 1.

对于逻辑斯蒂克模型，根据数据 x_{0}, x_{1}, x_{2} , 可以确定 r, K  
x_{1}=（1+r）x_{0}− r * K * ( x_{0} )^2,  
x_{2}=（1+r）x_{1}− r * K * ( x_{1} )^2
联立可以求解可得 r 和 K.

我们现在引入一个新的词语: 态(state), 大家去百度搜索，你会发现在量子力学等很多学
科都在使用这个词语。我这里只是简单把他用成描述某个研究对象的状态， 例如我要研究的
是一个质点，它在时间 t 的状态可以用  x_t =  [q_t; p_t] 来表示, 这里 q_t 表示
时间 t 时的该质点的位置, p_t 表示 时间 t 时该质点的 动量。 这里 [q_t; p_t] 表示
一个列向量, q_t 和 p_t 之间的分号是北太天元 把q_t 和 p_t 做竖直方向的连接。
因此，这里的 态(state) 就是一个向量。

再举一个态的例子: 我们仅仅考虑在草原上生活着两种动物 兔子 和 狼, 我们要研究这
两个物种构成的一个系统，这个系统在时间n的态(state) 可以用这个两个物种的数量
组成的向量来表示, 也就是   
   x_t = [  兔_t ; 狼_t ]
其中 兔_t 表示 t 时刻的兔子的数量， 狼_t 表示时间t时狼的数量。  
下面给出一个简单的捕食者-被捕食者(predator-prey)离散模型：
   兔_{t+1}  = (1+r) * 兔_{t}  -  f * 狼_{t} * 兔_{t}
   狼_{t+1}  = (1+w) * 狼_{t}  +  g * 狼_{t} * 兔_{t}
这里出现了四个模型参数, r 表示兔子的自然生长率，w是狼的自然生长率,
 f 和 g 则是反应 狼兔相逢 对兔 和 狼 数量的影响。这些参数的确定可以通过
观测两个时间的态[兔_0; 狼_0], [兔_1; 狼_1]来确定。 

下面我们给出北太天元写的一个包含子函数的脚本计算兔狼数量随时间变化:
%狼是捕食者，兔子是被捕食者
%
%
%
clf;
close all;
clear all;


x = [ 1000; 1];
r = 0.01; f = 1e-4; w = 2e-3; g = 3e-5;
para = [ r; w; f; g];
x1 = pp_step(x,para);
x2 = pp_step(x1,para);

A = [ x(1),  0,    -x(1)*x(2), 0;
      0,     x(2), 0,          x(1)*x(2);
     x1(1),  0,    -x1(1)*x1(2), 0;
      0,     x1(2), 0,          x1(1)*x1(2) ] ;
b = [ x1 - x;
        x2 - x1];
para = A\b; % 模型参数求解
N = 100;
p = zeros(2,N+1);
p(:,1) = x;
for n=1:N
    p(:,n+1) = pp_step(p(:,n), para);  
end

subplot(1,2,1)
plot(0:N, p(1,:),'r-','LineWidth', 3)
title('兔子')
subplot(1,2,2)
plot(0:N, 100*p(2,:),'b-','LineWidth', 3)
title('狼')



function x1 = pp_step(x, para)
    [r, w, f, g] =  deal(para(1), para(2), para(3), para(4));
    x1 = x;
    x1(1) = (1+r)*x(1) - f*x(2)*x(1);
    x1(2) = (1+w)*x(2) + g*x(2)*x(1);
end

离散模型还有一个非常重要的应用，是谷歌(goole)搜索的网页排名模型(page rank model),
这个值得单独讲一次，因此我们留在下一节。