17、微分方程求解、种群动态模拟与统计分析

sre5engineer

于 2025-11-18 16:59:51 发布

阅读量3

点赞数

CC 4.0 BY-SA版权

分类专栏： MATLAB数值方法入门文章标签：微分方程求解边界值问题种群动态

本文链接：https://blog.youkuaiyun.com/sre5engineer/article/details/155172125

MATLAB数值方法入门专栏收录该内容

23 篇文章 ¥499.90

订阅专栏¥69.90

会员秒杀 ¥9.9 重磅福利

超级会员免费看

微分方程求解、种群动态模拟与统计分析

1. 微分方程求解

1.1 边界值问题

对于一般的边界值问题，若为线性问题，可直接将整个问题转化为矩阵问题求解。下面通过具体例子说明。

例 1：求解微分方程 $\frac{d^2y}{dx^2} - y = 0$，边界条件为 $y(0) = 0$ 和 $y(1) = 1$

步骤 1：设置网格点

x = linspace(0,1,11);
y = zeros(size(x));
h = x(2)-x(1);

步骤 2：推导离散方程
在点 $x = x_j$ 处，近似可得 $\frac{y_{j + 1} - 2y_j + y_{j - 1}}{h^2} - y_j = 0$，整理得 $y_{j + 1} + y_j(-2 - h^2) + y_{j - 1} = 0$，此方程适用于内部点 $j = 2$ 到 $j = 10$。
步骤 3：考虑边界条件并构建矩阵
在 $x = 0$ 处，$y_1 = 0$；在 $x = 1$ 处，$y_{11} = 1$。可得到如下矩阵方程：
$\begin{pmatrix}
1 & 0 & 0 & 0 & 0 & 0 & 0 & 0 & 0 & 0 & 0 \
1 & \tau & 1 & 0 & 0 & 0 & 0 & 0 & 0 & 0 & 0 \
0 & 1 & \tau & 1 & 0 & 0 & 0 & 0 & 0 & 0 & 0 \
0 & 0 & 1 & \tau & 1 & 0 & 0 & 0 & 0 & 0 & 0 \
0 & 0 & 0 & 1 & \tau & 1 & 0 & 0 & 0 & 0 & 0 \
0 & 0 & 0 & 0 & 1 & \tau & 1 & 0 & 0 & 0 & 0 \
0 & 0 & 0 & 0 & 0 & 1 & \tau & 1 & 0 & 0 & 0 \
0 & 0 & 0 & 0 & 0 & 0 & 1 & \tau & 1 & 0 & 0 \
0 & 0 & 0 & 0 & 0 & 0 & 0 & 1 & \tau & 1 & 0 \
0 & 0 & 0 & 0 & 0 & 0 & 0 & 0 & 1 & \tau & 1 \
0 & 0 & 0 & 0 & 0 & 0 & 0 & 0 & 0 & 0 & 1
\end{pmatrix}
\begin{pmatrix}
y_1 \
y_2 \
y_3 \
y_4 \
y_5 \
y_6 \
y_7 \
y_8 \
y_9 \
y_{10} \
y_{11}
\end{pmatrix}
=
\begin{pmatrix}
0 \
0 \
0 \
0 \
0 \
0 \
0 \
0 \
0 \
0 \
1
\end{pmatrix}$
其中 $\tau = -2 + h^2$。
步骤 4：使用 MATLAB 求解矩阵方程

A = diag(ones(11,1)*(-2-h^2),0)+diag(ones(10,1),1) +diag(ones(10,1),-1);
A(1,:) = 0; A(1,1) = 1;
A(11,:) = 0; A(11,11) = 1;
b = zeros(11,1);
b(11) = 1;
y = A\b;

此方法求解结果误差约为 $10^{-4}$。

例 2：求解微分方程 $\frac{d^2y}{dx^2} - y = 0$，边界条件为 $y’(0) = 0$ 和 $y(1) = 1$

同样的微分方程，只需修改边界条件的设置。在 $x = 0$ 处，由 $\frac{dy}{dx}|_{x = x_1} \approx \frac{y_1 - y_0}{h} = 0$ 可得 $y_0 = y_1$，则矩阵第一行设置为：

A(1,:) = 0;
A(1,1) = 1; A(1,2) = -1;

此方法求解结果误差约为 $10^{-2}$。

1.2 误差分析

在例 1 中，使用了近似公式 $\frac{d^2y}{dx^2}| {x = x_j} \approx \frac{y {j + 1} - 2y_j + y_{j - 1}}{h^2}$。通过泰勒级数展开：
$y(a + h) = y(a) + h\frac{dy}{dx}| {x = a} + \frac{h^2}{2!}\frac{d^2y}{dx^2}| {x = a} + \frac{h^3}{3!}\frac{d^3y}{dx^3}|_{x = a} + \cdots$
将泰勒级数代入近似公式可得误差与 $h^2$ 和区间内四阶导数的最大值成正比。

在例 2 中，在边界处使用了 $y’(0) \approx \frac{y_2 - y_1}{h}$，误差与 $h$ 成正比，所以误差比例 1 大。

1.3 种群动态

1.3.1 逻辑斯谛模型

考虑种群增长的逻辑斯谛模型 $\frac{dN}{dt} = rN - \frac{rN^2}{K}$，其中 $r$ 为正，$K$ 表示环境承载能力。其解析解为 $N(t) = \frac{N_0K}{N_0 + (K - N_0)e^{-rt}}$，当 $t \to \infty$ 时，$N(t) \to K$。
使用简单的欧拉积分方案进行数值求解：

r = 0.5;
K = 2;
N0 = 1;
t = 0:0.1:10.0;
delta_t = t(2)-t(1);
nt = length(t);
N = zeros(size(t));
N(1) = N0;
for j = 1:(nt-1)
    N(j+1) = N(j) + delta_t*(r*N(j)-r*N(j)^2/K);
end
exact = N0*K./(N0+(K-N0)*exp(-r*t));

数值解与解析解匹配良好。

1.3.2 多维系统

考虑捕食者 - 猎物模型：
$\frac{dN}{dt} = N(a - bP)$
$\frac{dP}{dt} = P(cN - d)$
其中 $P$ 表示捕食者数量，$N$ 表示猎物数量，$a$、$b$、$c$、$d$ 反映两物种的繁殖和捕食习性。使用欧拉积分求解：

a = 1; b = 2;
c = 1; d = 2;
N0 = 4;
P0 = 0.2;
t = 0:0.1:10.0;
delta_t = t(2)-t(1);
nt = length(t);
N = zeros(size(t));
P = zeros(size(t));
N(1) = N0;
P(1) = P0;
for j = 1:(nt-1)
    N(j+1) = N(j)+delta_t*(N(j)*(a-b*P(j)));
    P(j+1) = P(j)+delta_t*(P(j)*(c*N(j)-d));
end

该模型的解呈现振荡特性，猎物数量增长后，捕食者数量随之增加，当捕食者过多时，猎物数量减少，随后猎物数量又会恢复，循环往复。

1.4 微分系统的特征值

考虑微分方程 $\frac{d^2y}{dx^2} + \lambda y = 0$，边界条件为 $y(0) = y(\pi) = 0$。将方程离散化：
$-\frac{y_{j + 1} - 2y_j + y_{j - 1}}{h^2} = \lambda y_j$，$j = 2, \cdots, N - 1$，且 $y_1 = y_N = 0$，可得到矩阵方程 $Ay = \lambda y$，$\lambda$ 为矩阵 $A$ 的特征值。

n = 200;
h = pi/(n-1);
A = diag(2*ones(n-2,1)/h^2,0) +diag(-ones(n-3,1)/h^2,1) +diag(-ones(n-3,1)/h^2,-1);
[V,D] = eigs(A,3,'SM');
y = linspace(0,pi,n);
for i = 1:3
    eigf = [0 V(:,i)' 0];
    plot(y,eigf)
    pause
end

此代码可得到前三个特征值和特征函数，理论上特征值应为 1、4 和 9，对应特征函数为 $\sin x$、$\sin 2x$ 和 $\sin 3x$。

1.5 任务

以下是一些相关任务及要求：
1. 任务 1 ：修改代码求解微分方程 $\frac{dy}{dt} = -y\sqrt{t}$，边界条件为 $y(0) = 1$，$t$ 从 0 到 1，并求出解析解进行对比。
2. 任务 2 ：计算微分方程 $\frac{dy}{dt} = ty$ 和 $\frac{dy}{dt} = -ty$ 在边界条件 $y(0) = 1$，$t$ 从 0 到 1，$\Delta t = \frac{1}{4}$ 时的误差，并进行比较。
3. 任务 3 ：修改代码求解微分方程 $\frac{dy}{dt} = \sin t + \sin y$，边界条件为 $y(0) = 0$，$t$ 从 0 到 $10\pi$，步长小于 1。
4. 任务 4 ：使用隐式欧拉方案求解微分方程 $\frac{dy}{dt} = -\frac{t}{2} - \frac{y}{3}$，边界条件为 $y(0) = 0$，$\Delta t = \frac{1}{3}$，$t$ 从 0 到 1，并与解析解和使用 MATLAB 代码（1001 个点）的结果进行比较。
5. 任务 5 ：修改代码使用 MATLAB 例程 ode45 求解微分方程 $\frac{dy}{dt} = t^2 - y^2$，边界条件为 $y(0) = 0$，$t$ 从 0 到 2。
6. 任务 6 ：比较微分方程 $\frac{dy}{dt} = -y + t^2$ 在边界条件 $y(0) = 1$，$t$ 从 0 到 2 时的数值解和解析解。
7. 任务 7 ：修改代码求解以下问题：
- $y’’ = x\cos x$，边界条件为 $y(0) = 1$ 和 $y(1) = 0$。
- $y’’ = x\cos x$，边界条件为 $y’(0) = 0$ 和 $y(1) = 0$。
- $y’’ + 2y’ = 0$，边界条件为 $y(0) = 0$ 和 $y(1) = 1$。
8. 任务 8 ：求解微分方程 $\frac{d^2y}{dt^2} + 3y = t$，边界条件为 $y(0) = y’(0) = 0$，$t$ 从 0 到 1，进行数值和解析求解。
9. 任务 9 ：求解微分方程 $\frac{d^2y}{dt^2} + ty = \sin t$，初始条件为 $y(0) = 1$ 和 $y’(0) = 0$，$t$ 从 0 到 2；再使用打靶法或直接求解法，边界条件为 $y(0) = y(2) = 0$。
10. 任务 10 ：求解一阶微分方程 $(t^2 + 1)y’ + 2ty = 0$，初始条件为 $y(0) = 1$，$t$ 从 0 到 5，进行数值和解析求解。
11. 任务 11 ：将三阶微分方程 $y’‘’ - 2y’’ - y’ + 2y = 0$ 转化为三个一阶方程，求解初始条件为 $y(0) = y’(0) = 0$ 和 $y’‘(0) = 1$ 的情况。
12. 任务 12 ：使用有限差分法确定微分方程 $y’’ + (\sin x - \lambda)y = 0$ 在边界条件 $y(0) = y(\pi) = 0$ 下的特征值 $\lambda$，考虑绝对值最小的三个特征值对应的特征函数。

2. 模拟与随机数

2.1 统计量

2.1.1 平均值

常见的平均值有均值、中位数和众数。
- 均值：$\bar{x} = \frac{1}{N}\sum_{i = 1}^{N}x_i$，在 MATLAB 中可使用 mean(x) 计算，也可使用代码 xbar = sum(x)/length(x); 计算。
- 中位数 ：将数据排序后，若数据个数为奇数，中位数为中间值；若为偶数，中位数为中间两个数的平均值。在 MATLAB 中可使用 median(x) 计算，也可使用以下代码：

xs = sort(x);
lx = length(x);
if mod(lx,2)==0
    median = (xs(lx/2)+xs(lx/2+1))/2;
else
    median = xs((lx+1)/2);
end

例如，对于数据 x = [3 4 5 3 2 1 3 5 6 3 2 5 6 8 2 0]; ，均值为 3.625，中位数为 3。

2.1.2 其他统计量

方差：$\frac{1}{N - 1}\sum_{i = 1}^{N}(x_i - \bar{x})^2$，可表示为 $E(X^2) - E(X)^2$。在 MATLAB 中可使用 var(x) 计算，也可使用代码：

N = length(x);
xbar = sum(x)/N;
var = sum((x-xbar).^2)/(N-1);

标准差 ：方差的平方根，在 MATLAB 中使用 std(x) 计算。
协方差 ：$\sigma_{XY} = E(XY) - E(X)E(Y)$，对于一组数据，$\sigma_{XY} = \frac{1}{n}\sum_{i = 1}^{n}x_iy_i - \frac{1}{n^2}(\sum_{i = 1}^{n}x_i)(\sum_{i = 1}^{n}y_i) = \frac{1}{n}\sum_{i = 1}^{n}(x_i - \bar{x})(y_i - \bar{y})$。
相关系数 ：$r_{XY} = \frac{\sigma_{XY}}{\sigma_X\sigma_Y} = \frac{\sum_{i = 1}^{n}(x_i - \bar{x})(y_i - \bar{y})}{\sqrt{\sum_{i = 1}^{n}(x_i - \bar{x})^2\sum_{i = 1}^{n}(y_i - \bar{y})^2}}$，相关系数取值范围为 $[-1, 1]$，1 表示正相关，-1 表示负相关。

以下是一个简单的流程图，展示统计量计算的大致流程：

graph LR
    A[输入数据 x] --> B[计算均值]
    B --> C[计算中位数]
    B --> D[计算方差]
    D --> E[计算标准差]
    A --> F[输入数据 y]
    B --> G[计算协方差]
    D --> H[计算相关系数]

通过以上内容，我们介绍了微分方程的求解方法、种群动态模型的模拟以及一些统计量的计算，这些知识在科学研究和工程应用中都有广泛的应用。希望大家通过实践这些任务，能更好地掌握相关知识和技能。

2.2 示例演示

2.2.1 计算均值、中位数、方差和标准差

下面通过具体代码展示如何计算给定数据的均值、中位数、方差和标准差。

x = [3 4 5 3 2 1 3 5 6 3 2 5 6 8 2 0];
% 计算均值
xbar = sum(x)/length(x);
disp(['均值: ', num2str(xbar)]);
% 计算中位数
xs = sort(x);
lx = length(x);
if mod(lx,2)==0
    median = (xs(lx/2)+xs(lx/2+1))/2;
else
    median = xs((lx+1)/2);
end
disp(['中位数: ', num2str(median)]);
% 计算方差
N = length(x);
xbar = sum(x)/N;
var = sum((x-xbar).^2)/(N-1);
disp(['方差: ', num2str(var)]);
% 计算标准差
std_dev = sqrt(var);
disp(['标准差: ', num2str(std_dev)]);

运行上述代码，输出结果如下：
|统计量|值|
| ---- | ---- |
|均值|3.625|
|中位数|3|
|方差|4.2232|
|标准差|2.055|

2.2.2 计算协方差和相关系数

假设有两组数据 x 和 y ，以下代码展示如何计算它们的协方差和相关系数。

x = [1 2 3 4 5];
y = [2 4 6 8 10];
n = length(x);
xbar = sum(x)/n;
ybar = sum(y)/n;
% 计算协方差
covariance = 1/n * sum((x - xbar).*(y - ybar));
disp(['协方差: ', num2str(covariance)]);
% 计算相关系数
std_x = sqrt(1/(n - 1) * sum((x - xbar).^2));
std_y = sqrt(1/(n - 1) * sum((y - ybar).^2));
correlation = covariance / (std_x * std_y);
disp(['相关系数: ', num2str(correlation)]);

运行上述代码，输出结果如下：
|统计量|值|
| ---- | ---- |
|协方差|2.5|
|相关系数|1|

2.3 任务实践

以下是一些与统计量计算相关的任务：
1. 任务 1 ：生成一组包含 100 个随机数的数据，计算其均值、中位数、方差和标准差。

% 生成 100 个随机数
x = rand(1, 100);
% 计算均值
mean_x = mean(x);
% 计算中位数
median_x = median(x);
% 计算方差
var_x = var(x);
% 计算标准差
std_x = std(x);
disp(['均值: ', num2str(mean_x)]);
disp(['中位数: ', num2str(median_x)]);
disp(['方差: ', num2str(var_x)]);
disp(['标准差: ', num2str(std_x)]);

任务 2 ：生成两组不同的随机数据，计算它们的协方差和相关系数，并分析结果。

% 生成两组随机数据
x = rand(1, 100);
y = rand(1, 100);
% 计算协方差
covariance = cov(x, y);
% 计算相关系数
correlation = corrcoef(x, y);
disp(['协方差: ', num2str(covariance(1,2))]);
disp(['相关系数: ', num2str(correlation(1,2))]);

任务 3 ：根据给定的数据，分析不同统计量之间的关系，例如均值和中位数的差异与数据分布的关系。
可以通过绘制数据的直方图来直观观察数据分布，以下是示例代码：

x = [3 4 5 3 2 1 3 5 6 3 2 5 6 8 2 0];
figure;
histogram(x);
title('数据直方图');
xlabel('数据值');
ylabel('频数');

2.4 总结

通过上述内容，我们详细介绍了微分方程求解的多种方法，包括边界值问题的求解、误差分析，以及种群动态模型的模拟。同时，还深入探讨了统计量的计算，如均值、中位数、方差、标准差、协方差和相关系数等，并通过具体的代码示例和任务实践帮助大家更好地理解和掌握这些知识。

在实际应用中，我们可以根据具体问题选择合适的方法进行求解和分析。例如，对于微分方程求解，要根据边界条件和方程类型选择合适的数值方法；对于统计分析，要根据数据特点和分析目的选择合适的统计量。

希望大家通过阅读本文，能够在实际工作和学习中灵活运用这些知识和技能，解决相关问题。

以下是一个综合流程图，展示从数据输入到微分方程求解和统计分析的整体流程：

graph LR
    A[输入数据] --> B{数据类型}
    B -->|微分方程数据| C[微分方程求解]
    C --> D[误差分析]
    B -->|统计数据| E[统计量计算]
    E --> F[均值计算]
    E --> G[中位数计算]
    E --> H[方差计算]
    E --> I[标准差计算]
    E --> J[协方差计算]
    E --> K[相关系数计算]
    F --> L[结果分析]
    G --> L
    H --> L
    I --> L
    J --> L
    K --> L
    D --> L

总之，微分方程求解和统计分析是科学研究和工程应用中非常重要的工具，掌握这些知识和技能将有助于我们更好地解决实际问题。