数学生物学-4-混乱系统(Chaotic Systems)

tRNA做科研

于 2024-08-23 11:05:13 发布

阅读量1.1k

点赞数 30

分类专栏：计算生物学文章标签：算法 r语言 python 数学建模混沌系统计算生物学

本文链接：https://blog.csdn.net/qq_51566832/article/details/141457634

版权

计算生物学专栏收录该内容

6 篇文章 0 订阅

订阅专栏

混沌系统（Chaotic Systems）是一种复杂且非线性的动态系统，其基本概念和特征有：

对初始条件的敏感依赖性：
混沌系统的核心特征之一是对初始条件的极端敏感性，这一特性通常被称为“蝴蝶效应”。这意味着即使是极其微小的差异在初始状态也可能随时间放大，导致系统行为的显著不同。例如，天气预报模型中微小误差的累积可能导致长期预测的显著偏差。这种敏感性使得混沌系统的确切长期行为变得难以预测。
内在随机性：
混沌系统表现出一种看似随机的行为，但实际上这种行为是由确定性方程所描述的。换句话说，尽管系统的演化遵循严格的数学规则，但其表现出来的模式却类似于随机过程。这种内在的随机性并不意味着外部噪声或不确定性对系统的影响，而是系统本身固有的性质。
长期预测的不可能性：
由于混沌系统对初始条件的敏感依赖性和内在的随机性，这些系统通常被认为在长期内是不可预测的。即使我们拥有系统的完整数学描述和精确的初始条件，随着时间的推移，预测误差的增长也会导致长期预测失去意义。因此，对于混沌系统，只能进行有限时间内的可靠预测。

Robert May's Model

罗伯特·梅模型（Robert May's Model）是由理论生态学家罗伯特·梅在20世纪70年代中期提出的一种生态学模型。该模型的数学表达式非常简单，但却能描述复杂的生态动态行为。

模型的方程如下：

其中，Nt表示在时间t时的人口规模，以承载能力的比例表示，通常介于0和1之间。参数a是一个正数，影响人口的增长率。

这个方程虽然简单，但包含了丰富的生态动态信息。项aNt模拟了人口的指数增长特性，而(1−Nt)因子引入了密度依赖效应，意味着当人口规模接近承载能力（即Nt接近1）时，增长率会下降。

参数a的值对模型的行为至关重要。它的取值范围必须在0到4之间。当a小于1时，人口将趋向于灭绝。如果a在1到3之间，人口将稳定在一个特定水平，显示出稳定的平衡状态。然而，当a大于3时，系统将变得混沌，人口规模随时间表现出不可预测的行为。

我们使用R语言进行可视化：

# 定义罗伯特·梅模型的函数
robert_may_model <- function(N, a, t_max) {
  Nt <- numeric(t_max)
  Nt[1] <- N
  for (t in 1:(t_max - 1)) {
    Nt[t + 1] <- a * Nt[t] * (1 - Nt[t])
  }
  return(Nt)
}

# 设置初始条件
N0 <- 0.5 # 初始人口比例
t_max <- 100 # 时间步长

# 设置不同的a值
a_values <- c(1.5, 2.2, 3.5, 3.8)

# 初始化绘图
plot(0, type = "n", xlim = c(1, t_max), ylim = c(0, 1),
     xlab = "Time (t)", ylab = "Population Size (N)",
     main = "Robert May's Model for Different Values of a")

# 对每个a值进行模拟并绘制结果
for (a in a_values) {
  Nt <- robert_may_model(N0, a, t_max)
  lines(1:t_max, Nt, type = "l", col = rainbow(length(a_values))[which(a == a_values)])
}

# 添加图例
legend("topright", legend = paste("a =", a_values), col = rainbow(length(a_values)), lty = 1)

当参数a取较小值时，模型的解与离散逻辑模型相似；而当a的值增加时，解开始在两值之间振荡；对于更大的a值，解的振荡似乎没有特定模式。

模型求解

我们使用python对模型进行逐步求解，以更细致角度：

# May's Chaotic system population model 
import numpy as np 
import matplotlib.pyplot as plt 
a = 1.5 
steps = 15 
N = np.zeros(steps+1) 
t = np.zeros(steps+1) 
N[0] = 0.05 
t[0] = 0.0 
for i in range(0, steps): 
    N[i+1] = a*N[i]*(1 -N[i]) 
    t[i+1] = t[i] + 1.0 
print (N) 
# plot results 
plt.plot(t,N) 
plt.xlabel('Time') 
plt.ylabel('N') 
plt.show()

我们还可以使用Matlab得到同样的结果：

 % Script to run the solution 
steps = 20; 
N = zeros(1,steps); 
t = zeros(1,steps); 
N(1) = 0.05; 
for i = 1:1:steps 
N(i+1) = maysfunc(N(i)); 
t(i+1) = t(i) + 1; 
end 
plot(t,N); 
title('Mays model solution for a = 1.5') 
xlabel('Time') 
ylabel('N') 
function [ nnew ] = maysfunc( n ) 
a = 1.5; 
nnew = a*n*(1-n); 
end

模型固定点

固定点是数学模型中的一个关键概念，特别是在研究动态系统如May模型时。

数学生物学-3-固定点、稳定性和蛛网图(Fixed Points, Stability, and Cobwebbing)-CSDN博客

固定点是指系统中的一个状态，在该状态下，如果系统从该点出发，它将保持不变。在May模型中，固定点可以通过解析方法确定，其方程为：

通过代数变换，我们可以找到固定点的表达式：

这表明，固定点N的值取决于参数a。例如，当a=3.8时，固定点N计算如下：

我们使用mathematica求解固定点：

这个量是一个有理数，但当用浮点变量表示时，不能用有限位数完美地表示。使用一个能考虑32位数字的计算器进行一些计算，得到的值是0.7368421052631578947368...，在小数点后第18位开始重复。一个典型的Python（或C语言）程序中的float或double类型的变量只能保留大约16位数字。因此，在Mathematica笔记本的第一部分（图6.5）中，解再次开始振荡并不令人惊讶，因为我们无法完美地表示这个固定点的值。然而，请注意笔记本的第二部分，其中初始值保持为分数28/38。Mathematica有一个有趣的特性，即如果你以分数形式输入一个值，该值将保持这种形式，而不是转换为浮点形式，从而在后续计算中保持其准确性。正如你在笔记本中看到的，对于这个案例，解保持在这个初始值。根据我们之前关于不动点和稳定性的讨论，我们会说在这种情况下，28/38是一个不动点，但一个不稳定的不动点，因为即使从N的这个值有非常轻微的偏差，也会导致解偏离不动点。希望你看作是这个混沌系统对初始条件高度敏感的另一种结果。