当前位置：首页 > news >正文

ODE-by-Matlab-01-人口增长模型

news 2025/8/15 12:14:05

博客地址：Matlab微分方程01-模型

马尔萨斯模型

马尔萨斯模型是人口增长模型中最简单的模型，它由英国牧师家马尔萨斯在1798年提出。
他利用在教堂工作的机会，收集英国100多年的人口数据，发现人口的相对增长率是常数。
在这个基础上，建立了一个描述人口增长的模型，也就是著名的“马尔萨斯人口模型”。

在这个模型中，最重要的额概念是相对增长率。

$u˙u=r\frac{\dot{u}}{u} = r$

这里就涉及到微分的概念，变量 $u$ 对时间 $t$ 的导数，也就是变化率(这里可以称之为增长率)，记作 $u˙\dot{u}$ 。

$u˙=dudt\dot{u} = \frac{du}{dt}$

那么“马尔萨斯人口模型”表达为微分方程是这样的形式。

$dudtu=α→dudt=αu\frac{\frac{du}{dt}}{u} = \alpha \rightarrow \frac{du}{dt} = \alpha u$

其中，时间为 $t$ ，人口 $u$ 为依赖于时间的函数，相对增长率是 $α\alpha$ （ $α>0\alpha >0$ ）。

这个方程的解很容易通过不定积分求解。

$∫duu=∫αdt→ln⁡u=αt+C→u(t)=eαt+C=neαt\int \frac{du}{u} = \int \alpha dt \rightarrow \ln u = \alpha t + C \rightarrow u(t) = e^{\alpha t + C} = n e^{\alpha t}$

这个解是一个指数函数！众所周知，指数函数的增长是非常快的。这在一定的程度上导致了社会主义国家考虑对人口增长进行控制。

n = 81;
t = linspace(0, 100, 100);
alpha = 0.02; % 人口增长率
% 马尔萨斯人口模型
u = n * exp(alpha * t);
figure;
plot(t, u , 'b-', 'LineWidth', 2);
xlabel('时间(年)');
ylabel('人口(亿)');
title('马尔萨斯人口模型')
grid on;
% 保存图像  
exportgraphics(gcf, 'chp01/malthus_population_model.png', 'Resolution', 300);

在这里插入图片描述

UN预测数据

我们可以从UN官网下载人口预测数据。
数据的格式是csv，第二表格包含了世界各国从2025年到2100年的人口预测数据。

我们通过websave命令下载数据，并用readtable命令读取数据。注意实际的数据位置、文件名可能会发生改变，需要根据实际情况修改。

UN_Projection_File = 'https://population.un.org/wpp/assets/Excel%20Files/1_Indicator%20(Standard)/EXCEL_FILES/1_General/WPP2024_GEN_F01_DEMOGRAPHIC_INDICATORS_COMPACT.xlsx';
% 下载UN人口预测数据
if ~exist('chp01/WPP2024_GEN_F01_DEMOGRAPHIC_INDICATORS_COMPACT.xlsx', 'file')websave('chp01/WPP2024_GEN_F01_DEMOGRAPHIC_INDICATORS_COMPACT.xlsx', UN_Projection_File);
end% 读取UN人口预测数据, sheet 'Medium variant', range 'L18:94': total population, K18:94: Years
if ~exist('chp01/WPP2024_GEN_F01_DEMOGRAPHIC_INDICATORS_COMPACT.xlsx', 'file')error('UN人口预测数据文件未找到，请检查下载路径。');
end
data = readtable('chp01/WPP2024_GEN_F01_DEMOGRAPHIC_INDICATORS_COMPACT.xlsx', ...'Sheet', 'Medium variant', 'Range', 'K18:L94');

然后我们来看看，马尔萨斯模型数据与UN预测数据的差异。

% 提取年份和人口数据
years = data{:, 1};
population = data{:, 2} * 1e3 * 1e-8; % 单位为千人-->亿% 绘制UN人口预测数据
figure; 
plot(years, population, 'k-', 'LineWidth', 2);
xlabel('年份');
ylabel('人口(亿)');
title('联合国人口预测');
grid on;
legends = {'UN Population Projection'};hold on;
% 绘制拟合的线性模型alphas = linspace(0.001, 0.01, 5); % 人口增长率范围for idx = 1:length(alphas)alphai = alphas(idx);% 使用马尔萨斯模型拟合UN人口预测数据iyears = years - years(1); % 将年份归一化u = population(1) * exp(alphai * iyears); % 初始人口为第一个年份的人口    % 绘制拟合曲线plot(years, u, 'DisplayName', ['alpha = ', num2str(alphai)], 'LineWidth', 1.0);legends{end + 1} = ['alpha = ', num2str(alphai)]; %#ok<*SAGROW>
endlegend(legends, 'Location', 'best');% 保存图像
exportgraphics(gcf, 'chp01/un_population_projection.png', 'Resolution', 300);

外链图片转存失败,源站可能有防盗链机制,建议将图片保存下来直接上传

增长率从0.001到0.01的马尔萨斯模型预测数据跟UN的预测有很大的区别，从图上看趋势都是不正确的。
联合国的预测中，人口会达到一个饱和值（峰值），并在其后缓慢下降。

本质上来讲，马尔萨斯模型仅仅考虑人口相对增长率的线性特征，没有考虑非线性的饱和特征。也就是，在人口较少时，人口的增长所受的限制很少，能够出现指数增长；而当人口达到一定的数量是，生存环境、生态资源、社会因素等都会对人口增长产生限制，导致人口增长率逐渐减小，最终趋近于0。

但是，我们并不能说马尔萨斯模型是错误的，因为马尔萨斯模型是线性的，它实际上能够很好的预测早期人口增长（例如英国工业革命带来的人口增长），实际上，从前面的UN数据和马尔萨斯模型拟合曲线中，我们也可以看到，马尔萨斯模型能够很好的预测早期人口增长。

这反映了微分方程模型的适用范围问题。在其适用条件和假设成立的范围内，模型能够准确描述系统的动态行为；但当系统偏离这些基本假设时，模型的预测能力将显著降低。这是数学建模中普遍存在的局限性。

基于以上分析，我们需要构建一个能够刻画人口增长非线性特征的数学模型，以更准确地描述人口动态变化规律。这就引出了下面将要介绍的Logistic模型。

Logistic模型

Logistic模型为什么叫做Logistic模型呢？因为它的解是一个Logistic函数。什么叫Logistic函数呢？它的形式是这样的：

$\frac{L}{1 + e^{-k(x - x_0)}}$

这里 $L$ 是函数的最大值， $k$ 是增长率， $x_0$ 是函数的中点。

figure;
% Plot Logistic function
t = linspace(0, 100, 100);
L = 10; % 最大人口
K = 0.1; % 资源限制plot(t, L ./ (1 + exp(-K * (t - 50))), 'r-', 'LineWidth', 2);
xlabel('时间(年)');
ylabel('人口(亿)');
title('Logistic Population Model');
grid on;% 保存图像
exportgraphics(gcf, 'chp01/logistic_population_model.png', 'Resolution', 300);

外链图片转存失败,源站可能有防盗链机制,建议将图片保存下来直接上传

至于这个函数为什么叫做Logistic函数呢？

因为法国数学家 Pierre François Verhulst就是这么命名的，他在1845年的论文（居然还有电子版可以看……）提出了这个函数来描述人口增长。这个论文现在还可以看到电子版，是法语的。在这论文电子版的地54也，作者比较了他的模型和马尔萨斯模型，并在图上给出了Logarithmique vs. Logistique的曲线命名。

在这里插入图片描述

Verhulst在论文里面写到：Nous donnerons le nom de logistique à la courbe (voyez la figure) caractérisée par l’équation précédente. 就是什么我们给前述方程给出的曲线（如图）命名为Logistic曲线之类……

说到底就是他给这个函数起了个名字……

我看了好多地方，都没有找到为什么叫logistique……

反正这个跟后勤（logistics）没有关系。国外网友猜测这个跟logistique的词源有关，logistique是法语，源自拉丁语logisticus，意思是“计算的”，而这个函数的计算……我编不出来了……

Verhulst的原始推导

现在我们来看看Verhulst是如何推导出Logistic模型的，这里我们遵循他在1845年论文中的原始推导。

Verhulst首先称，马尔萨斯的 $Mdppdt=l\frac{Mdp}{pdt}=l$ 并不正确，他给出一个修正的增长方程。

$Mdppdt=l−n(p−b)(2)\frac{M dp}{p dt} = l - n(p - b) \tag{2}$

这里 $M$ , $l$ , $n$ , $b$ 都是常数， $p$ 是人口， $t$ 是时间。

他接着写道：d’où, en posant, pour abréger, $m = l + nb$ （为了简化，设 $m = l + nb$ ），

$Mdppdt=m−np\frac{M dp}{p dt} = m - np$

然后得到：

$\frac{M dp}{mp - np^2} \tag{3}$

这个方程经过积分后（Cette équation étant intégrée donne），在观察到当 $t = 0$ 时对应 $p = b$ 的条件下，得到：

$\frac{1}{m} \log_e \left[ \frac{p(m - nb)}{b(m - np)} \right] \tag{4}$

这就是Verhulst在他的原始论文中推导出的积分形式的解。

从这个积分方程可以反解得到人口 $p$ 关于时间 $t$ 的显式表达式，这就是我们今天所知的Logistic函数。

传染病问题

% 编制一个在二维网格中的小老鼠染病模拟
% M只老鼠，N只有传染病，可以通过接触传染给健康的老鼠，
% 老鼠在[0,1]x[0,1]的二维网格中随机分布，并且随机运动
% 当二者距离小于0.01时，健康的老鼠有0.1的概率被感染
M = 500; % 老鼠数量
N = 1; % 传染病数量
% 初始化老鼠位置和状态
positions = rand(M, 2); % 老鼠位置
states = zeros(M, 1); % 0: 健康, 1: 感染
infected_indices = randperm(M, N); % 随机选择N只老鼠
states(infected_indices) = 1; % 设置感染状态% 模拟老鼠运动和传染
num_steps = 140; % 模拟步数history_infected = zeros(num_steps, M); % 记录每一步的感染状态
history_infected(1, :) = states; % 初始状态
position_history = zeros(num_steps, M, 2); % 记录每一步的位置
position_history(1, :, :) = positions; % 初始位置for step = 2:num_steps% 更新老鼠位置positions = positions + randn(M, 2) * 0.04; % 随机运动positions = mod(positions, 1); % 保持在[0,1]范围内% 检查感染传播for i = 1:Mif states(i) == 0 % 如果是健康的老鼠% 找到所有感染的老鼠infected_positions = positions(states == 1, :);distances = sqrt(sum((infected_positions - positions(i, :)).^2, 2));% 如果距离小于0.01且有感染概率，则感染if any(distances < 0.01) && rand < 0.8states(i) = 1; % 感染endendendhistory_infected(step, :) = states; % 记录当前状态position_history(step, :, :) = positions; % 记录当前位置信息
endfigure;
plot(1:num_steps, sum(history_infected, 2), 'b-', 'LineWidth', 2);
xlabel('时间步');
ylabel('感染老鼠数量');
title('老鼠传染病传播模拟');
grid on;exportgraphics(gcf, 'chp01/mouse_infection_simulation.png', 'Resolution', 300);

外链图片转存失败,源站可能有防盗链机制,建议将图片保存下来直接上传

实际上，无聊的话，还可以整一个老鼠传染病传播模拟，代码如下：

% 绘制初始状态
states = history_infected(1, :); % 获取初始状态
positions = squeeze(position_history(1, :, :)); % 获取初始位置
% 创建图形窗口
figure("Visible", 'off'); % 显示图形窗口
scatter(positions(:, 1), positions(:, 2), 50, states, 'filled');
title(['老鼠感染试验 - 时间 ', num2str(1), ' 感染老鼠数量: ', num2str(sum(states)), '/', num2str(M)]);
grid off;
box on; % 关闭网格线
axis off; % 关闭坐标轴显示
xlim([0 1]);
ylim([0 1]);
% 保存初始状态图像
fn = 'chp01/mouse_initial_distribution.gif';
if exist(fn, 'file')delete(fn); % 删除旧的文件
end
resolution = 100;
exportgraphics(gcf, fn , 'Resolution', resolution, 'Append', false);for step = 2:num_stepsstates = history_infected(step, :); % 获取当前步的状态positions = squeeze(position_history(step, :, :)); % 获取当前步的位置% 绘制当前状态scatter(positions(:, 1), positions(:, 2), 50, states, 'filled');title(['老鼠感染试验 - 时间 ', num2str(step), ' 感染老鼠数量: ', num2str(sum(states)), '/', num2str(M)]);grid on;box on; % 关闭网格线axis off; % 关闭坐标轴显示xlim([0 1]);ylim([0 1]);% 保存每一步的图像exportgraphics(gcf, fn, 'Resolution', resolution, 'Append', true);
end