微分方程、動力系統與混沌導論 第6章 高維線性系統[書摘]

第6章 高維線性系統函數

      在線性代數領域稍做停留後,如今該回到微分方程了,特別地,要回到求解具備常係數的高維線性系統的任務中來。和線性代數那一章同樣,咱們要討論大量的不一樣情形。學習

 

6.1 不一樣特徵值spa

      首先考慮線性系統$\boldsymbol X' = \boldsymbol {AX}$,其中$n \times n$矩陣$\boldsymbol A$具備$n$個不一樣的實特徵值$\lambda_1,\cdots,\lambda_n$。根據第5章中的結論,存在座標變換$\boldsymbol T$使得新系統$\boldsymbol Y' = (\boldsymbol T^{-1} \boldsymbol {AT})\boldsymbol Y$具備以下至關簡單的形式:3d

\[\begin{array}{l}{{y'}_1} = {\lambda _1}{y_1}\\\;\;\;\;\;\; \vdots \\{{y'}_n} = {\lambda _n}{y_n}\end{array}\]blog

其中的線性映射$T$將標準基向量$\boldsymbol E_j$映到屬於$\lambda_j$的特徵向量$\boldsymbol V_j$。顯然,下面形式的函數get

\[\boldsymbol Y(t)=\left( \begin{array}{l} c_1e^{\lambda_1t} \\  \vdots \\c_ne^{\lambda_nt}\end{array}\right)\]it

是$\boldsymbol Y' = (\boldsymbol T^{-1} \boldsymbol {AT})\boldsymbol Y$知足初值條件$\boldsymbol Y(0) = (c_1,\cdots,c_n)$的一個解。與第3章同樣,它是這樣的惟一解。變量

       由此可得$\boldsymbol X(t) = \boldsymbol {TX}(t)$就是$\boldsymbol X' = \boldsymbol {AX}$的通解,並且這個通解能夠寫成以下形式bfc

\[\boldsymbol X(t) = \sum \limits_{j=1}^nc_je^{\lambda_jt}\boldsymbol V_j.\]lambda

如今假設$\boldsymbol A$的特徵值中$\lambda_1,\cdots,\lambda_k$爲負的,而$\lambda_{k+1},\cdots,\lambda_n$爲正的。由於沒有零特徵值,故此時系統是雙曲的。首先,全部從$\boldsymbol V_1,\cdots,\boldsymbol V_k$張成的子空間出發的解都將永遠逗留在這個子空間中,這是由於此時$c_{k+1}=\cdots=c_n=0$。其次,每一個這樣的解在$t \to \infty$時都趨於原點。與平面系統中引入的術語相似,咱們稱這個子空間爲穩定子空間。相似地,$\boldsymbol V_{k+1},\cdots,\boldsymbol V_n$張成的子空間中的解將遠離原點。這個子空間稱爲不穩定子空間。全部其它的解在時間向後時都趨於穩定子空間,而在時間向前時都趨於不穩定子空間。於是這個系統對應於高維的鞍點

  考慮

\[\boldsymbol X' = \left( \begin{array}{l} 1&2&-1 \\ 0&3&-2\\0&2&-2 \end{array} \right) \boldsymbol X.\]

在5.2節中,咱們已經看到這個矩陣的特徵值爲$2,1,-1$,屬於它們的特徵向量分別爲$(3,2,1),(1,0,0),(0,1,2)$。因而矩陣

\[\boldsymbol T = \left( \begin{array}{l} 3&1&0 \\ 2&0&1 \\ 1&0&2 \end{array} \right) \]

就將 $\boldsymbol X' = \boldsymbol {AX}$化爲

\[\boldsymbol Y' = (\boldsymbol T^{-1} \boldsymbol {AT})\boldsymbol Y = \left( \begin{array}{l} 2&0&0 \\ 0&1&0 \\ 0&0&-1 \end{array} \right) \boldsymbol Y,\]

咱們能夠馬上獲得它的解,用$\boldsymbol T$乘以這個解就能夠獲得$\boldsymbol X' = \boldsymbol {AX}$的通解(若是單爲求通解的話,無需座標變換,由於前面已經求得了特徵值及對應的特徵向量,從而能夠直接寫出通解,此處做變換的目的是爲了化爲標準形)爲

\[\boldsymbol X(t) = {c_1}{e^{2t}}\left( \begin{array}{l}3\\2\\1\end{array} \right) + {c_2}{e^t}\left( \begin{array}{l}1\\0\\0\end{array} \right) + {c_3}{e^{ - t}}\left( \begin{array}{l}0\\1\\2\end{array} \right).\]

過原點以及$(0,1,2)$的直線就是穩定線,而由$(3,2,1),(1,0,0)$張成的平面是不穩定平面。該系統以及系統$\boldsymbol Y' = (\boldsymbol T^{-1} \boldsymbol {AT})\boldsymbol Y$的解如圖6.1所示。

image

  若是$3 \times 3$矩陣$\boldsymbol A$有3個不一樣的負特徵值,則咱們可做座標變換使得系統化成

\[\boldsymbol Y' = (\boldsymbol T^{-1} \boldsymbol {AT})\boldsymbol Y = \left( \begin{array}{l} \lambda_1 &0&0 \\ 0 & \lambda_2 & 0 \\ 0 & 0 &\lambda_3 \end{array} \right) \boldsymbol Y,\]

其中$\lambda_3 < \lambda_2 < \lambda_1 <0$。系統的全部解都趨向原點,這樣咱們就獲得了一個高維匯點(見圖6.2)。當初值條件$(x_0,y_0,z_0)$中的3個座標都非零時,對應的解將沿切於$x$軸(特徵值絕對值小的對應座標)的方向趨於原點。

image

       如今,假設$n \times n$矩陣$\boldsymbol A$具備$n$個不一樣特徵值,其中$k_1$個實的,$2k_2$個非實的,於是$n = k_1+2k_2$。因而,由第5章可知,存在座標變換使得系統化成

\[\begin{align} x'_j & =\lambda_jx_j \\ u'_l & = \alpha_lu_l + \beta_lv_l \\ v'_l & =-\beta_lu_l + \alpha_lv_l \end{align},\]

這裏$j=1,\cdots,k_1,l_1=1,\dots,k_2$。由第3章可知,咱們有以下的解:

\[\begin{align} x_j(t) & =c_je^{\lambda_jt} \\ u_l(t) & =p_l e^{\alpha_lt}\cos\beta_lt + q_l e^{\alpha_lt}\sin\beta_lt \\ v_l(t) & =-p_l e^{\alpha_lt}\sin\beta_lt + q_l e^{\alpha_lt}\cos\beta_lt. \end{align}\]

像前面同樣,咱們能夠直接驗證這就是通解。因而,咱們就證實了:

定理  考慮系統$\boldsymbol X' = \boldsymbol {AX}$,其中$\boldsymbol A$具備不一樣的特徵值,$\lambda_1,\cdots,\lambda_{k_1} \in \mathbb R$,而$\alpha_1+\text i\beta_1,\cdots,\alpha_{k_2}+\text i\beta_{k_2} \in \mathbb C$。假設$\boldsymbol T$將$\boldsymbol A$化爲標準形

\[\boldsymbol T^{-1}\boldsymbol {AT} = \left( \begin{array}{l} \lambda_1 \\&\ddots \\ && \lambda_{k_1} \\ &&& B_1 \\ &&&& \ddots \\ &&&&& B_{k_2} \end{array} \right),\]

其中

\[\boldsymbol B_j = \left( \begin{array}{l} \alpha_j & \beta_j \\ -\beta_j & \alpha_j \end{array} \right).\]

則$\boldsymbol X' = \boldsymbol {AX}$的通解爲$\boldsymbol {TY}(t)$,而

 

\[\boldsymbol Y(t) = \begin{pmatrix} {c_1}{e^{{\lambda _1}t}}\\ \vdots \\{c_{{k_1}}}{e^{{\lambda _{{k_1}}}t}}\\{a_1}{e^{{\alpha _1}t}}\cos {\beta _1}t + {b_1}{e^{{\alpha _1}t}}\sin {\beta _1}t\\ - {a_1}{e^{{\alpha _1}t}}\sin {\beta _1}t + {b_1}{e^{{\alpha _1}t}}\cos {\beta _1}t\\ \vdots \\{a_{{k_2}}}{e^{{\alpha _{{k_2}}}t}}\cos {\beta _{{k_2}}}t + {b_{{k_2}}}{e^{{\alpha _{{k_2}}}t}}\sin {\beta _{{k_2}}}t\\ - {a_{{k_2}}}{e^{{\alpha _{{k_2}}}t}}\sin {\beta _{{k_2}}}t + {b_{{k_2}}}{e^{{\alpha _{{k_2}}}t}}\cos {\beta _{{k_2}}}t \end{pmatrix}.\]

與一般同樣,矩陣$\boldsymbol T$的列向量由對應於每一個特徵值的特徵向量(或特徵向量的實部和虛部)構成。依然與前面同樣,具備負實部(正實部)的特徵值對應的特徵向量張成的子空間爲穩定(不穩定)子空間。

  考慮系統

\[\boldsymbol X' = \begin{pmatrix} 1&0&1\\ -1&0&0 \\ 0&0&-1 \end{pmatrix} \boldsymbol X,\]

這裏的矩陣已是標準形,其特徵值爲$\pm \text i,-1$。知足初值條件$(x_0,y_0,z_0)$的解爲

\[\boldsymbol Y(t) = x_0 \begin{pmatrix} \cos t \\ -\sin t \\ 0 \end{pmatrix} + y_0 \begin{pmatrix} \sin t \\ -\cos t \\ 0 \end{pmatrix} + z_0e^{-t} \begin{pmatrix} 0 \\ 0 \\ 1 \end{pmatrix},\]

這就是系統的通解。系統的相圖如圖6.3所示。穩定線爲沿$z$軸的直線,而$xy$平面上的全部解都繞中心爲原點的圓周旋轉。事實上,其它不在穩定線上的每個解都位於由「$x^2 + y^2$ =常數 」所肯定的$\mathbb R^3$中的圓柱面上,當初值條件$(x_0,y_0,z_0)$中的$z_0 \ne 0$時,這些解都盤旋地趨向$xy$平面上以$\sqrt {x_0^2 + y_0^2}$爲半徑的圓形解。

image 

  如今考慮系統 $\boldsymbol X' = \boldsymbol {AX}$,其中

\[\boldsymbol A = \begin{pmatrix} -0.1 & 0 &1 \\ -1 & 1 & –1.1 \\ -1 & 0 & –0.1 \end{pmatrix}.\]

顯然它的特徵方程爲

\[-\lambda^3 + 0.8\lambda^2 -0.81\lambda + 1.01 = 0,\]

求得特徵值分別爲$1,-0.1\pm \text i$。求解方程$(\boldsymbol A - (0.1 + \text i) \boldsymbol I ) \boldsymbol X = \boldsymbol 0$可得屬於$-0.1 + \text i$的特徵向量$(-\text i,1,1)$。記$\boldsymbol V_1 = \text {Re} (-\text i,1,1) = (0,1,1), \boldsymbol V_2 = \text {Im} (-\text i,1,1) = (-1,0,0)$。再求解方程$(\boldsymbol A - \boldsymbol I)\boldsymbol X = \boldsymbol 0$可得屬於$\lambda = 1$的特徵向量$\boldsymbol V_3 = (0,1,0)$。從而以$\boldsymbol V_i$爲列向量的矩陣

\[\boldsymbol T = \begin{pmatrix} 0&-1&0 \\ 1&0&1 \\ 1 & 0 &0  \end{pmatrix}\]

就將 $\boldsymbol X' = \boldsymbol {AX}$化爲標準形

\[\boldsymbol Y' = \begin{pmatrix} –0.1 & 1 & 0 \\ -1 & –0.1 & 0 \\ 0 & 0 & 1 \end{pmatrix} \boldsymbol Y.\]

這個系統的不穩定線爲$z$軸,而$xy$平面爲穩定平面。注意,穩定平面上的全部解都盤旋進入原點。咱們稱這個系統爲一個螺線鞍點(見圖6.4)。穩定平面外的典型解都將盤旋趨向於$z$軸,而其$z$座標增長或減少(原文和譯文彷佛都有點問題?)(見圖6.5)。

 image

image

 

 

6.2 調和振子

      考慮一對無阻尼的調和振子,它們的方程是:

\[\begin{array}{l} x''_1 = -\omega_1^2x_1 \\ x''_2 = -\omega_2^2x_2. \end{array} \]

經過檢驗腦海中的$\sin \omega t,\cos \omega t$這樣的函數,咱們就幾乎能夠解出這兩個方程。然而,讓咱們再進一步,首先對上一節中的定理在非實特徵值情形作些解釋,但更重要的是引入一些有趣的幾何。

      先引入兩個新變量$y_j=x'_j,j=1,2$,這樣方程組就能夠改寫成下面的系統

\[\begin{align} x'_1 &= y_j \\ y'_j &= -\omega_j^2x_j. \end{align} \]

該系統的矩陣形式爲$\boldsymbol X' = \boldsymbol {AX}$,其中$\boldsymbol X = (x_1,y_1,x_2,y_2)$,

\[\boldsymbol A = \begin{pmatrix} 0&1 && \\ -\omega_1^2 & 0 && \\ && 0 & 1 \\ && -\omega_2^2 & 0 \end{pmatrix}.\]

這個系統的特徵值爲$\pm \text i \omega_1, \pm \text i \omega_2$。屬於$ \text i \omega_1, \text i \omega_2$的特徵向量分別爲$\boldsymbol V_1 = (1,\text i \omega_1,0,0),\boldsymbol V_2 = (0,0,1,\text i \omega_2)$。記$\boldsymbol W_1,\boldsymbol W_2$分別爲$\boldsymbol V_1$的實部和虛部,$\boldsymbol W_3,\boldsymbol W_4$分別爲$\boldsymbol V_2$的實部和虛部。像一般業樣,令$\boldsymbol {TE}_j = \boldsymbol W_j$,則線性映射$\boldsymbol T$就將系統化爲標準形,對就的矩陣爲

\[\boldsymbol T_{-1}\boldsymbol {AT} = \begin{pmatrix} 0&\omega_1 && \\ -\omega_1 & 0 && \\ && 0 & \omega_2 \\ && -\omega_2 & 0 \end{pmatrix}.\]

這樣就能夠獲得通解(見上一節的定理)

\[\boldsymbol Y(t) = \begin{pmatrix} x_1(t) \\ y_1(t) \\ x_2(t) \\ y_2(t) \end{pmatrix} = \begin{pmatrix} a_1\cos\omega_1t+b_1\sin\omega_1t \\ -a_1\sin\omega_1t+b_1\cos\omega_1t \\ a_2\cos\omega_2t+b_2\sin\omega_2t \\ -a_2\sin\omega_2t+b_2\cos\omega_2t \end{pmatrix}, \]

它的形式正像咱們所想像的那樣。

      咱們也許會認爲事情就至此爲止了,由於解的表達式都已經獲得了,可是,咱們將往前進一步。

      顯然,每對解$(x_j(t),y_j(t)),j=1,2$都是單個方程的一個以$2\pi / \omega_j$爲週期的週期解,但這並不意味着完整的四維解也是一個週期解。事實上,完整的解是一個週期爲$\tau$的周期函數當且僅當存在整數$m,n$使得

\[\omega_1 \tau = m·2\pi, \;\;\;\; \omega_2 \tau = n·2\pi .\]

從而,爲了獲得週期性,咱們必須有

\[\tau = \frac{2\pi m}{\omega_1} =\frac {2\pi n}{\omega_2}, \]

這等價於

\[\frac{\omega_2}{\omega_1} = \frac{n}{m}. \]

即,調和振子的兩個頻率比必須是一個有理數。在圖6.6中,咱們畫出了當頻率比爲5/2時系統的一個特解$(x_1(t),x_2(t))$(這是一個平面相圖,能夠這麼想像,有一個沙漏沿$x_1$軸來回擺動,同時沿$x_2$軸來回擺動,合成的相圖即如圖所示)。

image

      當頻率之比爲無理數時,狀況很是不同。爲了理解這點,咱們做另一個(你們很是熟悉的)座標變換。用標準形來寫,如今的系統是

\[\begin{align} x'_j &= \omega_j y_j \\ y'_j &= -\omega_jx_j \end{align}. \]

咱們引入極座標$(r_j,\theta_j)$來代替變量$x_j,y_j$。對式子

\[r_j^2 = x_j^2 + y_j^2\]

兩邊求導得得

\[2r_j r'_j = 2x_jx'_j + 2y_jy'_j = 2x_jy_j\omega_j - 2x_jy_j\omega_j = 0.\]

因而,對每個$j$都有$r'_j = 0$。一樣地,對方程

\[\tan \theta_j  = \frac{y_j}{x_j}\]

求導可得

\[(\sec^2 \theta_j) \theta'_j = \frac{y'_jx_j - y_jx'_j}{x_j^2} = \frac{-\omega_jr_j^2}{r_j^2\cos^2 \theta_j},\]

由此可得

\[\theta'_j = -\omega_j.\]

從而在極座標下,這些方程就變得很是簡單:

\[\begin{align} r'_j & =0 \\ \theta'_j &= -\omega_j \end{align}. \]

第一個方程告訴咱們,$r_1,r_2$沿任何解都始終爲常數。從而,不管咱們如何選取$r_1,r_2$的初始值,關於$\theta_j$的方程都是同樣的。這樣咱們就能夠將注意力集中到$r_1 = r_2 =1$的情形。在$\mathbb R^4$中獲得的點集是一個環面,也就是油炸圈餅的表面(挺形象了,試想一隻螞蟻既能夠沿橫截面的一個圓做週期運動,也能夠沿縱截面的一個圓做週期運動)。雖然在四維空間中這有點難以想像,但無論怎樣,這個集合上有2個獨立的變量,即$\theta_1,\theta_2$,並且它們都以$2\pi$爲週期的。在這一點上,這相似於$\mathbb R^3$中熟知的環面能夠被2個獨立的循環走向所參數化(螞蟻能夠走到油炸圈餅的任何一點)。

      如今,方程限制在環面上就變成了:

\[\begin{align}\theta'_1  &= -\omega_1 \\ \theta'_2 & = -\omega_2. \end{align} \]

爲了方便,咱們可將$\theta_1,\theta_2$當作是邊長爲$2\pi$的正方形中的變量(只要將正方形的對邊$\theta_j = 0$和$\theta_j = 2\pi$粘合起來就獲得了環面——注:這已是一種很好的處理方式了,儘管很難將一個正方形變成四維中的環)。在正方形中,向量場的斜率爲常數

\[\frac{\theta'_2}{\theta'_1} = \frac{\omega_2}{\omega_1}.\]

因而,在正方形中,解都位於斜率爲$\omega_2 / \omega_1$的直線上。當解到達$\theta_1 = 2\pi$這條邊時(不妨設在$\theta_2 = c$處),它就在$\theta_1 = 0$的邊立刻從新出現(($\theta_2$的座標依然爲$c$),而後以斜率$\omega_2 / \omega_1$再繼續向前。在解到達$\theta_2 = 2\pi$時,將會出現相似的等同。

      這樣,咱們就對解在環面上的行爲有了一個簡化的幾何圖像(何時,國內的書能寫得如此簡明而透徹呢?)。將會出現什麼呢?答案依賴於$\omega_2 / \omega_1$的取值。若是這個值爲有理數,好比,$n/m$,則解從$(\theta_1(0),\theta_2(0))$出發將垂直通過環面$n$次、水平通過環面$m$次以後再回到出發點。這就是咱們上面觀察到的週期解。有時,直線解在$\theta_1 \theta_2$平面上的圖像和圖6.6中所畫出的$x_1x_2$平面上的圖像並不徹底同樣。

      在無理情形,將會出現一些很不一樣的狀況(見圖6.7)。爲了理解如今出現的狀況,讓咱們回到第1章討論過的龐加萊映射的概念。考慮圓周$\theta_1 = 0$,也就是環面正方形表示中的左邊。任給該圓周上的一個初值,好比,$\theta_2 = x_0$,咱們跟隨從該點出發的解一直到它下次再到達$\theta_1 = 2\pi$。因爲環面是等同正方形的對邊獲得的,解如今實際上就是回到了出發的圓周$\theta_1 = 0$。在這個進行過程當中,解有可能會穿過邊界$\theta_2 = 2\pi$幾回,但它最終會回到$\theta_1 = 0$。於是咱們能夠定義$\theta_1 = 0$上的龐加萊映射以下:龐加萊映射在圓周上$x_0$處的取值就是首次返回點的相應座標。假設首次返回出如今點$\theta_2(\tau)$處,其中$\tau$是知足$\theta_1(\tau) = 2\pi$的時刻。因爲$\theta_1(t) = \theta_1(0) + \omega_1t$,於是$\tau = 2\pi / \omega_1$。從而$\theta_2(t) = x_0 + \omega_2(2\pi / \omega_1)$。因而在這個圓周上的龐加萊映射能夠寫成

\[f({x_0}) = {x_0} + 2\pi ({\omega _2}/{\omega _1})\;\bmod \;2\pi \]

其中$x_0 = \theta_2(0)$是初值在這個圓周上的$\theta_2$座標(見圖6.8)(上式可理解爲初始位置在$x_0$,步長爲$2\pi ( \omega_2 / \omega_1)$的一種運動,只是要進行模$2\pi$,並依次記錄點的軌跡,所以也是一種迭代過程)。因而,該圓周上的龐加萊映射就是將圓周上的點旋轉角度$2\pi( \omega_2 / \omega_1)$這個函數。因爲$\omega_2 / \omega_1$是無理的,這個函數稱爲圓周上的無理旋轉

image

image

定義  點集$x_0,x_1 = f(x_0),x_2 = f(f(x_0)),\cdots, x_n = f(x_{n-1})$稱爲$x_0$在$f$迭代下的軌道

      $x_0$的軌道記錄了當時間增長時,解是如何穿過$\theta_1 = 2\pi$的(言下之意是另外一個變量的狀態是什麼樣的)。

命題  假設$\omega_2 / \omega_1$是無理數,則圓周$\theta_1 = 0$上任何初值$x_0$的軌道在該圓周上都是稠密的(意思是幾乎充滿整個圓周)。

證實詳見書本。

      由於$x_0$的軌道在圓周$\theta_1 = 0$上是稠密的,從而在正方形中將這些點鏈接起來的直線解也是稠密的,因而原來的解在其所在的環面上不稠密的。這就解釋了爲何當$\omega_2 / \omega_1 = \sqrt2$時,解如圖6.7所示到$x_1x_2$平面時會稠密地佈滿。

      回到振子的實際運動,咱們就獲得當$\omega_2 / \omega_1$爲無理數時,兩個質點的運動並非週期的。然而,因爲解在環面上的稠密性,當時間增長時,它們又的確會一次次地回到很是靠近初始位置的地方。這種形式的運動稱爲擬週期運動

 

6.3 重特徵值

      在上一章中,咱們已經看到,具備實重特徵值系統的求解歸結爲求解其矩陣含有以下形式的塊的系統(若是存在不惟一解的重特徵值呢?):

\[\left( \begin{array}{l} \lambda & 1 \\ & \lambda & 1 \\ & & \ddots \ddots \\ & & & \ddots &1 \\ & & & & \lambda \end{array} \right). \]

  令

\[\boldsymbol X' = \begin{pmatrix} \lambda & 1 & 0 \\ 0 & \lambda & 1 \\ 0 & 0 & \lambda \end{pmatrix} \boldsymbol X.\]

這個系統惟一的特徵值就是$\lambda$,並且惟一的特徵向量(1,0,0)。咱們像在第3章所作的那樣來求解這個系統。首先,因爲$x'_3 = \lambda x_3$,故必有

\[x_3(t) = c_3e^{\lambda t},\]

進而有,

\[x'_2 = \lambda x_2 + c_3e^{\lambda t}.\]

與第3章同樣,咱們猜想它有以下形式的解:

\[x_2(t) = c_2e^{\lambda t}+ \alpha t e^{\lambda t}.\]

將上式代入$x'_2$的微分方程,能夠定出$\alpha = c_3$,從而找到

\[x_2(t) = c_2e^{\lambda t}+ c_3 t e^{\lambda t}.\]

最後,方程

\[x'_1 = \lambda x_1 + c_2e^{\lambda t} + c_3 t e^{\lambda t}\]

讓人聯想起下面形式的解:

\[x_1(t) = c_1 e^{\lambda t} + \alpha t e^{\lambda t} + \beta t^2 e^{\lambda t}. \]

像剛纔同樣的求解可得

\[x_1(t) = c_1 e^{\lambda t} + c_2 t e^{\lambda t} +  c_3 \frac {t^2}{2} e^{\lambda t}. \]

總之,咱們找到了

\[\boldsymbol X(t) = c_1 e^{\lambda t} \begin{pmatrix} 1 \\ 0 \\ 0 \end{pmatrix} + c_2 e^{\lambda t} \begin{pmatrix} t \\ 1 \\ 0 \end{pmatrix} + c_3 e^{\lambda t} \begin{pmatrix} t^2/2 \\ t \\ 1 \end{pmatrix} \]

它就是系統的通解。雖然解的表達式中出現了多項式項,但當$\lambda < 0$時,指數將起主導做用,並且全部的解都趨於零。當$\lambda < 0$時,一些有表明性的解如圖6.9所示。注意,該系統只有惟一的直線解(關於直線解的概念,能夠參考2.3節和2.6節,對於實特徵值對應的特徵向量就是直線解),該解位於$x$軸上,並且$xy$平面是不變的,其上解的形態與平面重特徵值情形徹底相同。

image

  考慮下面的四維繫統:

\[\begin{align} x'_1 &= x_1 + x_2 - x_3 \\ x'_2 &= x_2 + x_4 \\ x'_3 &= x_3 + x_4 \\ x'_4 &= x_4 \end{align} \]

咱們能夠將該系統寫成以下的矩陣形式:

\[\boldsymbol X' = \boldsymbol {AX} = \begin{pmatrix} 1&1&-1&0 \\ 0&1&0&1 \\ 0&0&1&1 \\ 0&0&0&1 \end{pmatrix} \boldsymbol X.\]

由於$\boldsymbol A$爲上三角的,故它的全部特徵值都是1。求解$(\boldsymbol A - \boldsymbol I) \boldsymbol X = \boldsymbol 0$可找到兩個線性無關的特徵向量$\boldsymbol V_1 = (1,0,0,0),\boldsymbol W_1 = (0,1,1,0)$。這就將$\boldsymbol A$可能的標準形減小到兩種。進一步求解$(\boldsymbol A - \boldsymbol I) \boldsymbol X = \boldsymbol V_1$可得一個解$\boldsymbol V_2 = (0,1,0,0)$,求解$(\boldsymbol A - \boldsymbol I) \boldsymbol X = \boldsymbol W_1$可得一個解$\boldsymbol W_2 = (0,0,0,1)$(這種方法的理論可參考5.5節的內容)。由此,咱們知道系統$\boldsymbol X' = \boldsymbol {AX}$能夠變換成

\[\boldsymbol Y' = (\boldsymbol T^{-1} \boldsymbol {AT}) \boldsymbol Y = \begin{pmatrix} 1&1&0&0 \\ 0&1&0&0 \\ 0&0&1&1 \\ 0&0&0&1 \end{pmatrix} \boldsymbol Y,\]

其中矩陣$\boldsymbol T$由下式給出:

\[\boldsymbol T = \begin{pmatrix} 1&0&0&0 \\ 0&1&1&0 \\ 0&0&1&0 \\ 0&0&0&1 \end{pmatrix} \]

因而$\boldsymbol Y' = (\boldsymbol T^{-1} \boldsymbol {AT}) \boldsymbol Y$的解爲(重特徵值情形下求通解可參考3.3節的內容 )

\[\begin{align} y_1(t) &= c_1 e^t + c_2 t e^t \\ y_2(t) &= c_2 e^t \\ y_3(t) &= c_3 e^t + c_4 t e^t  \\ y_4(t) &= c_4 e^t. \end{align} \]

再用座標變換$\boldsymbol T$做用,咱們就獲得了原系統的通解爲

\[\begin{align} x_1(t) &= c_1 e^t + c_2 t e^t \\ x_2(t) &= c_2 e^t + c_3 e^t + c_4 t e^t \\ x_3(t) &= c_3 e^t + c_4 t e^t  \\ x_4(t) &= c_4 e^t. \end{align} \]

 

6.4 矩陣指數

       如今,咱們轉到矩陣指數求解線性系統,這是解線性系統的一個不一樣但簡練的方法。從某種意義上講,這是攻克這類系統天然的方法。

      咱們先回憶1.1節中是如何求解線性方程$x' = ax$的$1 \times 1$「系統」,其中咱們的矩陣就是簡單的$(a)$。那時咱們並無進行找特徵值和特徵向量的過程,(然而,事實上,咱們進行了,但過程實在是太簡單了。)而是直接取矩陣$(a)$的指數函數就找到了通解$x(t)  = c \text {exp} (at)$。事實上,這個過程對通常$n \times n$矩陣$\boldsymbol A$也是行得通的。咱們惟一須要知道的就是如何求一個矩陣的指數。

      辦法以下。由微積分可知,指數函數能夠表示成一個無窮級數

\[e^x = \sum \limits_{k=0}^\infty \frac {x^k}{k!},\]

並且咱們知道對每一個$x \in \mathbb R$這個級數都收斂。因爲咱們能夠將矩陣相加,能夠取它們的$k$次方,還能夠將矩陣的每一個元素乘上$1/k!$,這將意味着咱們同樣能夠用上面的級數去取它們的指數。

定義  設$\boldsymbol A$爲一個$n \times n$矩陣。咱們定義$\boldsymbol A$的指數是由下式給出的矩陣:

\[\text {exp} (\boldsymbol A) = \sum \limits_{k=0}^\infty \frac {\boldsymbol A^k}{k!}.\]

      固然,咱們還要擔憂上面矩陣和式收斂的含義。

  令

\[\boldsymbol A = \begin{pmatrix} \lambda & 0 \\ 0 & \mu \end{pmatrix} . \]

則咱們有

\[\boldsymbol A^k =\begin{pmatrix} \lambda^k & 0 \\ 0 & \mu^k \end{pmatrix},\]

因而

\[\text {exp}(\boldsymbol A) = \begin{pmatrix} \sum \limits_{k=0}^ \infty \lambda^k / k!  & 0 \\ 0 & \sum \limits_{k=0}^\infty \mu^k / k! \end{pmatrix}  = \begin{pmatrix} e^\lambda & 0 \\ 0 & e^\mu \end{pmatrix},\]

這個結果你們可能已經事先猜到。

  下面來看一個稍等複雜些的例子。令

\[\boldsymbol A = \begin{pmatrix} 0 & \beta \\ -\beta & 0 \end{pmatrix}.\]

由計算可得

\[\text {exp} (\boldsymbol A) =\begin{pmatrix} \sum \limits_{k=0}^\infty (-1)^k \frac{\beta^{2k}}{(2k)!} & \sum \limits_{k=0}^\infty (-1)^k \frac{\beta^{2k+1}}{(2k+1)!} \\ -\sum \limits_{k=0}^\infty (-1)^k \frac{\beta^{2k+1}}{(2k+1)!} & \sum \limits_{k=0}^\infty (-1)^k \frac{\beta^{2k}}{(2k)!} \end{pmatrix} = \begin{pmatrix} \cos \beta & \sin \beta \\ - \sin\beta & \cos\beta \end{pmatrix}.\]

  如今令

\[\boldsymbol A = \begin{pmatrix} \lambda & 1 \\ 0 & \lambda \end{pmatrix},\]

其中$\lambda \ne 0$。爲了照顧隨後的應用,咱們來計算$\text {exp} (t \boldsymbol A)$,而不是$\text {exp} A)$。咱們有

\[(t \boldsymbol A)^k = \begin{pmatrix} (t\lambda)^k & kt^k\lambda^{k-1} \\ 0 & (t\lambda)^k \end{pmatrix}.\]

因而咱們獲得

\[\text {exp}(t \boldsymbol A)  = \begin{pmatrix} \sum \limits_{k=0}^\infty \frac{(t\lambda)^k}{k!} & t\sum \limits_{k=0}^\infty \frac{(t\lambda)^k}{k!}  \\ 0 & \sum \limits_{k=0}^\infty \frac{(t\lambda)^k}{k!} \end{pmatrix}= \begin{pmatrix} e^{t\lambda} & t e^{t\lambda} \\ 0 & e^{t\lambda}.\end{pmatrix} \]

      能夠看到,在上面3個例子中,矩陣$\text {exp}(\boldsymbol A)$的元素都是一些無窮級數。所以,咱們稱矩陣無窮級數$\text {exp}(\boldsymbol A)$絕對收斂,若是每個單獨的項絕對收斂。在上面的幾種情形中,收斂性是顯然的。遺憾的是,對於通常的矩陣$\boldsymbol A$,收斂性並不清楚。爲了證實這裏的收斂性,咱們須要作一些努力。

      記$a_{ij}(k)$爲$\boldsymbol A^k$的$ij$元素。令$a = \text {max}|a_{kj}|$。咱們有

\[|a_{ij}(k)| \le n^{k-1}a^k. \]

因而,咱們獲得了$n \times n$矩陣$\text {exp}(\boldsymbol A)$的$ij$元素的一個上界:

\[\Bigg| \sum \limits_{k=0}^\infty \frac{a_{ij}(k)}{k!} \Bigg| \le \sum \limits_{k=0}^\infty \frac{|a_{ij}(k)|}{k!} \le \sum \limits_{k=0}^\infty \frac{n^{k-1}a^k}{k!} \le \text {exp}(na),\]

所以,根據比較判別法,這個級數絕對收斂。這樣,對任何$\boldsymbol A \in L(\mathbb R^n),\text {exp}(\boldsymbol A)$有意義。

      下面的結果代表,矩陣指數也具備一般指數函數的許多熟悉的性質。

命題  設$\boldsymbol A,\boldsymbol B,\boldsymbol T$都是$n \times n$矩陣,則
      (1)若是$\boldsymbol B=\boldsymbol T^{-1} \boldsymbol {AT}$,則$\text {exp}(\boldsymbol B) = \boldsymbol T^{-1} \text {exp}(\boldsymbol A) \boldsymbol T$。
      (2)若是$\boldsymbol {AB} = \boldsymbol {BA}$,則$\text {exp}(\boldsymbol A+\boldsymbol B) = \text {exp} (\boldsymbol A)\text {exp} (\boldsymbol B)$。
      (3)$\text {exp}(- \boldsymbol A) = (\text {exp} (\boldsymbol A))^{-1}$。

證實詳見書本。注意,命題中論斷(3)意味着,對每個矩陣$\boldsymbol A,\text {exp}(\boldsymbol A)$均可逆。這一點相似於對每一個實數$a,e^a \ne 0$。

      在$\boldsymbol A$和$\text {exp}(\boldsymbol A)$的特徵向量之間有一個很簡單的關係:

命題  若是$\boldsymbol V \in \mathbb R^n$是$\boldsymbol A$的屬於特徵值$\lambda$的特徵向量,則$\boldsymbol V$也是$\text {exp} (\boldsymbol A)$的特徵向量,它屬於$e^\lambda$。

       如今咱們回到微分方程系統的討論。令$\boldsymbol A$爲$n \times n$矩陣,考慮系統$\boldsymbol X' = \boldsymbol {AX}$。回想咱們之前用$L(\mathbb R^n)$記全體$n \times n$矩陣的集合。把每一個$t \in \mathbb R$映到$\text {exp}(t\boldsymbol A)$就獲得了一個函數$\mathbb R \to L(\mathbb R^n)$。因爲咱們將$L(\mathbb R^n)$等同於$\mathbb R^{n^2}$,於是談論這個函數的導數是有意義的。

命題

\[\frac{\text d}{\text d t} \text {exp}(t \boldsymbol A) = \boldsymbol A \text {exp}(t \boldsymbol A) = \text {exp}(t \boldsymbol A) \boldsymbol A. \]

換句話說,矩陣值函數$t \to \text {exp}(t \boldsymbol A)$的導數是另外一個矩陣值函數$\boldsymbol A \text {exp} (t \boldsymbol A)$。

      如今咱們回到求解微分方程系統。下面的結論能夠當作是具備常係數的線性微分方程組的基本定理。

定理  設$\boldsymbol A$爲$n \times n$矩陣,則初值問題$\boldsymbol X' = \boldsymbol {AX},\boldsymbol X(0) = \boldsymbol X_0$的解爲$\boldsymbol X(t) = \text {exp}(t \boldsymbol A)\boldsymbol X_0$,並且它是惟一解。

  考慮系統

\[\boldsymbol X' = \begin{pmatrix} \lambda &1 \\ 0 & \lambda \end{pmatrix} \boldsymbol X.\]

根據上面的定理,它的通解爲

\[\boldsymbol X(t) = \text {exp}(t \boldsymbol A) \boldsymbol X_0 = \text {exp} \begin{pmatrix} t \lambda & t \\ 0 & t \lambda \end{pmatrix} \boldsymbol X_0.  \]

這其中的矩陣指數在前面已經算過,因而咱們有

\[\boldsymbol X(t) = \begin{pmatrix} e^{t\lambda} & t e^{t\lambda} \\ 0 & e^{t\lambda}\end{pmatrix} \boldsymbol X_0.\]

能夠看出,這與第3章中的計算結果相吻合。

 

6.5 非自治線性系統

      直到如今,幾乎全部碰到的線性微分方程系統都是自治的。可是,在應用中仍是常常會出現一些形式的非自治系統。其中之一就是以下形式的系統

\[\boldsymbol X' = \boldsymbol A(t) \boldsymbol X,\] 

其中$\boldsymbol A(t) = [a_{ij}(t)]$爲$n \times n$矩陣,它連續地依賴於時間。當之後章節中遇到變分方程時,咱們會再深刻研究這類系統。

      如今咱們只關係另外一個不一樣的非自治線性系統:

\[\boldsymbol X' = \boldsymbol {AX} + \boldsymbol G(t),\] 

其中$\boldsymbol A$爲常值的$n \times n$矩陣,$\boldsymbol G:\mathbb R \to \mathbb R^n$顯式地依賴於$t$,咱們稱之爲強迫項。這是一個一階線性非自治方程系統。

  (受迫調和振子)若是調和振子系統中有一個外力做用,決定運動的微分方程就變成

\[x'' + bx' + kx = f(t),\]

其中$f(t)$爲外力的大小。一個重要的特殊情形是其中的外力是時間的周期函數,這可對應於,比方說,將質點-彈簧裝置所在的桌子先後週期地移動。做爲一個系統,受迫調和振子方程變成

\[\boldsymbol X' = \begin{pmatrix} 0 & 1 \\ -k & -b \end{pmatrix} \boldsymbol X + \boldsymbol G(t),\;\;\; \boldsymbol G(t) = \begin{pmatrix} 0 \\ f(t) \end{pmatrix}.\]

      對一個非齊次的系統,若是將其中的時間項扔掉後,所得的系統$\boldsymbol X' = \boldsymbol {AX}$稱爲齊次方程。咱們已經知道如何找到這類系統的通解。利用上節中的記號,齊次系統知足初值條件$\boldsymbol X(0)= \boldsymbol X_0$的解爲

\[\boldsymbol X(t) = \text {exp}(t \boldsymbol A) \boldsymbol X_0,\]

而這也是齊次方程的通解。

      爲了找到非齊次方程的通解,假設咱們已經有了該方程的一個特解$\boldsymbol Z(t)$。因而$\boldsymbol Z'(t) = \boldsymbol {AZ}(t) + \boldsymbol G(t)$。當$\boldsymbol X(t)$爲齊次方程的任一個解時,則函數$\boldsymbol Y(t) = \boldsymbol X(t) + \boldsymbol Z(t)$就是非齊次方程的另外一個解。這點可由以下的計算得出:

\[\begin{align} \boldsymbol Y' = \boldsymbol X' + \boldsymbol Z' &=\boldsymbol {AX} + \boldsymbol {AZ} +\boldsymbol G(t) \\ &=\boldsymbol A(\boldsymbol X+\boldsymbol Z) +\boldsymbol G(t) \\ &=\boldsymbol {AY}+\boldsymbol G(t).\end{align} \]

於是,因爲已經了齊次方程的全部解,一旦找到了非齊次方程的哪怕只是一個特解,咱們就立刻能夠找到非齊次方程的通解。一般,咱們能夠經過猜解獲得這樣一個解(微積分中,這種方法一般稱爲待定係數法)。然而 ,猜解並不老是可行。下面的參數變易法卻老是可行的,固然,這並不能保證咱們總能計算其中的積分。

定理   (參數變易法)考慮非齊次方程

\[\boldsymbol X' = \boldsymbol {AX} + \boldsymbol G(t),\] 

其中$\boldsymbol A$爲常值的$n \times n$矩陣,$\boldsymbol G(t)$爲$t$的連續函數。則

\[\boldsymbol X(t) = \exp (t \boldsymbol A)\left( {{\boldsymbol X_0} + \int_0^t {\exp ( - s\boldsymbol A) \boldsymbol G(s) \text ds} } \right)\]

是該方程知足$\boldsymbol X(0) = X_0$的一個解。(對方程兩邊求導便可驗證)

      下面咱們來給出這個結果在週期受迫調和振子中的幾個應用。首先,假設有個外力爲$\cos t$的阻尼振子,此時外力項的週期爲$2\pi$。這個系統爲

\[\boldsymbol X' = \boldsymbol {AX} + \boldsymbol G(t),\]  

其中$\boldsymbol G(t) = (0,\cos t), \boldsymbol A$爲矩陣

\[\boldsymbol A = \begin{pmatrix} 0&1 \\ -k & -b \end{pmatrix},\]

其中$b,k>0$。咱們斷言:這個系統有惟一一個以$2\pi$爲週期的週期解。爲了證實這個斷言,咱們首先須要找到一個知足條件$\boldsymbol X(0) = \boldsymbol X_0 =\boldsymbol X(2\pi)$的解。由參數變易法,咱們須要找到$\boldsymbol X_)$使得

\[{\boldsymbol X_0} = \exp (2\pi \boldsymbol A){\boldsymbol X_0} + \exp (2\pi \boldsymbol A)\int_0^{2\pi } {\exp ( - s\boldsymbol A)\boldsymbol G(s)\text ds.} \]

如今其中一項

\[\exp (2\pi \boldsymbol A)\int_0^{2\pi }  \exp ( - s\boldsymbol A)\boldsymbol G(s)\text ds\]

爲常值向量,將它記爲$\boldsymbol W$。這樣咱們就要求解方程

\[(\exp (2\pi \boldsymbol A) - \boldsymbol I)\boldsymbol X_0 = -\boldsymbol W.\]

因爲$\exp (2\pi \boldsymbol A) - \boldsymbol I$可逆,這個方程有惟一解。爲了證實這個矩陣的可逆性,假設它是不可逆的,則有非零向量$\boldsymbol V$使得

\[(\exp (2\pi \boldsymbol A) - \boldsymbol I) \boldsymbol V = \boldsymbol 0,\]

或者,換句話說,矩陣$\exp (2\pi \boldsymbol A)$有一個特徵值爲1。可是,由上節可知,$\exp (2\pi \boldsymbol A)$的特徵值都由$\exp (2\pi \lambda_j)$給出,其中$\lambda_j$爲$\boldsymbol A$的特徵值。因爲每一個$\lambda_j$的實部都小於0,從而$\exp (2\pi \lambda_j)$小於1,因而矩陣$\exp (2\pi \boldsymbol A) - \boldsymbol I$的確是可逆的。這樣致使$2\pi$週期解的惟一初值(這多少有點出人意料,是否是週期解不只依賴強迫項,同時還依賴初值,而且這樣的初值仍是惟一的!!)就是

\[\boldsymbol X_0 = (\exp (2\pi \boldsymbol A) - \boldsymbol I)^{-1}(-\boldsymbol W).\]

      令$\boldsymbol X(t)$爲知足$\boldsymbol X(0) = \boldsymbol X_0$的週期解。這個解稱爲穩態解。若是$\boldsymbol Y_0$是其它的初值條件,咱們可將它寫成$\boldsymbol Y_0 = (\boldsymbol Y_0 -\boldsymbol X_0) +\boldsymbol X_0$,從而過$\boldsymbol Y_0$的解爲

\[\begin{align} \boldsymbol Y(t) &= \exp (t \boldsymbol A)({ \boldsymbol Y_0} - { \boldsymbol X_0}) + \exp (t \boldsymbol A){ \boldsymbol X_0} + \exp (t \boldsymbol A)\int_0^t {\exp ( - s \boldsymbol A) \boldsymbol G(s)\text ds} &= \exp (t \boldsymbol A)({ \boldsymbol Y_0} - { \boldsymbol X_0}) + \boldsymbol X(t). \end{align}\]

對於表達式中的第一項,因爲它是齊次方程的一個解,於是在$t \to \infty$時,趨於0(只是針對本例而言,由於是帶阻尼的) 。因而,這個系統的每一個解在$t \to \infty$時都趨於穩態解。在物理上,這一點很明顯的:帶阻尼的(非受迫)調和振子的運動都趨於平衡態,剩下的運動就是由週期外力引發的。這樣,咱們就證實了:

定理   考慮受迫帶阻尼的調和振子方程

\[x'' + bx' + kx = \cos t,\]

其中$,b>0$。則這個方程的全部解都趨於一個以$2\pi$爲週期的穩態解。

      如今考慮一個受迫的無阻尼調和振子的特例:

\[\boldsymbol X' = \begin{pmatrix} 0 & 1 \\ -1 & 0 \end{pmatrix} \boldsymbol X + \begin{pmatrix} 0 \\ \cos\omega t \end{pmatrix},\]

其中,如今外力的週期爲$2\pi/\omega,\omega \ne \pm 1$。記

\[\boldsymbol X = \begin{pmatrix} 0 & 1 \\ -1 & 0 \end{pmatrix}.\]

齊次方程的解爲

\[\boldsymbol X(t) = \exp (t \boldsymbol A) \boldsymbol X_0 = \begin{pmatrix} \cos t & \sin t \\ -\sin t& \cos t \end{pmatrix} \boldsymbol X_0.\]

由參數變易法可知,非齊次方程從原點出發的解爲

\[\begin{align} \boldsymbol Y(t) &= \exp (t\boldsymbol A)\int_0^t \exp ( - s\boldsymbol A) \begin{pmatrix} 0\\ \cos \omega s\end{pmatrix}\text ds \\ &= \exp (t\boldsymbol A)\int_0^t \begin{pmatrix} \cos s &- \sin s \\ \sin s &\cos s \end{pmatrix} \begin{pmatrix} 0\\ \cos \omega s \end{pmatrix} \text ds\\ &= \exp (t\boldsymbol A) \int_0^t \begin{pmatrix} - \sin s\cos \omega s \\ \cos s\cos \omega s \end{pmatrix} \text ds \\ &= \frac{1}{2} \exp (t\boldsymbol A)\int_0^t \begin{pmatrix} \sin (\omega - 1)s - \sin (\omega + 1)s \\ \cos (\omega - 1)s + \cos (\omega + 1)s \end{pmatrix} \text ds. \end{align} \]

回想到

\[\exp (t \boldsymbol A) = \begin{pmatrix} \cos t & \sin t \\ -\sin t& \cos t \end{pmatrix},\]

而後利用$\omega \ne \pm 1$的事實,算出其中的積分值再通過一個長的計算過程後就獲得

\[\begin{align} \boldsymbol Y(t) &= \frac{1}{2}\exp (t\boldsymbol A) \begin{pmatrix}\frac{{ - \cos (\omega - 1)t}}{{\omega - 1}} + \frac{{\cos (\omega + 1)t}}{{\omega + 1}}\\\frac{{sin(\omega - 1)t}}{{\omega - 1}} + \frac{{sin(\omega + 1)t}}{{\omega + 1}}\end{pmatrix} + \exp (t\boldsymbol A) \begin{pmatrix}{({\omega ^2} - 1)^{ - 1}}\\0\end{pmatrix} \\&= \frac{1}{{{\omega ^2} - 1}} \begin{pmatrix} - \cos \omega t\\\omega \sin \omega t\end{pmatrix} + \exp (t\boldsymbol A) \begin{pmatrix}{({\omega ^2} - 1)^{ - 1}}\\0\end{pmatrix} . \end{align}\]

這樣就獲得了原方程的通解爲(由於上式爲初值爲零時,原方程的通解,加上初值爲$X_0$時,相應齊次方程的通解,即獲得初值爲$X_0$的非齊次方程的通解)

\[\boldsymbol Y(t) = \exp (t\boldsymbol A) \left(\boldsymbol X_0+ \begin{pmatrix}{({\omega ^2} - 1)^{ - 1}}\\0 \end{pmatrix} \right) + \frac{1}{{{\omega ^2} - 1}} \begin{pmatrix} - \cos \omega t\\\omega \sin \omega t\end{pmatrix}. \]

這個表達式中的第一項是以$2\pi$爲週期的,而第二項是以$2\pi / \omega$爲週期的。與有阻尼的情形不一樣,如今的解不必定會有周期運動。事實上,這個解爲週期的當且僅當$\omega$是一個有理數。當$\omega$爲無理數時,這個運動是擬週期的,這正好與6.2節中看到的情形是同樣的。

小結:本章主要討論了高維線性系統解的特性,經過調和振子的例子,引入了一些新的概念,如無理旋轉,擬週期。着重討論了重特徵值情形下解的形式和性質,須要結合5.5節中相關內容進行學習,即重特徵值情形下如何求特徵向量。最後討論了一個非自治線性系統的例子,用以說明強迫項、穩態解等。

相關文章
相關標籤/搜索