利用微分算子法求解常系数非齐次线性微分方程的特解
概述
$$x^{(n)}(t)+a_{1}x^{(n-1)}(t)+\ldots+a_{n-1}x'(t)+a_{n}x(t)=f(t)\cdots\cdots(eq1)$$
高阶常系数非齐次线性微分方程特解的计算是高阶线性微分方程部分的重要内容, 目前大多数教材都是以比较系数法（待定系数法）为核心方法, 部分教材介绍了拉普拉斯变换的方法, 实际上, 还有一个比较小众但是效率很高的求解方法---微分算子法一直被忽视, 下面我对这个方法作简单的介绍, 并辅以几个例题供大家加深对此方法的理解.
微分算子基本概念
首先我们定义微分算子$D=\dfrac{{\rm d}}{{\rm d}t}$, $\ldots$, $D^{n}=\dfrac{{\rm d}^{n}}{{\rm d}t^{n}}$, 这样一来我们可以将$(eq1)$重新改写为
$$D^{n}x(t)+D^{n-1}(a_{1}x(t))+\ldots+D(a_{n-1}x(t))+a_{n}x(t)=f(t).$$
由于求导运算是满足线性性质的, 所以不难知道对于任意的常数$a$, $b$和$k$阶可导函数$x(t),y(t)$, 都有
$$D^{k}(a x(t)+b y(t))=a D^{k}x(t)+b D^{k} y(t).$$
所以$(eq1)$可以重新记为：
$$(D^{n}+a_{1}D^{n-1}+\ldots+a_{n-1}D+a_{n})x(t)=f(t).$$
进一步, 我们定义
$$P(D)=D^{n}+a_{1}D^{n-1}+\ldots+a_{n-1}D+a_{n},$$
则可将$(eq1)$表示为：$P(D)x(t)=f(t)$.
此时此刻我们都很希望直接得到$x(t)=\dfrac{1}{P(D)}f(t)$的结果, 但需要注意的是此时的$\dfrac{1}{P(D)}$又是一个新的符号, 我们可以用这种形式给它下定义：若$P(D)x(t)=f(t)$, 则$x(t)=\dfrac{1}{P(D)}f(t)$. 换言之, $\dfrac{1}{P(D)}$只是一个符号, 它表示$P(D)$的逆运算, 为了实际计算$\dfrac{1}{P(D)}f(t)$, 我们必须要分析$P(D)$和$\dfrac{1}{P(D)}$的运算性质.
微分算子的运算性质
下面我们给出微分算子的几个常用性质, 这些性质的证明可以直接根据求导的相关运算得到, 这里不给出详细的证明过程了.
$(1)$ $D^{n}{\rm e}^{at}=a^{n}{\rm e}^{at}$, 注意, 这里的$a$可以是复数.
$(2)$ $D^{n}\sin(t)=\sin(t+\dfrac{n{\rm \pi}}{2})$, $D^{n}\cos(t)=\cos(t+\dfrac{n{\rm \pi}}{2})$.
以上这几个性质是针对具体函数的运算结果, 下面这个性质具有一般性, 也是利用微分算子求解常微分方程极其重要的性质：
$(3)$ $D^{n}[{\rm e}^{at}u(t)]={\rm e}^{at}(D+a)^{n}u(t)$.
性质$(3)$的证明可借助数学归纳法完成, 首先当$n=1$时直接计算可以得到相关结果, 然后假设当$n=k-1$时等式成立, 接着讨论$n=k$的情形即可.
$(4)$ $(D-\alpha)(D-\beta)x(t)=(D-\beta)(D-\alpha)x(t)=[D^{2}-(\alpha+\beta)D+\alpha\beta]x(t)$.
上述4个性质是接下来$P(D)$的运算性质的基础.
$P(D)$的运算性质
下面我们给出$P(D)$的运算性质, 首先我们给出$P(D)$的线性性质：
$(1)$ $P(D)[\alpha x(t)+\beta y(t)]=\alpha P(D)x(t)+\beta P(D)y(t)$.
这个性质直接根据微分算子的线性性质可以得到. 事实上, $P(D)$是关于$D$的多项式, 它也遵循下述的多项式性质.
$(2)$ 若$P(\lambda)=0$有$n$个根(不计重数)$\lambda_{i}$, $i=1,2,\ldots,n$, 则$P(D)=\prod\limits_{i=1}^{n}(D-\lambda_{i})$.
这个性质可以根据微分算子的性质$(4)$得到. 结合$P(D)$的性质$(1)$和$(2)$以及微分算子的性质$(1),(2)$我们不难得到如下性质：
$(3)$ $P(D){\rm e}^{at}=P(a){\rm e}^{at}$.
$(4)$ $P(D^{2})\sin(\beta t)=P(-\beta^{2})\sin(\beta t)$, $P(D^{2})\cos(\beta t)=P(-\beta^{2})\cos(\beta t)$.
这里性质$(4)$的证明可以取$a={\rm i}\beta$代入性质$(3)$, 分别考虑实部虚部即可. 微分算子的性质$(3)$也可以直接推广到$P(D)$上：
$(5)$ $P(D)[{\rm e}^{at}u(t)]={\rm e}^{at}P(D+a)u(t)$（移位性质）.
有了$P(D)$的性质我们可以根据$\dfrac{1}{P(D)}$的定义得到它的运算性质.
$\dfrac{1}{P(D)}$的运算性质
$\dfrac{1}{P(D)}$本身只是一个符号, 重要的是这个符号有什么样的性质可以用来支撑我们的计算, 下面列出几个常用的性质, 需要熟练掌握方可在求解$(eq1)$的特解时事半功倍.
$(1)$ $\dfrac{1}{D}x(t)=\int x(t){\rm d}t$, $\dfrac{1}{D^{n}}x(t)=\int \dfrac{1}{D^{n-1}}x(t){\rm d}t$.
在这里, 计算特解的过程不需要加上任意常数, 也就是说右端的不定积分符号表示被积函数的某一个原函数, 例如$\dfrac{1}{D}t^2=\dfrac{t^{3}}{3}$即可.
$(2)$ 若$P(a)\neq0$, 则$\dfrac{1}{P(D)}{\rm e}^{at}=\dfrac{1}{P(a)}{\rm e}^{at}$.
$(3)$ 若$P(-\beta^{2})\neq0$, 则$\dfrac{1}{P(D^{2})}\sin(\beta t)=\dfrac{1}{P(-\beta^{2})}\sin(\beta t)$, $\dfrac{1}{P(D^{2})}\cos(\beta t)=\dfrac{1}{P(-\beta^{2})}\cos(\beta t)$.
当$(eq1)$的右端函数$f(t)$只包含指数函数或者三角函数时, 性质$(2)$和$(3)$给我们的计算带来了极大的便利, 参看例1, 例2, 例3, 例4. 不难发现, 这里比较强的限制就是$P(a)\neq0$和$P(-\beta^{2})\neq0$, 如果遇到$P(a)=0$或者$P(-\beta^{2})=0$该如何解决呢？下面的性质可以部分解决这个问题.
$(4)$ 若$a$是$P(\lambda)=0$的$k$重根, 则$\dfrac{1}{P(D)}{\rm e}^{at}=\dfrac{t^{k}}{P^{(k)}(a)}{\rm e}^{at}$.
相应的如果右端是三角函数的形式, 我们可以借助例4的解法二(即实部虚部分离的方法)求解, 参看例5, 例6. 到此为止右端函数只包含指数函数或者三角函数我们都可以解决了, 但是如果右端出现了多项式似乎还解决不了, 下面两个性质结合起来则可以解决右端出现多项式的情形.
$(5)$ $\dfrac{1}{P(D)}[{\rm e}^{at}u(t)]={\rm e}^{at}\dfrac{1}{P(D+a)}u(t)$.
这个性质的证明可以直接将算子$P(D)$作用到右端函数${\rm e}^{at}\dfrac{1}{P(D+a)}u(t)$上, 进而根据移位性质可得结果. 较为常见的情形是$u(t)$是多项式函数, 可以看到即使利用这个性质, 最终还是有$\dfrac{1}{P(D+a)}u(t)$需要面对, 到目前为止我们只会计算$\dfrac{1}{D^{k}}x(t)$, 即积分$k$次, 也就是说即使是$\dfrac{1}{D-1}t^{2}$这样简单的表达我们都无法计算. 事实上这个式子的计算方法不唯一, 这边介绍一种幂级数展开的方法, 即如下性质.
$(6)$ $\dfrac{1}{1-D}=1+D+D^{2}+D^{3}+\ldots$, $\dfrac{1}{1+D}=1-D+D^{2}-D^{3}+\ldots$.
如果$u(t)$是$k$次多项式, 性质(6)的展开只需要到$D^{k}$即可, 因为更高阶的导数都是零. 综合前面这些性质, 对于多项式和指数函数的乘积类型的$f(t)$, 微分算子法可以彻底解决了. 参看例7, 例8, 例9.
例题（求特解）
相应的齐次方程的通解(基本解组)的求解可以根据特征方程法计算, 难度相对小一点, 所以这里的例题只要求计算非齐次线性方程的特解.
例1 $x'(t)-3x(t)={\rm e}^{t}$.
解 原方程可改写为$(D-3)x(t)={\rm e}^{t}$, 所以$x(t)=\dfrac{{\rm e}^{t}}{D-3}=\dfrac{{\rm e}^{t}}{1-3}=-\dfrac{{\rm e}^{t}}{2}$.
例2 $x^{(4)}(t)+2x^{(3)}(t)+2x(t)={\rm e}^{-t}$.
解 原方程可改写为$(D^{4}+2D^{3}+2)x(t)={\rm e}^{-t}$, 所以$x(t)=\dfrac{{\rm e}^{-t}}{D^{4}+2D^{3}+2}={\rm e}^{-t}$.
例3 $x'(t)-x(t)=\sin(2t)$.
解 原方程可改写为$(D-1)x(t)=\sin(2t)$, 所以$x(t)=\dfrac{\sin(2t)}{D-1}=\dfrac{D+1}{D^{2}-1}\sin(2t)$. 由于$(D+1)\sin(2t)=2\cos(2t)+\sin(2t)$, 所以$x(t)=\dfrac{2\cos(2t)}{D^{2}-1}+\dfrac{\sin(2t)}{D^{2}-1}=-\dfrac{1}{5}(2\cos(2t)+\sin(2t))$.
例4 $x''(t)+2x'(t)+2x(t)=\cos(t)$.
解法一 易知$x(t)=\dfrac{\cos(t)}{D^{2}+2D+2}=\dfrac{\cos(t)}{1+2D}=\dfrac{2D-1}{4D^2-1}\cos(t)=\dfrac{1}{5}(2\sin(t)+\cos(t))$. 这里我们可以看到, 如果确定了$P(-\beta^{2})\neq0$, 那么分母上的所有$D^2$都可以直接换成$-\beta^2$, 最后剩下的关于$D$的一次式可根据平方差公式凑出$D^{2}$, 再替换. 对于右端是三角函数的方程我们也可以考虑用指数函数替换, 再取实部或者虚部.
解法二 考虑方程$z''(t)+2z'(t)+2z(t)={\rm e}^{{\rm i}t}$, 则$x(t)$可取$z(t)$的实部. 易知$z(t)=\dfrac{{\rm e}^{{\rm i}t}}{D^{2}+2D+2}$, 根据右端是指数函数时$\dfrac{1}{P(D)}$的性质可得$z(t)=\dfrac{{\rm e}^{{\rm i}t}}{2{\rm i}+1}=-\dfrac{1}{5}(2{\rm i}-1)(\cos t+{\rm i}\sin(t))$, 化简可得$z(t)$的实部$x(t)=\dfrac{1}{5}(2\sin(t)+\cos(t))$.
例5 $x'''(t)-2x''(t)+x'(t)={\rm e}^{t}$.
解 易知$x(t)=\dfrac{{\rm e}^{t}}{D(D-1)^{2}}$, 这里$a=1$是$P(\lambda)=\lambda(\lambda-1)^2$的$2$重根, 且$P^{(2)}(1)=2$, 所以$x(t)=\dfrac{t^{2}{\rm e}^{t}}{2}$.
例6 $x^{(4)}(t)+2x''(t)+x(t)=\sin(t)$.
解 考虑方程$z^{(4)}(t)+2z''(t)+z(t)={\rm e}^{{\rm i}t}$, 则$x(t)$可取$z(t)$的虚部. 易知$z(t)=\dfrac{{\rm e}^{{\rm i}t}}{(D^{2}+1)^{2}}$, 这里$a={\rm i}$是$P(\lambda)=(\lambda^{2}+1)^2$的$2$重根, 且$P^{(2)}({\rm i})=-8$, 所以$z(t)=-\dfrac{t^{2}{\rm e}^{{\rm i}t}}{8}$. 化简可得$z(t)$的虚部$x(t)=-\dfrac{t^{2}}{8}\sin(t)$.
例7 $x'''-3x''+3x'-1={\rm e}^{t}(t^2-1)$.
解 显然, $x(t)=\dfrac{1}{(D-1)^{3}}{\rm e}^{t}(t^2-1)={\rm e}^{t}\dfrac{1}{(D+1-1)^{3}}(t^2-1)={\rm e}^{t}\dfrac{1}{D^{3}}(t^2-1)$. 对函数$(t^{2}-1)$作三次积分(忽略常数项)可得$\dfrac{1}{60}t^{3}(t^2-10)$. 所以$x(t)=\dfrac{1}{60}{\rm e}^{t}t^{3}(t^2-10)$.
注 本题如果考虑比较系数法, 则应该假设$x(t)=t^{3}(At^2+Bt+C){\rm e}^{t}$, 然后代回方程比较系数求出$A,B,C$, 希望读者能自己尝试计算一下, 比较微分算子法和经典的比较系数法的区别.
例8 $x''(t)-3x'(t)+2x(t)={\rm e}^{t}(t^2-1)+{\rm e}^{2t}+{\rm e}^{3t}$.
解 显然, $x(t)=\dfrac{1}{D^{2}-3D+2}{\rm e}^{t}(t^2-1)+\dfrac{1}{D^{2}-3D+2}{\rm e}^{2t}+\dfrac{1}{D^{2}-3D+2}{\rm e}^{3t}$. 下面分别计算后面三项. 为方便起见, 我们令$P(D)=D^{2}-3D+2=(D-1)(D-2)$.
所以$\dfrac{1}{D^{2}-3D+2}{\rm e}^{3t}=\dfrac{1}{P(3)}{\rm e}^{3t}=\dfrac{1}{2}{\rm e}^{3t}$. 且不难得到$\dfrac{1}{D^{2}-3D+2}{\rm e}^{2t}=\dfrac{t{\rm e}^{2t}}{P'(2)}=t{\rm e}^{2t}$.
$\dfrac{1}{D^{2}-3D+2}{\rm e}^{t}(t^2-1)=\dfrac{1}{(D-1)(D-2)}{\rm e}^{t}(t^2-1)={\rm e}^{t}\dfrac{1}{(D+1-1)(D+1-2)}(t^2-1)={\rm e}^{t}\dfrac{1}{D(D-1)}(t^2-1)$.
因为$\dfrac{1}{D-1}(t^{2}-1)=-(1+D+D^{2})(t^{2}-1)=-t^{2}-2t-1$, 所以$\dfrac{1}{D(D-1)}(t^2-1)=-\dfrac{1}{D}(t^{2}+2t+1)=-\dfrac{1}{3}t^3-t^2-t$.
综上所述, $x(t)=-t(\dfrac{1}{3}t^2+t+1){\rm e}^{t}+t{\rm e}^{2t}+\dfrac{1}{2}{\rm e}^{3t}$.
例9 $x''(t)-2x'(t)+2x(t)=t{\rm e}^{t}\cos(t)$.
解 考虑方程$z''(t)-2z'(t)+2z(t)=t{\rm e}^{(1+{\rm i})t}$, 则$x(t)$可取$z(t)$的实部.
易知$z(t)=\dfrac{1}{D^2-2D+2}t{\rm e}^{(1+{\rm i})t}={\rm e}^{(1+{\rm i})t}\dfrac{1}{(D+{\rm i}+1)^2-2(D+{\rm i}+1)+2}t$.
因为$(D+{\rm i}+1)^2-2(D+{\rm i}+1)+2=D^2+2{\rm i}D$.
所以$z(t)={\rm e}^{(1+{\rm i})t}\dfrac{1}{D^2+2{\rm i}D}t={\rm e}^{(1+{\rm i})t}\dfrac{1}{D}\dfrac{1}{D+2{\rm i}}t$.
根据微分算子幂级数展开的性质可知：
$\dfrac{1}{D+2{\rm i}}t=\dfrac{1}{2{\rm i}}\dfrac{1}{1+\dfrac{1}{2{\rm i}}D}t=\dfrac{1}{2{\rm i}}(1-\dfrac{1}{2{\rm i}}D)t=-\dfrac{1}{2}{\rm i}t+\dfrac{1}{4}$. 而$\dfrac{1}{D}(-\dfrac{1}{2}{\rm i}t+\dfrac{1}{4})=-\dfrac{1}{4}{\rm i}t^2+\dfrac{1}{4}t$.
综上, $z(t)={\rm e}^{t}(\cos(t)+{\rm i}\sin(t))\left(-\dfrac{1}{4}{\rm i}t^2+\dfrac{1}{4}t\right)$.
所以$z(t)$的实部$x(t)=\dfrac{1}{4}t{\rm e}^{t}\cos(t)+\dfrac{1}{4}t^2{\rm e}^{t}\sin(t)$.
注 本题过程较为繁琐, 读者可尝试用比较系数法计算, 即设$x(t)={\rm e}^{t}\left([at+b]\cos t+[ct+d]\sin t\right)$, 再代回方程比较系数求出$a,b,c,d$, 两相比较自然会感受到微分算子法的高效率.
注 如果出现$\dfrac{1}{D(D+a)}x(t)$, 我们一般写成$\dfrac{1}{D}\dfrac{1}{(D+a)}x(t)$形式计算, 即后积分, 否则中间过程用幂级数展开的方法会出现不必要的高阶导$D^{k}$, 使得最终结果出现一项常数, 这一项是否存在不会影响最终的通解, 但是会使得计算过程变得麻烦一些, 读者可尝试将例9中的$\dfrac{1}{D}\dfrac{1}{D+2{\rm i}}t$改成$\dfrac{1}{D+2{\rm i}}\dfrac{1}{D}t$计算, 再与例9比较.
总结
以上是我对微分算子法的基本概念、运算性质以及如何求特解的一个梳理, 如果读者切实动笔计算, 会发现如果右端函数出现了指数函数和三角函数, 总的来说用微分算子法会很大程度提高计算效率, 但是如果右端只有多项式, 则可以考虑比较系数法, 因为多项式的求导还是很方便的. 因此在实际的问题中需要灵活运用这几种不同的方法：比较系数法、微分算子法、拉普拉斯变换三种方法能做的类型仅限于右端非齐次函数是指数函数、三角函数与多项式的乘积或者求和类型, 如果右端不是这几种类型, 大家不能忘记最原始、最全能的常数变易法, 正如全科大夫在某个具体病症上的诊治水平往往没有专科大夫高一样, 常数变易法虽然理论上右端任意的$f(t)$都可以计算, 但是计算过程相对繁琐一点. 能力一般, 水平有限, 希望简单的梳理对大家理解和学习微分算子法能起到作用.
王飞
2023年5月11日
