chapter23.ipynb

{
 "cells": [
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "# 第二十三讲：微分方程和$e^{At}$\n",
    "\n",
    "## 微分方程$\\frac{\\mathrm{d}u}{\\mathrm{d}t}=Au$\n",
    "\n",
    "本讲主要讲解解一阶方程（first-order system）一阶倒数（first derivative）常系数（constant coefficient）线性方程，上一讲介绍了如何计算矩阵的幂，本讲将进一步涉及矩阵的指数形式。我们通过解一个例子来详细介绍计算方法。\n",
    "\n",
    "有方程组$\\begin{cases}\\frac{\\mathrm{d}u_1}{\\mathrm{d}t}&=-u_1+2u_2\\\\\\frac{\\mathrm{d}u_2}{\\mathrm{d}t}&=u_1-2u_2\\end{cases}$，则系数矩阵是$A=\\begin{bmatrix}-1&2\\\\1&-2\\end{bmatrix}$，设初始条件为在$0$时刻$u(0)=\\begin{bmatrix}u_1\\\\u_2\\end{bmatrix}=\\begin{bmatrix}1\\\\0\\end{bmatrix}$。\n",
    "\n",
    "* 这个初始条件的意义可以看做在开始时一切都在$u_1$中，但随着时间的推移，将有$\\frac{\\mathrm{d}u_2}{\\mathrm{d}t}>0$，因为$u_1$项初始为正，$u_1$中的事物会流向$u_2$。随着时间的发展我们可以追踪流动的变化。\n",
    "\n",
    "* 根据上一讲所学的知识，我们知道第一步需要找到特征值与特征向量。$A=\\begin{bmatrix}-1&2\\\\1&-2\\end{bmatrix}$，很明显这是一个奇异矩阵，所以第一个特征值是$\\lambda_1=0$，另一个特征向量可以从迹得到$tr(A)=-3$。当然我们也可以用一般方法计算$\\left|A-\\lambda I\\right|=\\begin{vmatrix}-1-\\lambda&2\\\\1&-2-\\lambda\\end{vmatrix}=\\lambda^2+3\\lambda=0$。\n",
    "\n",
    "    （教授提前剧透，特征值$\\lambda_2=-3$将会逐渐消失，因为答案中将会有一项为$e^{-3t}$，该项会随着时间的推移趋近于$0$。答案的另一部分将有一项为$e^{0t}$，该项是一个常数，其值为$1$，并不随时间而改变。通常含有$0$特征值的矩阵会随着时间的推移达到稳态。）\n",
    "\n",
    "* 求特征向量，$\\lambda_1=0$时，即求$A$的零空间，很明显$x_1=\\begin{bmatrix}2\\\\1\\end{bmatrix}$；$\\lambda_2=-3$时，求$A+3I$的零空间，$\\begin{bmatrix}2&2\\\\1&1\\end{bmatrix}$的零空间为$x_2=\\begin{bmatrix}1\\\\-1\\end{bmatrix}$。\n",
    "\n",
    "* 则方程组的通解为：$u(t)=c_1e^{\\lambda_1t}x_1+c_2e^{\\lambda_2t}x_2$，通解的前后两部分都是该方程组的纯解，即方程组的通解就是两个与特征值、特征向量相关的纯解的线性组合。我们来验证一下，比如取$u=e^{\\lambda_1t}x_1$带入$\\frac{\\mathrm{d}u}{\\mathrm{d}t}=Au$，对时间求导得到$\\lambda_1e^{\\lambda_1t}x_1=Ae^{\\lambda_1t}x_1$，化简得$\\lambda_1x_1=Ax_1$。\n",
    "\n",
    "    对比上一讲，解$u_{k+1}=Au_k$时得到$u_k=c_1\\lambda^kx_1+c_2\\lambda^kx_2$，而解$\\frac{\\mathrm{d}u}{\\mathrm{d}t}=Au$我们得到$u(t)=c_1e^{\\lambda_1t}x_1+c_2e^{\\lambda_2t}x_2$。\n",
    "    \n",
    "* 继续求$c_1,c_2$，$u(t)=c_1\\cdot 1\\cdot\\begin{bmatrix}2\\\\1\\end{bmatrix}+c_2\\cdot e^{-3t}\\cdot\\begin{bmatrix}1\\\\-1\\end{bmatrix}$，已知$t=0$时，$\\begin{bmatrix}1\\\\0\\end{bmatrix}=c_1\\begin{bmatrix}2\\\\1\\end{bmatrix}+c_2\\begin{bmatrix}1\\\\-1\\end{bmatrix}$（$Sc=u(0)$），所以$c_1=\\frac{1}{3}, c_2=\\frac{1}{3}$。\n",
    "\n",
    "* 于是我们写出最终结果，$u(t)=\\frac{1}{3}\\begin{bmatrix}2\\\\1\\end{bmatrix}+\\frac{1}{3}e^{-3t}\\begin{bmatrix}1\\\\-1\\end{bmatrix}$。\n",
    "\n",
    "稳定性：这个流动过程从$u(0)=\\begin{bmatrix}1\\\\0\\end{bmatrix}$开始，初始值$1$的一部分流入初始值$0$中，经过无限的时间最终达到稳态$u(\\infty)=\\begin{bmatrix}\\frac{2}{3}\\\\\\frac{1}{3}\\end{bmatrix}$。所以，要使得$u(t)\\to 0$，则需要负的特征值。但如果特征值为复数呢？如$\\lambda=-3+6i$，我们来计算$\\left|e^{(-3+6i)t}\\right|$，其中的$\\left|e^{6it}\\right|$部分为$\\left|\\cos 6t+i\\sin 6t\\right|=1$，因为这部分的模为$\\cos^2\\alpha+\\sin^2\\alpha=1$，这个虚部就在单位圆上转悠。所以只有实数部分才是重要的。所以我们可以把前面的结论改为**需要实部为负数的特征值**。实部会决定最终结果趋近于$0$或$\\infty$，虚部不过是一些小杂音。\n",
    "\n",
    "收敛态：需要其中一个特征值实部为$0$，而其他特征值的实部皆小于$0$。\n",
    "\n",
    "发散态：如果某个特征值实部大于$0$。上面的例子中，如果将$A$变为$-A$，特征值也会变号，结果发散。\n",
    "\n",
    "再进一步，我们想知道如何从直接判断任意二阶矩阵的特征值是否均小于零。对于二阶矩阵$A=\\begin{bmatrix}a&b\\\\c&d\\end{bmatrix}$，矩阵的迹为$a+d=\\lambda_1+\\lambda_2$，如果矩阵稳定，则迹应为负数。但是这个条件还不够，有反例迹小于$0$依然发散：$\\begin{bmatrix}-2&0\\\\0&1\\end{bmatrix}$，迹为$-1$但是仍然发散。还需要加上一个条件，因为$\\det A=\\lambda_1\\cdot\\lambda_2$，所以还需要行列式为正数。\n",
    "\n",
    "总结：原方程组有两个相互耦合的未知函数，$u_1, u_2$相互耦合，而特征值和特征向量的作则就是解耦，也就是对角化（diagonalize）。回到原方程组$\\frac{\\mathrm{d}u}{\\mathrm{d}t}=Au$，将$u$表示为特征向量的线性组合$u=Sv$，代入原方程有$S\\frac{\\mathrm{d}v}{\\mathrm{d}t}=ASv$，两边同乘以$S^{-1}$得$\\frac{\\mathrm{d}v}{\\mathrm{d}t}=S^{-1}ASv=\\Lambda v$。以特征向量为基，将$u$表示为$Sv$，得到关于$v$的对角化方程组，新方程组不存在耦合，此时$\\begin{cases}\\frac{\\mathrm{d}v_1}{\\mathrm{d}t}&=\\lambda_1v_1\\\\\\frac{\\mathrm{d}v_2}{\\mathrm{d}t}&=\\lambda_2v_2\\\\\\vdots&\\vdots\\\\\\frac{\\mathrm{d}v_n}{\\mathrm{d}t}&=\\lambda_nv_n\\end{cases}$，这是一个各未知函数间没有联系的方程组，它们的解的一般形式为$v(t)=e^{\\Lambda t}v(0)$，则原方程组的解的一般形式为$u(t)=e^{At}u(0)=Se^{\\Lambda t}S^{-1}u(0)$。这里引入了指数部分为矩阵的形式。\n",
    "\n",
    "## 指数矩阵$e^{At}$\n",
    "\n",
    "在上面的结论中，我们见到了$e^{At}$。这种指数部分带有矩阵的情况称为指数矩阵（exponential matrix）。\n",
    "\n",
    "理解指数矩阵的关键在于，将指数形式展开称为幂基数形式，就像$e^x=1+\\frac{x^2}{2}+\\frac{x^3}{6}+\\cdots$一样，将$e^{At}$展开成幂级数的形式为：\n",
    "\n",
    "$$e^{At}=I+At+\\frac{(At)^2}{2}+\\frac{(At)^3}{6}+\\cdots+\\frac{(At)^n}{n!}+\\cdots$$\n",
    "\n",
    "再说些题外话，有两个极具美感的泰勒级数：$e^x=\\sum \\frac{x^n}{n!}$与$\\frac{1}{1-x}=\\sum x^n$，如果把第二个泰勒级数写成指数矩阵形式，有$(I-At)^{-1}=I+At+(At)^2+(At)^3+\\cdots$，这个式子在$t$非常小的时候，后面的高次项近似等于零，所以可以用来近似$I-At$的逆矩阵，通常近似为$I+At$，当然也可以再加几项。第一个级数对我们而言比第二个级数好，因为第一个级数总会收敛于某个值，所以$e^x$总会有意义，而第二个级数需要$A$特征值的绝对值小于$1$（因为涉及矩阵的幂运算）。我们看到这些泰勒级数的公式对矩阵同样适用。\n",
    "\n",
    "回到正题，我们需要证明$Se^{\\Lambda t}S^{-1}=e^{At}$，继续使用泰勒级数：\n",
    "\n",
    "$$\n",
    "e^{At}=I+At+\\frac{(At)^2}{2}+\\frac{(At)^3}{6}+\\cdots+\\frac{(At)^n}{n!}+\\cdots\\\\\n",
    "e^{At}=SS^{-1}+S\\Lambda S^{-1}t+\\frac{S\\Lambda^2S^{-1}}{2}t^2+\\frac{S\\Lambda^3S^{-1}}{6}t^3+\\cdots+\\frac{S\\Lambda^nS^{-1}}{n!}t^n+\\cdots\\\\\n",
    "e^{At}=S\\left(I+\\Lambda t+\\frac{\\Lambda^2t^2}{2}+\\frac{\\Lambda^3t^3}{3}+\\cdots+\\frac{\\Lambda^nt^n}{n}+\\cdots\\right)S^{-1}\\\\\n",
    "e^{At}=Se^{\\Lambda t}S^{-1}\n",
    "$$\n",
    "\n",
    "需要注意的是，$e^{At}$的泰勒级数展开是恒成立的，但我们推出的版本却需要**矩阵可对角化**这个前提条件。\n",
    "\n",
    "最后，我们来看看什么是$e^{\\Lambda t}$，我们将$e^{At}$变为对角矩阵就是因为对角矩阵简单、没有耦合，$e^{\\Lambda t}=\\begin{bmatrix}e^{\\lambda_1t}&0&\\cdots&0\\\\0&e^{\\lambda_2t}&\\cdots&0\\\\\\vdots&\\vdots&\\ddots&\\vdots\\\\0&0&\\cdots&e^{\\lambda_nt}\\end{bmatrix}$。\n",
    "\n",
    "有了$u(t)=Se^{\\Lambda t}S^{-1}u(0)$，再来看矩阵的稳定性可知，所有特征值的实部均为负数时矩阵收敛，此时对角线上的指数收敛为$0$。如果我们画出复平面，则要使微分方程存在稳定解，则特征值存在于复平面的左侧（即实部为负）；要使矩阵的幂收敛于$0$，则特征值存在于单位圆内部（即模小于$1$），这是幂稳定区域。（上一讲的差分方程需要计算矩阵的幂。）\n",
    "\n",
    "同差分方程一样，我们来看二阶情况如何计算，有$y''+by'+k=0$。我们也模仿差分方程的情形，构造方程组$\\begin{cases}y''&=-by'-ky\\\\y'&=y'\\end{cases}$，写成矩阵形式有$\\begin{bmatrix}y''\\\\y'\\end{bmatrix}=\\begin{bmatrix}-b&-k\\\\1&0\\end{bmatrix}\\begin{bmatrix}y'\\\\y\\end{bmatrix}$，令$u'=\\begin{bmatrix}y''\\\\y'\\end{bmatrix}, \\ u=\\begin{bmatrix}y'\\\\y\\end{bmatrix}$。\n",
    "\n",
    "继续推广，对于$5$阶微分方程$y'''''+by''''+cy'''+dy''+ey'+f=0$，则可以写作$\\begin{bmatrix}y'''''\\\\y''''\\\\y'''\\\\y''\\\\y'\\end{bmatrix}=\\begin{bmatrix}-b&-c&-d&-e&-f\\\\1&0&0&0&0\\\\0&1&0&0&0\\\\0&0&1&0&0\\\\0&0&0&1&0\\end{bmatrix}\\begin{bmatrix}y''''\\\\y'''\\\\y''\\\\y'\\\\y\\end{bmatrix}$，这样我们就把一个五阶微分方程化为$5\\times 5$一阶方程组了，然后就是求特征值、特征向量了步骤了。"
   ]
  }
 ],
 "metadata": {
  "kernelspec": {
   "display_name": "Python 3",
   "language": "python",
   "name": "python3"
  },
  "language_info": {
   "codemirror_mode": {
    "name": "ipython",
    "version": 3
   },
   "file_extension": ".py",
   "mimetype": "text/x-python",
   "name": "python",
   "nbconvert_exporter": "python",
   "pygments_lexer": "ipython3",
   "version": "3.5.1"
  }
 },
 "nbformat": 4,
 "nbformat_minor": 0
}