微分方程数值解

Introduction

本学期还选修了数学建模(Mathematical Modeling),也算是给自己的大学数学生涯画上一个句号吧。所以也会更新一些数学建模的相关的知识，其实自己感觉这些基本的数学知识还是蛮重要的，虽然深度学习之中使用的高深的数学的地方似乎并不多，但是不管怎么说，搞学术研究不懂数学还是说不过去，不知道自己以后有没有精力写一写关于基本数学的东西呀。看后续的精力和反馈吧(现在还是自娱自乐阶段。)

这个网站的数学建模的东西应该不会涉及到特别深，可能很多地方会点到为止，因为自己精力实在是有限(毕竟还要兼顾吃吃喝喝打游戏拍照看演出看电影看电视看小说看漫画)，不过大家如果有什么问题想要和我深入探讨的话，我也会好好再去学习一下哒。

既然是数学建模，就从实际的问题开始吧，有一个模型叫做是食饵捕食者模型(Prey Predator Model)，大概的意思就是在一个系统之中，有两个种群，分别是甲和乙，其中甲依靠丰富的自然资源生长，而乙则依靠捕食甲而生存。

模型有如下的假设，甲和乙在$t$时刻的数量分别是$x(t)$以及$y(t)$，其中甲单独生存的时候的自然的增长率是$r$，但是由于乙的存在，甲的增长率也会相应的减少，这个减少的幅度与乙的数量是成正比的，假设这样的一个正比是$a$，$a$所反映的就是捕食者的捕食的能力的强弱。而对于捕食者乙而言，在无食饵的时候自然的死亡率为$d$,但是由于食饵的存在，捕食者的死亡率降低，并且使得其数量的介绍的程度与食饵的数量成正比的关系，假设这个比率是$b$,也就是大概可以说明食饵的营养的程度。我们用数学的公式进行翻译一下哈

下面的就是食饵的种群的数量的变化(碎碎念一个，在Web端写LaTeX的话，改微分号还是挺麻烦的事情，自己以后就偷个懒了大家看到这里如果发现$d$是斜的，一定要记住真正的d是正的。)

$$\frac{d x}{d t} = x(t)r’ = x(t)(r – a y(t)) = rx -axy$$

而对于捕食者的种群满足

$$\frac{d y}{d t} = -y(t)d’ = -y(t)(d – b x(t)) = -dy +bxy$$

所以根据上面的情况就可以列出微分方程组

$$\begin{cases}\frac{d x}{d t} =  rx -axy \\    \frac{d y}{d t} = -y(t)d’ = -dy +bxy   \end{cases}$$

两个方程相除，可以得到一个方程

$$\frac{d x}{d y}  = \frac{rx-axy}{-dy + bxy}$$

这个的解称为 相轨线 ，名字骚里骚气的，其实就是把解空间的点给描出来。大家如果对于L1正则和L2正则有所了解的话，L1正则就是把解空间给丢进一个方形空间了。不过这里是后话了，大家如果没有了解过自然在这里的理解也是没啥问题。

数值解

对于这样的一个微分方程的解还是挺复杂的问题，毕竟很多时候我们根本没有办法找到精确的解，而只能通过模拟而获得一个近似的解。怎么说呢？打个比方，大家可能都用那种很透明的纸放在书上去临摹一些什么东西，这就是求数值解的Intuition，我们无法得到书上的一模一样的内容，只能是通过自己的临摹得到一个类似的样子。而临摹的像不像，在数学里面也有一个单独的名词，叫做精度。好了例子打完了，还是要面对数学的

假设有一个长这样的微分方程

$$\begin{cases}y’ =  f(x,y ) \quad  a\leq x \leq b \\    y(a) = y_0 \end{cases}$$

假设这样的一个微分方程的解$y(x)$存在并且唯一，其解析解很难求的情况下，我们转而利用一系列的离散点$x_0 < x_1 <\cdots <x_n < \cdots$求它的近似值，然后将近似值作为实际值就行了。对于$y’ =  f(x,y )$的两边取积分

$$\int_{x_k}^{x_{k+1}} y’ d x= \int_{x_k}^{x_{k+1}}  f(x,y ) dx$$

可以得到

$$y(x_{k+1})  = y(x_k) + \int_{x_k}^{x_{k+1}}f(x,y) dx$$

可以发现，在获得了$y(x_k)$的情况下，计算$y(x_{k+1})$需要的仅仅是计算$\int_{x_k}^{x_{k+1}}f(x,y) dx$，那么这个鬼东西怎么计算呢？这个鬼东西是非常难计算的，那么我们就要弄一些门道，什么门道呢？大家都学过微积分吧，微积分的基本的思想之一就是以直代曲，将面积分成一个一个的小块的区域然后累加起来就得到了微积分的公式，那么这个呢，就不分成一个一个的小块的区域了，而是直接用一个大矩形表示。下面的图我就暂时不贴了，如果有人看的话，在评论区说一下我再找来贴上去。

矩形公式如下所示

$$\int_{x_k}^{x_{k+1}}f(x,y) dx    \simeq h f(x_{k},y_{k})\simeq  h f(x_{k+1},y_{k+1}) $$

梯形公式是对于矩形公式的一种改进

$$\int_{x_k}^{x_{k+1}}f(x,y) dx    \simeq \frac{h}{2} [f(x_{k},y_{k})+ f(x_{k+1},y_{k+1})] $$

而辛普森公式，不是辛普森一家的那个辛普森哈

$$\begin{split}&\int_{x_k}^{x_{k+1}}f(x,y) dx    \\ & \simeq  \frac{h}{6} [f(x_{k},y_{k})+4f(x_{k+\frac{h}{2}},y_{k+\frac{h}{2}})+ f(x_{k+1},y_{k+1})] \end{split}$$

计算出来这个恶心的积分之后，就是用

$$y(x_{k+1})  = y(x_k) + \text{恶心的积分}$$

来迭代下个$y_{k+1}$了

如果是矩形公式的话那么就是

$$y(x_{k+1})  \simeq y(x_k) + h f(x_{k},y_{k})$$

有一个好听的名字，叫做是向前Euler求解公式，没错，又是欧拉这个勤奋的家伙(这个梗大家知道可以在评论区打出来)。

如果是

$$y(x_{k+1})  \simeq y(x_k) + h f(x_{k+1},y_{k+1})$$

就是向后欧拉公式。

有什么不同呢？

这个就是不同之处，其中红色的点是前向欧拉公式所计算出的，而绿色的线则是后向的计算的过程，蓝色的线是正常的函数。可以发现前向欧拉公式会倾向于包住函数本身，这个和函数的凹凸性有没有关系呢？留作思考题，大家可以在评论区扣出来。

而梯形公式就也是很简单的一个替换的过程

$$y(x_{k+1})  \simeq y(x_k) + \frac{h}{2} [f(x_{k},y_{k})+ f(x_{k+1},y_{k+1})] $$

但是不知道他家有没有注意到一个问题就是，不管是向后的欧拉公式，还是矩形公式，我们的$y_{k+1}$都是未知的，如果是未知的我们要怎么将这个值带入呢？很简单哈，首先使用向前欧拉公式计算出一个伪$y_{k+1}$然后用这个伪数带入这两个公式之中。有人说这不是脱裤子放屁吗？我都有了一个数字了，干嘛还要费尽去做这个事情。事实上这样会使得精度大大的提升，也就是说这个屁确实放的有道理。矩形公式的精度是1阶，而梯形就到2阶了。之后的辛普森花式脱裤子放屁还能到4阶。

而数学家除了探索一些奇奇怪怪的定理之外，还喜欢做的一个事情就是总结，把几个有着相似的性质的东西都用一套理论总结起来。比如说我们的龙哥Runge和库哥Kutta，他们两个把这些东西总结了一下，就不叫欧拉公式和辛普森公式了，全部都是叫Runge-Kutta方法，大概的就是：你这个公式很好，可惜下一秒就是我的了。

梯形公式又叫做是二阶Rugge-Kutta法，而辛普森公式叫做三阶Rugge-Kutta方法。当然大胆假设还有四阶的龙哥方法，阶数越高代表越精确。四阶的就是再多插一个点进去，用四个点进行拟合，五阶六阶都是换汤不换药，我这也不浪费大家的时间了。

PP Model

好了，兜兜转转这么久，我们又回到了食饵捕食者模型，我们的微分方程组的形式是

$$\begin{cases}\frac{d x}{d t} =  rx -axy \\    \frac{d y}{d t} = -y(t)d’ = -dy +bxy   \end{cases}$$

其中的$r,d,a,b$都是常数，而两个种群的初始的数量是$x_0$以及$y_0$,我们的使用的软件是MATLAB，为什么不用Python呢？因为不得不说MATLAB的确是有很多很方便的地方，建议还是得学。

MATLAB处理矩阵有天生的优势，里面的所有的运算都是基于矩阵所进行的，所以我们需要将微分方程组转换成为矩阵的形式

$$\begin{pmatrix} x_1′  \\  x_2′ \end{pmatrix} = \begin{pmatrix} rx_1 – ax_1x_2  \\  -dx_2 + bx_1x_2 \end{pmatrix}$$

MATLAB的代码的文件由于自己在网站上暂时没有什么用代码的Plugin，而且WordPress的网站总是没有连接，所以暂时就不贴漂漂亮亮的代码了，就贴个图先吧，如果有想要代码的可以在下面给我留言，我再给你发。

自己突然挺想学点PDE的，快来阻止我！