高阶线性方程

本节的研究对象是下面的 `n` 阶线性方程: `sum_(k=0)^n a_k(t) x^((k))(t) = f(t)`, `quad a_n = 1`. 我们重点掌握常系数方程的解法.

常数变易法

常数变易法是齐次方程到非齐次方程的桥梁. 具体来说, 要解方程 , 假设已知对应齐次方程 (令 `f = 0`) 的基本解组为 `{x_k(t)}_(k=1)^n`, 那么原方程的通解为 `x = sum_(k=1)^n c_k(t) x_k(t)`, 其中函数 `c_k(t)` 的导数满足线性方程组 `[ x_1, x_2, cdots, x_n; x_1', x_2', cdots, x_n'; vdots, vdots, , vdots; x_1^((n-1)), x_2^((n-1)), cdots, x_n^((n-1)) ]` `[c_1'; c_2'; vdots; c_n']` `= [0; 0; vdots; f(t)]`.

直接验证: `x = sum c_j x_j`,
`x' = sum c_j x_j' + c_j' x_j = sum c_j x_j'`,
`x'' = sum c_j x_j''`,
`cdots`
`x^((n)) = sum c_j x_j^((n)) + c_j' x_j^((n-1))` `= sum c_j x_j^((n)) + f`. 求和得 `sum_(k=0)^n a_k x^((k))` `= sum_(k=0)^n a_k sum_(j=1)^n c_j x_j^((k)) + f` `= f`.

受迫振动 解如下的常微分方程初值问题: `{ x'' + a^2 x = f(t); x"|"_(t=0) = x_0; x'"|"_(t=0) = x'_0; :}`

`{cos at, sin at}` 构成齐方程的基本解组. 设 `x = c_1 cos at + c_2 sin at`, 解线性方程组 `{ c_1' cos at + c_2' sin at = 0; c_1' sin at - c_2' cos at = -f//a; :}` 得 `(c_1', c_2') = f//a (-sin at, cos at)`. 由初值条件有 `c_2|_(t=0) = x_0'//a`, `quad c_1|_(t=0) = x_0`, 故 `c_1 = -1/a int_0^t f(tau) sin a tau "d"tau + x_0`,
`c_2 = 1/a (int_0^t f(tau) cos a tau "d"tau + x_0')`. 所以解为 `x = x_0 cos at + (x_0')/a sin at + 1/a int_0^t f(tau) sin[a(t-tau)] "d"tau`.

常系数齐次线性方程

常系数齐次线性方程 现在设 中的系数都是实常数, 且方程右边 `f = 0`. 方程化为 `x^((n)) + a_(n-1) x^((n-1)) + cdots + a_0 x = 0`, `quad a_k in RR`.

将猜测的解 `x = "e"^(lambda t)` 代入原方程, 得到方程的特征多项式: `lambda^n + a_(n-1) lambda^(n-1) + cdots + a_1 lambda + a_0`. 设特征根为 `{lambda_i}_(i=1)^l`, 重数分别为 `n_i`, 且 `sum_(i=1)^l n_i = n`. 则方程的基本解组为 `t^j "e"^(lambda_i t)`, `quad j = 0, 1, cdots, n_i-1`, `quad i = 1, 2, cdots, l`. 若 `lambda_i = alpha + "i" beta` 为复根, 则 `bar(lambda_i) = alpha - "i" beta` 也必为特征根, 可以将基本解组中的一对共轭复根 `"e"^((alpha +- "i" beta) t)` 换成 `"e"^(alpha t) cos beta t` 与 `"e"^(alpha t) sin beta t`.

Euler 方程 `t^n x^((n)) + a_(n-1) t^(n-1) x^((n-1)) + cdots + a_0 x = 0`.

将 `x = t^lambda` 代入方程, 得到 `n` 个特征根. 基本解组为 `ln^j |t| * t^(lambda_i)`.

常系数非齐次线性方程

非齐次方程 `x^((n)) + a_(n-1) x^((n-1)) + cdots + a_0 x = f(t)`, `quad a_k in RR` 的通解就是齐次方程的通解加上自身的一个特解; 因此我们着重介绍特解的求法.

多项式法 求解 `y''+p y'+q y = P_m(x) "e"^(lambda x)`, 其中 `p`, `q`, `lambda` 是常数, `P_m(x)` 是 `m` 次多项式.

令 `y = z "e"^(lambda z)`, 方程化为 `(F''(lambda))/(2!) z'' + (F'(lambda))/(1!) z' + F(lambda) z = P_m(x)`, 这里 `F(lambda) = lambda^2 + p lambda + q` 为方程的特征多项式. 新的方程有多项式形式的特解, 可用待定系数法求解.

升阶法 求解 `y''+p y'+q y = P_m(x)`, 右边是 `m` 次多项式.

记 `P_m(x) = sum_(k=0)^m a_k x^k`. 方程两边同时对 `x` 求导 `m` 次得 `y''' + p y'' + q y' = P_m'(x)`,
`cdots`
`y^((m+1)) + p y^((m)) + q y^((m-1)) = a_m m! x + a_(m-1) (m-1)!`,
`y^((m+2)) + p y^((m+1)) + q y^((m)) = a_m m!`. `q != 0` 时, 取最后一个方程的特解 `y^((m)) = q^(-1) a_m m!`, 和 `y^((m+1)) = 0` 一起代入倒数第二个方程求得 `y^((m-1))`, 如此依次代入前式, 最后可得方程的一个特解 `y(x)`.

一般而言, 多项式法与升阶法计算比较繁琐, 不如比较系数法.

    比较系数法 考虑以下形式的非齐次项:
  1. `f(t) = P_m(t) "e"^(lambda t)`, `P_m` 为 `m` 次多项式.
  2. `f(t) = P_m(t) "e"^(alpha t) cos beta t` 或 `f(t) = P_m(t) "e"^(alpha t) sin beta t`, `P_m` 为 `m` 次多项式.
    3. 利用以上两种类型的非齐次项, 结合线性叠加原理, 可写出更多非齐次项对应的特解.
  1. 此时有特解 `A_m(t) t^k "e"^(lambda t)`, `k` 为特征根 `lambda` 的重数, `A_m(t)` 是系数待定的 `m` 次多项式.
  2. 此时有特解 `(A_m(t) cos beta t + B_m(t) sin beta t) t^k "e"^(alpha t)`, `k` 为特征根 `lambda = alpha + "i" beta` 的重数, `A_m(t), B_m(t)` 为系数待定的 `m` 次多项式.
    3. 注意事项:

Lanchester 方程 传统战争中 `A, B` 两军作战, 设 `A` 对 `B` 的杀伤力为 `alpha^2`, `B` 对 `A` 的杀伤力为 `beta^2`, 则两军人数变化可以用下面的微分方程表示: `{ ("d"A)/dt = -beta^2 B; ("d"B)/dt = -alpha^2 A; :}` 其中 `A(0) = A_0`, `B(0) = B_0`.

令 `f' = B`, 且 `f(0) = 0`. 原方程组化为 `{ ("d"A)/dt = -beta^2 ("d"f)/dt; ("d"^2 f)/dt^2 = -alpha^2 A; :}` 第一式两边积分得 `A = -beta^2 f + A_0`, 代入第二式就有 `f'' = alpha^2 beta^2 f - alpha^2 A_0`, 它的解设为 `f = c_1 "e"^(alpha beta t) + c_2 "e"^(-alpha beta t) + c_3`. 代入前式得 `c_3 = A_0/beta^2`. 又 `B = f' = alpha beta (c_1 "e"^(alpha beta t) - c_2 "e"^(-alpha beta t))`, 由 `f(0) = 0` 和 `B(0) = B_0` 有 `{ c_1 + c_2 = -A_0/beta^2; c_1 - c_2 = B_0/(alpha beta); :}` 解得 `c_(1, 2) = 1/(2 beta) (-A_0/beta +- B_0/alpha)`. 最终 `A = (A_0 alpha cosh(alpha beta t) - B_0 beta sinh(alpha beta t))//alpha`,
`B = (-A_0 alpha sinh(alpha beta t) + B_0 beta cosh(alpha beta t))//beta`.

从相平面上看, 有 `("d"A)/("d"B) = (beta^2 B)/(alpha^2 A)`, 解为 `alpha^2 A^2 = beta^2 B^2 + c`, 是双曲线.

高阶方程的更多解法

Laplace 变换 也是求解微分方程的重要方法. 主要思想是将微分运算化为代数运算. 它的优点是无需求出通解, 可以将初值条件直接代入计算. 详情参见积分变换.

    降阶法 考虑齐次线性微分方程 `sum_(i=0)^n a_i(t) x^((i))(t) = 0`, 在以下几种情形可以降低方程阶数:
  1. 方程不显含小于 `k` 阶的导数时;
  2. 方程不显含自变元时;
  3. 已知方程的 `k` 个非零特解时.
  1. 令 `y = x^((k))`, 方程降低 k 阶.
  2. 令 `y = x'`, 降低 1 阶.
  3. 若有非零特解 `x_1, cdots, x_k`, 则令 `x = x_k int z dt`, 化为 `z` 的 `n-1` 阶方程, 附带 `k-1` 个特解 `z_i = (x_i//x_k)'`, `quad i = 1, 2, cdots, k-1`.

线性方程组

本节研究方程 `x' = A(t) x + f(t)`, 其中 `A` 为矩阵, `x, x', f` 为向量.

常数变易法

线性方程组的常数变易公式 解方程 `x' = A(t) x + f(t)`, `quad x(t_0) = eta`. 假设已知对应齐次方程 `x' = A(t) x` 的基解矩阵为 `Phi(t)`.

由已知 `Phi'(t) = A(t) Phi(t)`. 将 `x(t) = Phi(t) c(t)` 代入原方程得 `Phi(t) c'(t) = f(t)`. 因为 `Phi(t)` 可逆, 所以令 `c(t) = int_(t_0)^t Phi^-1(s) f(s) "d"s + Phi^-1(t_0) eta`, 就得到满足初值条件的解.

    高阶方程化为方程组
  1. 将微分方程 `sum_(i=0)^n a_i(t) x^((n-i)) = f(t)`, `quad a_0(t) = 1`,
    `x^((i))(t_0) = eta_(i+1)`, `quad i = 0, 1, cdots, n-1`     化为等价的方程组 `x' = A x + vec(f)(t)`, `quad x(t_0) = eta`.
  2. 设 1. 中方程相应的齐次方程有基本解组 `{x_k}_(k=1)^n`, 求原方程的解.
  1. `x = (x_1, x_2, cdots, x_n)^T` `= (x, x', cdots, x^((n-1)))^T`,
    `eta = (eta_1, eta_2, cdots, eta_n)^T`,
    `vec(f)(t) = (0, 0, cdots, 0, f(t))^T`,
    `A = [ 0, 1, 0, cdots, 0; 0, 0, 1, cdots, 0; vdots, vdots, vdots, , vdots; 0, 0, 0, cdots, 1; -a_n, -a_(n-1), -a_(n-2), cdots, -a_1; ]`.
  2. 由线性方程组的常数变易公式知, 原方程有特解 `x(t) = sum_(k=1)^n x_k(t) int_(t_0)^t (W_k(s))/(W(s)) f(s) "d"s`, 其中 `W(s) = | x_1(s), x_2(s), cdots, x_n(s); x_1'(s), x_2'(s), cdots, x_n'(s); vdots, vdots, , vdots; x_1^((n-1))(s), x_2^((n-1))(s), cdots, x_n^((n-1))(s); |`,
    `W_k(s)` 是 `W(s)` 中 `x_k^((n-1))(s)` 的代数余子式.     特别当 `n = 2`, 上式化为 `varphi(t) = int_(t_0)^t |x_1(s), x_2(s); x_1(t), x_2(t)| / |x_1(s), x_2(s); x_1'(s), x_2'(s)| f(s) ds`.

常系数线性方程组

常系数线性方程组 `x' = Ax`.

  1. 此方程有基解矩阵 `exp At`. 显然 `exp A t|_(t=0) = E`.
  2. 若 `lambda` 是 `A` 的特征值, `bm v` 为对应的特征向量, 则 `"e"^(lambda t) bm v` 是一个解. 特别地, 若 `A` 有 `n` 个线性无关的特征向量, 对应的特征值为 `lambda_1, cdots, lambda_n`, 则方程有基解矩阵 `["e"^(lambda_1 t) bm v_1, cdots, "e"^(lambda_n t) bm v_n]`.
  3. 若已知一个基解矩阵 `Phi(t)`, 则存在一常数矩阵 `C`, 使得 `exp At = Phi(t) C`, 再由 `exp A t|_(t=0) = E` 知 `exp At = Phi(t) Phi^-1(0)`.

利用根子空间求矩阵指数 矩阵指数是无穷级数, 直接计算比较困难: `exp bm A := sum_(k ge 0) bm A^k/(k!)`. 但如果存在某个 `k`, 使得 `bm A^k = bm O`, 则只需计算有限项. 此方法采用类似的思想, 只不过每次只求矩阵 `"exp"bm A` 的一列.

设 `lambda` 是 `bm A` 的特征值, 它的根子空间定义为 `V_lambda = { bm x in V | EE r_lambda in ZZ^+, (bm A - lambda bm E)^(r_lambda) bm v = 0 }`. 根据线性代数知识, 线性空间 `V` 在根子空间上具有直和分解, 即任意 `n` 维向量 `bm x` 可以唯一写成 `bm x = sum_lambda bm x_lambda`, `quad bm x_lambda in V_lambda`. 上式表示对 `bm A` 的不同特征值求和. 从而 `(exp bm A) bm x` `= sum_lambda (exp bm A) bm x_lambda`
`= sum_lambda "e"^lambda exp(bm A-lambda bm E) bm x_lambda`
`= sum_lambda "e"^lambda (sum_(k lt r_lambda) (bm A-lambda bm E)^k // k!) bm x_lambda`. 特别取 `bm x = bm epsi_j`, 就得到矩阵 `exp bm A` 的第 `j` 列.

求 `exp bm A t`, `bm A = [cos a, -sin a; sin a, cos a]`. 推论: `exp[, -a; a,] = [cos a, -sin a; sin a, cos a]`.

`bm A` 的特征值是 `lambda_(1,2) = cos a +- "i" sin a`, 先设 `lambda_1 = lambda_2 = cos a`, 则 `bm A = lambda bm E`, `exp(bm A - lambda bm E) = exp bm O = bm E`. 从而 `[1;0] exp bm A t = [1;0] "e"^(lambda t)`,
`[0;1] exp bm A t = [0;1] "e"^(lambda t)`,
`exp bm A t = "e"^(t cos a) bm E`. 再设 `lambda_1 != lambda_2`, 对应的特征向量为 `bm v_1 = ["i";1]`, `bm v_2 = [1;"i"]`. 将单位向量分解到特征子空间, `[1;0] = 1/2 [1;"i"] - "i"/2 ["i";1]`,
`[0;1] = -"i"/2 [1;"i"] + 1/2 ["i";1]`. 由于特征值都是一重的, 指数的无穷级数中实际上只有一项: `exp(bm A - lambda bm E) = bm E`. 从而 `[1;0] exp bm A t` `= 1/2 [1;"i"] "e"^(lambda_2 t) - "i"/2 ["i";1] "e"^(lambda_1 t)` `= 1/2 ["e"^(lambda_2 t) + "e"^(lambda_1 t); "i"("e"^(lambda_2 t) - "e"^(lambda_1 t))]` `= "e"^(t cos a) [cos(t sin a); sin(t sin a)]`. 同理可求得第二列, 于是 `exp bm A t` `= "e"^(t cos a) [cos(t sin a), -sin(t sin a); sin(t sin a), cos(t sin a)]`. 易知上式对重根的情形也成立.

用线性微分方程组的知识解 Lanchester 方程 `{ ("d"A)/dt = -beta^2 B; ("d"B)/dt = -alpha^2 A; :}`

矩阵 `A = [0, -beta^2; -alpha^2, 0]` 的特征值为 `+-alphabeta`, 对应的特征向量为 `(beta, ∓ alpha)^T`. 代入公式计算, `exp At (1;0) = 1/(2beta)("e"^(alpha beta t) (beta;-alpha) + "e"^(-alpha beta t) (beta;alpha))`,
`exp At (0;1) = 1/(2alpha)("e"^(alpha beta t) (-beta;alpha) + "e"^(-alpha beta t) (beta;alpha))`. 所以基解矩阵 `Phi(t) = exp At = [ cosh(alpha beta t), -beta/alpha sinh(alpha beta t); -alpha/beta sinh(alpha beta t), cosh(alpha beta t); ]`. 而 `Phi^(-1)(0) = [ cosh(alpha beta t), beta/alpha sinh(alpha beta t); alpha/beta sinh(alpha beta t), cosh(alpha beta t); ](0) = E`. 故 Lanchester 方程满足初值条件 `(A, B)^T(0) = (A_0, B_0)^T` 的解为 `Phi(t)Phi^(-1)(0) (A_0, B_0)^T = Phi(t)(A_0, B_0)^T`.

寻找级数 `sum_(n ge 0) x^(3n)/((3n)!)` 的闭形式.

记 `x_i = sum_(n ge 0) t^(3n+i)/((3n+i)!)`, `i = 0, 1, 2`, 于是 `x_0' = x_2`, `x_2' = x_1`, `x_1' = x_0`. 写成向量 `bm x = (x_0, x_1, x_2)^T`, 有 `bm x' = bm (A x)`, `quad bm A = [0, 0, 1; 1, 0, 0; 0, 1, 0]`. 容易看出 `exp bm A t = bm E x_0(t) + bm A x_1(t) + bm A^2 x_2(t)`, 但这只会使我们回到原点, 因此我们寻找它的另一种表示. `bm A` 的特征多项式为 `lambda^3 - 1`, 三个根为 `1, omega, omega^2`, 对应的特征向量为 `bm v_1 = [1; 1; 1]`, `bm v_2 = [1; omega^2; omega]`, `bm v_3 = [1; omega; omega^2]`. 注意到 `x_0(0) = 1, x_1(0) = x_2(0) = 0`, 因此方程组的初值条件恰为 单位向量 `bm eta = (1, 0, 0)^T`. 将这个向量分解为 `bm eta = 1/3 (bm v_1 + bm v_2 + bm v_3)`, 从而解得 `bm x` `= bm eta exp bm A t` `= 1/3 (bm v_1 "e"^t + bm v_2 "e"^(omega t) + bm v_3 "e"^(omega^2 t))` `= 1/3 [ "e"^t + 2 "e"^(t/2) cos((sqrt 3)/2 t); "e"^t + 2 "e"^(t/2) cos((sqrt 3)/2 t - (2pi)/3); "e"^t + 2 "e"^(t/2) cos((sqrt 3)/2 t + (2pi)/3) ]`.