关于常系数高阶线性方程中所使用的辅助方程法
2022年01月07日
这是一篇小报告，在上交之后我发表在这里。
对于常系数高阶线性方程的第二类待定系数法解，其的计算实在过于麻烦。于是[1]中提及一种简便方法来得到结果，那就是辅助方程法。但不管是网络上还是其他书籍里都少有相关方法的介绍，所以我在这里做一个简单的解释和总结。
1 原理
简单来说这是对欧拉恒等式的灵活运用，当在实数域时，一个a+bi可以看作Rank=2，但在复数域内则有Rank=1。
辅助方程法的一个出发点是叠加原理，如下所示：
定理 1 (叠加原理). 设L[y]为线性微分算子，若：
y_1(x)是L[y]=f_1(y)的解；
y_2(x)是L[y]=f_2(y)的解；
则有y_1(x)+y_2(x)是L[y]=f_1(x)+f_2(x)的解。
通过叠加原理我们可以组合不同的方程来得到一个组合解，与下面的独立性原理组合就能导出辅助方程方法。
定理 2 (独立性). 对于a+bi(a,b\in\mathbb R)，其中a与b线性无关。
对于常微分方程这门课来说，我们利用复函数的目的是导出实函数解。考虑第二类其实是由第一类转换而得，显然我们可以通过第一类待定系数法来简化计算。
2 辅助方程法的步骤与构造方法
我们以y''+y=\frac12\cos x为例。
例 1. 求解y''+y=\frac12\cos x。
解. 对于\frac12\cos x，我们可以对\frac12e^{ix}展开得到。那么我们就设一方程为y''+y=\frac12e^{ix}。
首先求得特征根\lambda_{1,2}=\pm i，也即i是一重特征根，由第一类方法即得解 \phi_1(x)=-\frac i4xe^{ix}=\frac x4\sin x-\frac i4x\cos x 由独立性，显然y''+y=\frac12e^{ix}=\frac12\cos x+\frac i2\sin x在\mathbb R和\mathbb C\backslash\mathbb R\cup \{0\}上有分别的解。这样就可以得出\frac x4\sin x是y''+y=\frac12\cos x的特解，而\frac x4\cos x是y''+y=\frac12\sin x的特解。 ◻
对于同样的例子，如果我们使用经典的第二类待定系数法来解，则会得到以下过程：
由一重特征根设特解为\phi(x)=x(A\cos x+B\sin x)，对式子求导得\phi''(x)=(2B-Ax)\cos x-(A+Bx)\sin x。由是得\phi''+\phi=2B\cos x+(Bx-A-Bx)\sin x=\frac12\cos x，由\cos x与\sin x的正交性可得方程组 \begin{aligned} 2B&=\frac12\\ (-A)&=0\end{aligned} 由是解得B=\frac14，即有同上的解。
上面的方程是一个简单的例子，我们看一个稍微复杂的例子来进行分析。通过这个例子我们将更加明了地认识到辅助方程法的优越性。
我们以y''-2y'+2y=xe^x\cos x为例。
例 2. 求解y''-2y'+2y=xe^x\cos x。
解. 对于xe^x\cos x，很显然地我们找到xe^{(1+i)x}，他通过Euler等式即可化为xe^x\cos x+ixe^x\sin x。
由是我们解方程y''-2y'+2y=xe^{(1+i)x}。
求解特征根得\lambda_{1,2}=1\pm i，由是1+i是一重特征根，即得解的形式为\phi(x)=x(ax+b)e^{(1+i)x}，代入得解为 \begin{aligned} \phi(x)&=x\left(-\frac i4x+\frac14\right)e^{(1+i)x}\\&=\frac x4e^x\left(\cos x+x\sin x+i(\sin x-x\cos x)\right) \end{aligned} 而xe^x\cos x对应的解即是上面根的实部。所以原方程特解为 \phi(x)=\frac x4e^x(\cos x+x\sin x) ◻
对于该方程我们用第二类待定系数法去解：
由于1+i是一重特征根，则我们设解为\phi(x)=xe^x[(Ax+b)\cos x+(Cx+d)\sin x]，则有 \begin{gathered} \phi'(x)=e^x ((-A x^2 - b x + C x^2 + 2 C x + d x + d)\sin x + \\(A x^2 + 2 A x + b x + b + C x^2 + d x)\cos x)\end{gathered} \begin{gathered} \phi''(x)=-2 e^x ((A x^2 + 2 A x + b x + b - 2 C x - C - d)\sin x + \\(-2 A x - A - b - C x^2 - 2 C x - d x - d)\cos x)\end{gathered} 由是有 \begin{gathered} \phi''(x)-2\phi'(x)+2\phi(x)=\\ 2 e^x ((-2 A x - b + C)\sin x + (A + 2 C x + d)\cos x)=xe^x\cos x\end{gathered} 可列方程组 \left\{ \begin{aligned} -4A&=0\\ 2(-b+C)&=0\\ 2(A+d)&=0\\ 4C&=1 \end{aligned} \right. 由是得解。
3 求解过程的复杂程度分析
我们可以用简单的迭代知识求得对一个式子朴素的求导次数。通过这种方法我们可以进一步考虑辅助方程法的优越性。
3.1 对于第一类换元法
我们记有\psi_{p,k}=\frac{\mathrm d{\psi_{p,k-1}}}{\mathrm dx},
对于项\psi_{m,0}=Q_{m,0}(x)e^{\mu x}，通过求导的乘法公式可以得出我们需要做m次求导操作与m-1次加法，得到一个新的式子 \frac{\mathrm d{\psi_{m,0}}}{\mathrm dx}=\psi_{m,1}=Q_{m,1}(x)e^{\mu x}.
我们设Q_{m,0}(x)=\sum_{k=0}^m{a_kx^k}，则有 \begin{aligned} Q_{m,0}(x)&=\sum_{k=0}^m{a_kx^k}\\ Q_{m,1}(x)&=\sum_{k=0}^{m-1}{(a_k+\mu a_{k+1})x^k}+a_mx^m\\ Q_{m,2}(x)&=Q_{m,1}'+Q_{m,1}=Q_{m,0}''+2Q_{m,0}'+Q_{m,0}\\ \cdots\\ Q_{m,k}(x)&=\sum_{j=0}^k\binom{k}{j}Q_{m,0}^{(j)}\\ \cdots\end{aligned} 对于Q_{m,0}^{(j)}，其的每一项的系数为\frac{k!}{(k-j+1)!},k=j,\cdots,m。统共要做j(m-j+1)次乘法。对于系数又要做\min\{j,m-j\}次乘法。则求Q_{m,k}(x)要做\min\{j,m-j\}m(m-j+1)次乘法，合并\sum_{j=1}^mj(m-j)=\frac16m(1+m)(m-1)次同类项。最坏的情况下我们需要计算\sum_{l=0}^{m}b_lQ_{m,l},b_l\not=0。则我们要做\sum_{j=0}^m\min\{j,m-j\}m(m-j+1)次乘法，合并\frac16m^3(1+m)(m-1)次同类项。
最后解方程组，若通过高斯消元法则需要O(m^3)的时间复杂度消元。
如使m=3，则要做33+36+27次基本操作，假设每次操作有.1\%可能性出错。则P_{WRONG}=1-(.999)^{96}=.092。
3.2 对于第二类换元法
当使用第二类换元法时，由于要设两组Q_{m,0}，则m看作翻倍，但Q_{m,k}\equiv0,k>\frac m2，如是我们可以修正得如下表。
| 第一类 | 第二类 | |
|---|---|---|
| 乘法 | \sum_{j=0}^m\min\{j,m-j\}m(m-j+1) | \sum_{j=0}^m\min\{j,2m-j\}2m(2m-j+1) |
| 合并同类项 | \frac16m^3(1+m)(m-1) | \frac13m^3(1+2m)(2m-1) |
| 解方程 | O(m^3) | O(m^3) |
取m=3，则要做162+315+27次基本操作，同样假设每次操作有.1\%可能性出错。
则P_{WRONG}=1-(.999)^{505}=.396，容易犯错。
4 同方法在其他范围内的应用
由e^x的级数定义与Taylor展开很容易证明 \frac{\mathrm d{e^{\omega x}}}{\mathrm dx}=\omega e^{\omega x}, x\in\mathbb{R}, \omega\in\mathbb C 所以我们在对e^{\omega x}的运算中可以把他当作是实函数去计算以简化操作。
4.1 概率论的特征函数
参考[2]，实际上这是对原密度函数做Fourier变换，也即 \phi(t)=\int_{-\infty}^{\infty}e^{itx}p(x)dx
例 3. 求Exp(\lambda)的特征函数.
Proof. 考虑[expon]，显然地有\frac{\mathrm d{\frac{e^{\omega x}}{\omega}}}{\mathrm dx}=e^{\omega x}, x\in\mathbb{R}, \omega\in\mathbb C。那么就有不定积分\int e^{\omega x}\mathrm dx=\frac{e^{\omega x}}{\omega}+C。
而Exp(\lambda)的PDF为 p(x)=\left\{\begin{aligned} \lambda e^{-\lambda x}&,x\ge 0,\\ 0&,x<0, \end{aligned}\right. 由是有 \begin{aligned} \phi(t)&=\int_0^{+\infty}e^{itx}\lambda e^{-\lambda x}\mathrm dx\\ &=\lambda\int_0^{+\infty}e^{(it-\lambda)x}\mathrm dx\\ &=-\frac\lambda{it-\lambda} \end{aligned} 1 附上教材解法： \begin{aligned} \phi(t)&=\int_0^{+\infty}e^{itx}\lambda e^{-\lambda x}\mathrm dx\\ &=\lambda\left[\int_0^{+\infty}\cos{(tx)}e^{-\lambda x}\mathrm dx+i\int_0^{+\infty}\sin{(tx)}e^{-\lambda x}\mathrm dx\right]\\ &=\lambda\left(\frac\lambda{\lambda^2+t^2}+i\frac t{\lambda^2+t^2}\right)=-\frac\lambda{it-\lambda} \end{aligned} ◻
正态分布是一个更好的例子，我们通过变换形式而避免将函数展开。
例 4. 求N(0,1)的特征函数.
Proof. N(0,1)的PDF为 p(x)=\frac1{\sqrt{2\pi}}e^{-\frac{x^2}2}, -\infty<x<+\infty, 由是其特征函数为 \begin{aligned} \phi(t)&=\int_{-\infty}^{+\infty}e^{itx}\frac1{\sqrt{2\pi}}e^{-\frac{x^2}2}\mathrm dx\\ &=e^{-\frac{t^2}2}\int_{-\infty}^{+\infty}\frac1{\sqrt{2\pi}}e^{-\frac{(x-it)^2}2}\mathrm d(x-it)\\ &=e^{-\frac{t^2}2} \end{aligned} 2 同样附上教材解法： \begin{aligned} \phi(t)&=\int_{-\infty}^{+\infty}e^{itx}\frac1{\sqrt{2\pi}}e^{-\frac{x^2}2}\mathrm dx\\ &=\frac{1}{\sqrt{2\pi}} \int_{-\infty}^{+\infty} \sum_{n=0}^\infty \frac{(itx)^n}{n!} e^{-\frac{x^2}{2}}\mathrm dx\\ &=\sum_{n=0}^\infty\frac{(it)^n}{n!}\left[\frac{1}{\sqrt{2\pi}}\int_{-\infty}^{+\infty}x^ne^{-\frac{x^2}{2}}\mathrm dx\right] \end{aligned} 方括号内是标准正态分布的n阶矩E(X^n)，当n为奇数时由对称性E(X^n)=0；当n为偶数时，取n=2m，则有 E(X^n)=E(X^{2m})=(2m-1)!!=\frac{(2m)!}{2^m\cdot m!} 由是可得 \phi(t)=\sum_{m=0}^\infty\frac{(it)^2m}{(2m)!}\frac{(2m)!}{2^m\cdot m!}=\sum_{m=0}^\infty\left(-\frac{t^2}{2}\right)^m\frac1{m!}=e^{-\frac{t^2}2} ◻
4.2 Laplace变换具有类似形式
Laplace变换和Fourier变换在形式上类似，我们可以通过类似方法计算。
例 5. 计算f(t)=e^{at}的Laplace变换。
Proof. \begin{aligned} \int_0^{+\infty}e^{-st}e^{at}\mathrm dt&=\int_0^{+\infty}e^{(a-s)t}\mathrm dt\\ &=\frac{1}{a-s}\int_0^{+\infty}e^{(a-s)t}\mathrm d(a-s)t\\ &=\frac{1}{s-a} \end{aligned} 3 ◻