Python数学建模极简入门（三）简单动力系统

这一期讲动力系统

在准备这一期教程的过程中我找到了一个挺有意思的课件，推荐给大家:
动力系统建模-唐云（清华大学教授）
里面讲到了一个很吸引目光的名词：
蝴蝶效应
大家应该都知道蝴蝶效应是个什么玩意儿，但是蝴蝶效应是怎么产生的呢？
这次先不讲蝴蝶效应，我们先就这个动力系统来讲一讲,算是个引子：

---

什么叫动力系统呢，先说最简单的动力系统，讲起来也是个很熟的公式：

a_n+1 = r × a_n + b(其中r≠1）

这不还是递推公式吗？

和上期的差分方程那个递推公式不同的是，这个动力系统有个被称为平衡点的概念。

废话少说，先看一个最简单的实例：

---

实例1--血液中药物水平试验

有一个医生给病人开药，什么药我不记得了，也不重要，就叫A吧，要求病人血液中的钙离子浓度稳定在一个有效的值，求如何开药？
前提：
A在病人体内原本就有。
不服药的情况下，每一周的A浓度水平是前一周的一半。
列式：
a_n+1 = 0.5 × a_n + b
这里为了方便画图我们先假设b = 1，先画个图看看。
那么式子变成：a_n+1 = 0.5 × a_n + 1
接下来我们分别取
a₀ = 1
a₀ = 2
a₀ = 3
画图

python代码如下：
<code>
先画a₀ = 1：
import matplotlib.pyplot as plt
def drawCurve(Reagentperday,CurrentLevel):
RegentLevel = []
for i in range(15):
RegentLevel.append(CurrentLevel)
CurrentLevel = CurrentLevel * 0.5 + Reagentperday
Time = [i for i in range(15)]
print(Time,RegentLevel)
plt.plot(Time, RegentLevel)
for i in range(1,4):
drawCurve(1,i)
</code>

画出的图形是这样的：

药物浓度水平随着时间变化图像

这个时候就出现了一个很神奇的事情：
无论初始水平是多少，采用这种服药方式的话，最终药物水平都会收敛于2.0这个值，即所谓平衡点

是这样的吗，我们再多画几条：

药物初始浓度为1~9时水平变化

看来确实有这样的规律。

那么变一下每天服药的数量之后呢，我们将每天的服药量作为变量来画一下：

固定初始值变动每日服药量时

似乎对应每个服药方案都会有一个平衡点。而且我们可以看到，对于每日服药量±1的变动会引起平衡点±2的变化。
是不是呢，我们再画一下不同服药方案在不同初始浓度下的图像，得到下面这个酷炫的图形：

每日服药量与初始值都发生变化时

这个地方可能是我显卡有问题，截出来的图不是很清晰

究竟是为什么呢？

---

下一期我们讲讲为什么，主要是依据文章开头贴出来的动力系统建模-唐云（清华大学教授）的ppt来用Python代码讲解。

就这样结束了吗？？
熟悉高数的人肯定大呼上当。
因为
a_n+1 = r × a_n + b
这个式子可以写成：
a_n+1 + b/(r-1) = r × a_n + b + b/(r-1)
a_n+1 + b/(r-1) = r × [a_n +b/(r-1)]

可以判定[a_n + b/(r-1)]是一个公比为r的等比数列，当r<1，且首项≠0时，必有a收敛于[a₀+b/(r-1)]/(1-r)

动力系统就这么简单吗？

这只是一次的动力系统，那么二次的呢……