线性ODE（三）——初探高阶线性方程

一些让步

直接解释方程计算办法，看似很实在，但是反而会让问题变得更复杂，因为对结构很不清楚，只剩下一个空穴来风的公式。但是，接下来我们不得不假设我们能很轻松地算出一个方程的解(就当是超能力吧！)——因为有的方程本来就没有什么普遍的解决办法。这篇博客的目的在于，让读者认识到，就算没有具体的解决办法，也能宏观上对方程进行一些合适的处理。这篇博客中也不是解释如何解方程，而是关于如何利用矩阵和一些常见的数学理论，分析未知的问题。解方程的具体办法会在下一篇博客。

解结构的分析

从一阶一次线性方程中得到的启示

最开始我们分析的是这种方程

$$y'+a(x)y=b(x)$$

如果我们有一个解$y_1(x)$使得

$$y_1'(x)+a(x)y_1(x)=0$$

而且有一个特殊的解$y_2(x)$使得

$$y_2'(x)+a(x)y_2'(x)=b(x)$$

那么全体解就可以表示成$Cy_1(x)+y_2(x)$，有

$$(Cy_1(x)+y_2(x))'+a(x)(Cy_1(x)+y_2(x))=b(x)$$

所以如果是$n$个未知数呢？会不会是得到$n$个解，然后再有一个特解？但是一定得是$n$个吗？会不会是更多，或者更少个？如果是$n$个，会不会和一般研究的$n$维欧几里得空间建立联系？比如处理一个$2$元方程，是不是可以投影到二维平面？幸运的是，我们可以很轻松地把欧几里得空间和线性方程的解建立起联系。

线性方程与欧几里得空间的联系

这篇博客（和接下来的博客）中，加粗的数学符号都代表列向量或者矩阵。如果觉得不适应，可以自己写成向量和矩阵方便理解。

解的存在性和唯一性

线性方程组

$$\frac{d\mathbf{y}}{dx}=\mathbf{A}(x)\mathbf{y}+\mathbf{f}(x)$$

的满足条件

$$\mathbf{y}(x_0)=\mathbf{y_0}$$

的解是存在且唯一的(这里要求$\mathbf{f}$和$\mathbf{A}$中的函数都是连续的)。这个问题的证明是一个比较麻烦的事情，涉及到Banach不动点问题，在这里就不做证明了。但是也是很好理解的：这里的$\mathbf{A}$和$\mathbf{f}$都是连续的，都不是“病态”的，所以这个方程不会有“病态的”形式，解是很自然的。

和欧氏空间的一一对应

我们宣称，线性方程组

$$\frac{d\mathbf{y}}{dx}=\mathbf{A}(x)\mathbf{y}$$

的解和$n$维欧氏空间是一一对应的。也就是说，对每个$\mathbf{y}(x)$，都对应$n$维欧氏空间的唯一一个向量；反之，对$n$维欧氏空间的每一个向量，都对应着一个方程组的唯一一个解。

这其实就是所谓解的存在性和唯一性的直接应用。在定义域固定一个$x_0$，然后我们选择$\mathbb{R}^n$的自然基，也就是

$$\{\mathbf{e_1},\mathbf{e_2},\cdots,\mathbf{e_n}\}$$

其中$\mathbf{e_1}=(1,0,\cdots,0)$, $\mathbf{e_n}=(0,\cdots,0,1)$。

那么满足

$$\begin{aligned} \begin{cases} \mathbf{y_i}'(x)&=\mathbf{A}(x)\mathbf{y_i} \\ \mathbf{y_i}(x_0)&=\mathbf{e_i} \end{cases} \end{aligned}$$

的解就是存在且唯一的。如果这个对应，所以我们也选择了

$$\{\mathbf{y_1}(x),\mathbf{y_2}(x),\cdots,\mathbf{y_n}(x)\}$$

这样一组“基”(这里我们假设我们已经有一下子求出解的超能力)，那么能不能建立一个$\mathbb{R}^n$的平行版？答案是可以的。对于$\mathbb{R}^n$ ，我们肯定有

$$C_1\mathbf{e_1}+C_2\mathbf{e_2} \in \mathbb{R}^n$$

我们设这一节中方程的全体解集合为$\mathfrak{R}^n$，那么我们已经有了$\mathbf{y_i} \in \mathfrak{R}^n$，我们希望

$$C_1\mathbf{y_1}+C_2\mathbf{y_2} \in \mathfrak{R}^n$$

这其实利用求导的性质和矩阵的性质，就可以解决。因为

$$\frac{d(C_1\mathbf{y_1}+C_2\mathbf{y_2})}{dx}=C_1\frac{d\mathbf{y_1}}{dx}+C_2\frac{d\mathbf{y_2}}{dx}=C_1\mathbf{A}(x)\mathbf{y_1}+C_2\mathbf{A}(x)\mathbf{y_2}=\mathbf{A}(x)(C_1\mathbf{y_1}+C_2\mathbf{y_2})$$

而且不难验证，$C_1\mathbb{y_1}+C_2\mathbf{y_2}$和$C_1\mathbf{e_1}+C_2\mathbf{e_2}$也是一一对应的。所以，我们可以很自然地构建一个和$\mathbb{R}^n$一一对应的平行的空间，尽管这个空间我们可能没法可视化。

我们知道，$\{\mathbf{e_1},\mathbf{e_2},\cdots,\mathbf{e_n}\}$这组向量，能很自然地用加法和数乘表示出$\mathbb{R}^n$中任意一个向量，而我们建立的等价关系也已经保证，$\{\mathbf{y_1}(x),\mathbf{y_2}(x),\cdots,\mathbf{y_n}(x)\}$也能用加法和数乘，表示出$\mathfrak{R}^n$中的任意一个解。换句话说，本节的方程的全体解$\mathfrak{R}^n$，其中的任意一个解，都可以表示成

$$\mathbf{y_0}=\sum_{k=1}^{n}C_k\mathbf{y_k}$$

其中$C_i$为常数。这里其实是一个同构(isomorphism)关系，如果了解这个概念，可以很轻松的写出严格的证明，但在这里就不深入了。

注意，我们这里已经强迫$\mathbf{f}=\mathbf{0}$. 这样一个$\mathbf{y_0}$，就能保证$\mathbf{y’=Ay+f}$的普遍性。如果我们已经得到了一个$\mathbf{y’=Ay+f}$的特殊解$\mathbf{y_1}$，保证了准确性，那么，$\mathbf{y_0+y_1}$，总共$n+1$个函数，就能表示出所有解了。$\mathbf{f}=0$时，我们的$\mathbf{y_1}=0$也是很和谐的。我们在解决多元微分方程时，处理的问题和一元一次情况下，是很类似的。

尝试解决高阶方程

我们现在还有那个超能力，能一下子解决$\frac{d\mathbf{y}}{dx}=\mathbf{A}(x)\mathbf{y}$这样的方程，能求出一组基对应的解。如果我们令

$$\mathbf{A}= \begin{pmatrix} 0&1&0&\cdots&0 \\ 0&0&1&\cdots&0 \\ \vdots&\vdots&\vdots&\ddots&\vdots\\ 0&0&0&\cdots&1\\ -a_n(x)&-a_{n-1}(x)&-a_{n-2}(x)&\cdots&-a_1(x) \end{pmatrix}$$

那么我们就可以得到一个方程

$$y^{(n)}+a_1(x)y^{(n-1)}+\cdots+a_n(x)y=0$$

所以我们已经得到了$n$个不同的解，不妨设这组解是$y_1,y_2,\cdots,y_n$。分别对这些$y_i$求导$n-1$次，就能得到上一节中的$n$维列向量。

那么如果是

$$y^{(n)}+a_1(x)y^{(n-1)}+\cdots+a_n(x)y=f(x)$$

这种方程，该怎么求解呢？我们还是通过类似于之前一阶线性方程的办法，通过求导运算，间接求出来。假设有一个函数$y$满足这个方程。

我们不妨设

$$\mathbf{y}=\begin{bmatrix}y_1\\ \vdots \\ y_n\end{bmatrix}$$

以及

$$\mathbf{c}=(c_1,c_2,\cdots,c_n)$$

其中满足$y_p=\mathbf{c \cdot y}$(利用行向量、列向量的乘法规则)。所以我们接下来的工作就是求出这样一个$\mathbf{c}$。我们会尝试用线性代数的语言进行描述，可以看到，表述变得非常简单。在这里，对向量求导、求积分的意思是

$$\mathbf{y}'=\begin{bmatrix}y_1'\\ \vdots \\ y_n'\end{bmatrix} \\ \int\mathbf{y}dx=\begin{bmatrix}\int y_1dx\\ \vdots \\\int y_ndx\end{bmatrix}$$

但是直接进行下去也似乎漫无边际。我们希望，能把这个$c$给稳定下来，所以我们会尝试加一些限制条件。这样的限制应该是$n-1$条。这就好比，如果要在平面中确定，和一条直线垂直的直线，那么只需要一个向量就行了。如果要在空间中确定，和一个平面垂直的一条直线，确定两个向量就行了——所谓两条相交直线确定一个平面。

我们来分析一下怎么确定会更好。对$y$求导，我们有

$$y_p'=\mathbf{c'\cdot y+c\cdot y'}$$

如果我们令$c’y=0$，那么就有

$$y''_p=(\mathbf{c \cdot y'})'=\mathbf{c' \cdot y'}+\mathbf{c \cdot y''}$$

如果再令$\mathbf{c’ \cdot y’}=0$，继续下去，甚至是$\mathbf{c’ \cdot y^{(k)}}=0$，其中$0 \leq k \leq n-2$，那么对应就有$y_p^{(k)}=cy^{(k)}$。而且接下来求导就变得比较简单。但是，如果我们令$\mathbf{c \cdot y^{(p)}}=0$，就很麻烦了，因为接下来多了像是$c’,c’’,\cdots$这种未知量。这会让等式越来越长，没法表达。所以，我们选择限制条件

$$\mathbf{c' \cdot y^{(k)}}=0\quad 0 \leq k \leq n-2$$

已经选择了$n-1$个限制条件，那么能不能凑个整，得到第$n$个式子？那样的话，有$n$个未知数，有$n$个等式，像是再搞矩阵了。我们还有一点根本没用上：最开始的方程。对于$0 \leq k \leq n-1$，我们已经已经有(可以手动验算一下)：

$$y_{p}^{(k)}=\mathbf{c \cdot y^{(k)}}$$

所以就有

$$y_p^{(n)}=\mathbf{c' \cdot y^{(n-1)}}+\mathbf{c \cdot y^{(n)}}$$

这个时候，回到最开始的方程，代入这$n$个等式，有

$$\begin{aligned} y_p^{(n)}+a_1y_p^{(n-1)}+\cdots+a_ny_p&=f\\ \mathbf{c' \cdot y^{(n-1)}}+\mathbf{c \cdot y^{(n)}}+a_1\mathbf{c\cdot y^{(n-1)}}+\cdots+a_n\mathbf{c \cdot y}&=f\\ \mathbf{c' \cdot y^{(n-1)}}+\mathbf{c \cdot}\left(\mathbf{y^{(n)}}+a_1\mathbf{y^{(n-1)}}+\cdots+a_n\mathbf{y}\right)&=f \\ \mathbf{c' \cdot y^{(n-1)}}&=f \end{aligned}$$

最后一个等号成立是因为，在$\mathbb{y}$中，每个分量$y_k$都满足$y_k+a_1y_k^{(n-1)}+\cdot+a_ny_k=0$。所以这些向量加起来是一个$0$向量。

所以我们有了关于$\mathbf{c’}$的$n$个等式，也就是说

$$\begin{cases} c_1'y_1+\cdots+c_n'y_n&=&0 \\ c_1'y_1'+\cdots+c_n'y_n'&=&0 \\ \cdots\cdots\cdots\cdots\\ c_1'y_1^{(n-2)}+\cdots+c_n'y_{n}^{(n-2)}&=&0 \\ c_1'y_1^{(n-1)}+\cdots+c_n'y_{n}^{(n-1)}&=&f \end{cases}$$

如果我们写成矩阵的形式，就有

$$\begin{bmatrix} y_1 & y_2 & \cdots&y_{n} \\ y_1' & y_2' & \cdots&y_{n}' \\ \vdots&\vdots&\ddots&\vdots \\ y_1^{(n-2)} & y_2^{(n-2)} & \cdots&y_{n}^{(n-2)} \\ y_1^{(n-1)} & y_2^{(n-1)} & \cdots&y_{n}^{(n-1)} \\ \end{bmatrix} \begin{bmatrix} c_1'\\ c_2'\\ \vdots \\ c_{n-1}' \\ c_n' \end{bmatrix}= \begin{bmatrix} 0 \\ 0 \\ \vdots \\ 0 \\ f \end{bmatrix}$$

如果你学习了Cramer法则，那么怎么解决这个问题已经浮现在眼前了。没错，我们可以直接利用Cramer法则，求出$\mathbf{c’}$，然后再求积分。我们用线性代数，将这个结构厘清，又简单利用微积分，将结果得出。

对于每个$c_k’$，我们有

$$c_k'=\frac{W_k}{W}$$

从而$c_k = \int\frac{W_k}{W}dx$，重复$n$次，就得到了整个$y_p$，从而得到我们想要的解。其中$W$为上式中由$y_1,y_2,\cdots,y_n$构成的矩阵的行列式，而$W_k$是将$W$中第$k$列换成$(0,0,\cdots,0,f)^T$之后得到的矩阵的行列式。具体而言，我们有

$$y=\sum_{k=1}^{n}C_ky_k+\sum_{k=1}^{n}y_k\int\frac{W_k(x)}{W(x)}dx$$

并没有什么诡异又空穴来风的公式需要记忆，这一切都来得很自然。可能会有人注意到，如果这里的$W(x)=0$该怎么办？这其实说明，选择的解，并没有代表性，其中的一个解，可以用另外几个解表示出来。比如方程$y’’=y$，应该选择解$e^{x},e^{-x}$，不应该选择$e^{x}$和$2e^{x}$，因为这个时候$e^{x}$可以很直接地表示出$2e^{x}$来。

一个例子

这篇博客里并没有直接解释怎么直接解出高阶常微分方程，实际上这个问题并不见得有多简单。但是，如果我们有了超能力，能解出方程右边等于$0$时候的解，就能得到所有实际要求的解了。我们可以试着解一个方程。

$$y''+y-2y=x$$

如果要解$y’’+y-2y=0$，可以得到一组具有代表性的解：$e^{x}$和$e^{-2x}$。这个时候的$W$和$W_k$就可以直接求出来了。

$$W(x)=\begin{vmatrix}e^x&e^{-2x} \\e^x&-2e^{-2x}\end{vmatrix}=-3e^{-x}$$ $$W_1(x)=\begin{vmatrix}0 & e^{-2x} \\x&-2e^{-2x}\end{vmatrix}=-xe^{-2x}$$ $$W_2(x)=\begin{vmatrix}e^x&0\\e^x&x\end{vmatrix}=xe^x$$

从而

$$y=e^x(C_1+\int_0^x\frac{s}{3}e^{-s}ds)+e^{-2x}(C_2+\int_0^x\frac{-s}{3}e^{2s}ds)$$

假设我们知道在原点处的取值情况，例如$\begin{cases}y(0)&=0\\ y’(0)&=1\end{cases}$，我们就能得到方程组

$$\begin{cases} C_1+C_2&=0 \\ C_1-2C_2&=1 \end{cases}$$

解出$C_1=\frac{1}{3}$，$C_2=-\frac{1}{3}$.