简单实用！3个德国人创造的线性迭代法，超越了一个时代

加星标，才能不错过每日推送！方法见文末动图

今年12月26日是德国数学家高斯发明历史上第一个线性迭代法200周年。在前不久发表的《怎样迭代求解线性方程组？》的一文中，我们从基础概念出发，厘清了迭代法求解线性方程组的准备知识和具体过程。这次我们将探讨历史上最有名气、也最简单实用的两种一阶定常迭代法：雅可比迭代法和高斯-赛德尔迭代法的基本思想和收敛性。

撰文 | 丁玖（美国南密西西比大学数学系教授）
我在前一篇《返朴》文章《怎样迭代求解线性方程组？》中，为了能与适用于一般完备距离空间的压缩映射定理挂上钩，只给出了迭代法求解线性方程组的一个简单的收敛性充分条件，即若要迭代格式x_k = Mx_k-1 + c, k = 1, 2, 3, …对所有的初始列向量x₀都收敛，一个对迭代矩阵简单易懂的要求是：Rⁿ上的向量2-范数所诱导出的矩阵2-范数||M||₂小于1。这时，迭代向量序列{x_k}当k趋向于无穷大时将收敛到线性方程组x = Mx + c的唯一解。在上述假设下，一个直接的推论是I - M为一可逆矩阵。
然而，这里又产生一个技术问题：M的欧几里得2-范数极难计算，至少它不能从M各行各列的具体元素中立即通过简单代数运算算出。事实上，线性代数中的二次型理论告诉我们，实数方阵M的2-范数等于M的转置矩阵M^T与M的乘积M^TM这个所谓“半正定矩阵”（意指二次型x^TM^TMx = （||Mx||₂）² 对所有的实列向量x都是非负实数）的最大特征值之平方根（因为M^TM的所有特征值均为非负数，故平方根存在）。上句话里包含了好几个数学概念，可想而知计算出||M||₂绝对不是一件轻而易举的事情。

条条大路通罗马好在我们还有其他的范数可用。作为通常直线段长度概念的直接推广，n维欧几里得空间Rⁿ中向量的2-范数只是一种范数，在多元微积分里用得最多，但它是更一般的q-范数取q = 2时的特例。选取一个大于或等于1的实数q，对于分量为a₁, a₂, …, a_n的任意n维向量x，它的q-范数被定义为
||x||_q = {|a₁|^q + … + |a_n|^q}^1/q。
这些范数都满足前一篇文章《怎样迭代求解线性方程组？》里列出的欧几里得范数符合的那三项“范数公理”(i)-(iii)。这样，Rⁿ就有无穷多个q-范数，此外还有无穷多个不属于这个范畴的其他范数。如此数目众多的范数，会不会令人看得眼花缭乱？不同的范数会影响线性迭代向量序列的收敛性吗？幸好有一个关于“有穷维线性赋范空间”的基本事实让我们放心：Rⁿ上的任意两个范数||•||和||•||'都是所谓“彼此等价”的范数，即存在两个正数α和β，使得对Rⁿ中的所有向量x都有
α ||x|| ≤ ||x||' ≤ β ||x||。
范数的等价性有什么妙用呢？答案是它们各自定义出的序列收敛性彼此等价，也就是说Rⁿ中的向量序列{x_k}在由范数||•||所定义的距离下收敛于一个向量，当且仅当它在由范数||•||'所定义的距离下收敛于同一个向量。我们给出这个事实的证明，因为它不难：设当k趋向于无穷大时||x_k - x|| → 0，则||x_k - x||' ≤ β ||x_k - x|| → 0；反之，若||x_k - x||' → 0，则||x_k - x|| ≤ (1/α) ||x_k - x||' → 0。顺便提及，该性质是有穷维空间的特征，即一旦线性空间的维数变成无穷大，那么赋予它的两个范数不一定导出同等的收敛性。所以，从“有穷”到“无穷”，它们之间常有“万里之遥”。
这样一来，n维向量空间Rⁿ在随便哪个范数||•||下都变成一个完备的距离空间。给定一个n阶的方阵M，这个被选用的向量范数诱导出对应的矩阵范数||M||，它等于连续函数||Mx||在(n – 1)维单位球面上的最大值。对同一个矩阵，不同的向量范数诱导出不同的矩阵范数，因而如果某一个向量范数诱导出的矩阵范数小于1，那么根据巴拿赫不动点定理，线性迭代法x_k = Mx_k-1 + c, k = 1, 2, 3, …对任一初始点x₀都收敛到同一个极限。
有两个对于计算矩阵范数而言是最方便的向量范数，其中之一是上面的q-范数取q = 1时的1-范数，即n维向量x = (a₁, a₂, …, a_n)^T的1-范数为
||x||₁ = |a₁| + … + |a_n|。
另一个常用的向量范数是q-范数当q趋向于无穷大时的极限范数： ，称为∞-范数。不难证明，n维向量x的∞-范数为
||x||_∞ = max{|a₁|, …, |a_n|}。
可以验证，但细节比较繁琐，不在这里表出，Rⁿ上的向量1-范数所诱导出的n阶方阵M = [m_ij] 的1-范数为
||M||₁ = max {|m₁₁| + … + |m_n1|, …, |m_1n| + … + |m_nn|}，
而向量∞-范数所诱导出的n阶方阵M = [m_ij]的∞-范数为
||M||_∞ = max {|m₁₁| + … + |m_1n|, …, |m_n1| + … + |m_nn|}。
所以，只要迭代矩阵M满足||M||₁ < 1或||M||_∞ < 1，则迭代法对任意的初始向量都收敛。

特征值与谱半径然而，有一个收敛的充分必要条件等待着我们去挖掘，这个条件不再用到矩阵的范数，而是矩阵的谱。谱给出了矩阵不依赖于矩阵范数的一个内蕴性质。方阵的谱是它所有特征值全体组成的一个有限的复数集合。给定n阶方阵M，我们称复数λ为M的一个特征值，如果复矩阵M - λI不是可逆的，即存在一个复的非零列向量x使得Mx = λx。这个x被称为M对应于特征值λ的一个特征向量。我们曾经提到，一个矩阵是奇异的当且仅当它的行列式等于零，故λ是M的特征值当且仅当det(M - λI) = 0，其中符号det表示行列式。如果把这个等式左边中的λ看成是变元，根据行列式的定义，det(M - λI)的展开结果是关于λ的一个n阶多项式，所以一个n阶方阵M顶多有n个相异的特征值。我们把M的所有特征值绝对值中的最大值称作 M的谱半径，记为ρ(M)。
研究线性迭代的主要目的是数值求解线性方程组Ax = b，其中的系数矩阵A为非奇异的，这样保证对所有的右端常向量b，该方程组有并且仅有一个解，它就是p = A^-1b。为了设计一个迭代法，首先将A分裂成N - P的形式，其中N也是非奇异的。然后原方程组等价于不动点线性方程组x = Mx + c，其中M = N^-1P和c = N^-1b。基本的A非奇异性假设推出I - M也是非奇异的，故由F(x) = Mx + c定义的仿射向量函数F: Rⁿ → Rⁿ有并且仅有一个不动点向量p。现在我们可以证明：线性迭代
x_k = Mx_k-1 + c, k = 1, 2, 3, …
对所有初始列向量x₀都收敛的充分必要条件是ρ(M) < 1。
为了证明这个断言，先设迭代法对所有的初始列向量x₀都收敛，则极限列向量满足不动点方程x = Mx + c。又因为如上所述I - M是非奇异矩阵，这个极限向量必定是F唯一的不动点p。令e_k = x_k - p为迭代到第k步时的误差向量，那么对所有的列向量e₀，极限 e_k = Me_k-1 = M^ke₀ → 0都成立。如果ρ(M) ≥ 1，则存在一个特征值λ满足|λ| ≥ 1，也存在它所对应的一个特征向量e₀。因为e₀是非零向量，然后我们就会发现，e_k = M^ke₀ = λ^ke₀当k趋向于无穷大时不可能趋向于0，这与假设矛盾。所以ρ(M) < 1，必要性得证。
反过来，如果ρ(M) < 1，则存在一个足够小的正数ε 使得ρ(M) + ε < 1。这一看似简单的事实是实数的基本性质。另一方面，对于M的每一个特征值λ和Rⁿ上的任一个范数||•||，设x为λ所对应的一个特征向量，由不等式||M|| ||x|| ≥ ||Mx|| = ||λx|| = |λ| ||x||以及x ≠ 0可知，刚刚获得的不等式两边可以除以正数||x||，而得到不等式||M|| ≥ |λ|，从而可知：矩阵M的谱半径ρ(M)是所有向量范数诱导出的矩阵范数在M的值||M||的一个下界。事实上，这个下界是对应于固定矩阵M的全部范数值||M||构成的实数集合的“下确界”，即这个集合的所有下界中的“最大下界”。“上、下确界”的概念是微积分学的最基础、最重要的概念，真正学通了它，极限理论以及随之而来的连续性、导数、积分、级数等等概念学起来甚至可以“无师自通”了。在稍微高等一点的矩阵理论中，可以利用M的“若尔当标准型”，构造Rⁿ上的一个范数||•||，使得它诱导出的矩阵范数满足条件||M|| ≤  ρ(M)  + ε < 1（这个对任意正数ε的矩阵范数构造，恰恰证明了上述的“谱半径是所有矩阵范数值的下确界”断言）。故对任何初始列向量e₀有
||e_k|| = ||M^ke₀|| ≤ ||M||^k ||e₀|| ≤ (ρ(M) + ε)^k ||e₀|| → 0。
这就证明了充分性。
有一个充分体现矩阵谱半径与矩阵范数亲密关系的等式：。谱半径是矩阵固有的内蕴性质，而范数却是“强加于”矩阵身上的一个外在标尺，通过极限之桥梁，本质透过外观而表达。更神的是，这个关系式与维数无关，对无穷维巴拿赫空间上的有界线性算子同样正确。
另外可以证明，迭代法对所有初始向量都收敛的另一个充分必要条件是：矩阵幂次序列{M^k}当k趋向于无穷大时收敛到零矩阵，但这个等价条件远不及迭代矩阵谱半径小于1的条件简单实用。
我们想要提醒读者注意的是，上面的收敛性谱半径等价条件为真的一个前提是逆矩阵(I - M)^-1存在，否则的话，反例多得很，比如令M为一的对角矩阵，其对角元素为1和1/2，则它的谱半径不小于1。易知I - M是一个对角线元素为0和1/2的奇异矩阵。另一方面，如果让c = (0, 1)^T，则对任意的初始向量x₀ = (a, b)^T，迭代x_k = Mx_k-1 + c, k = 1, 2, …产生的向量序列x_k = (a, b/2^k + 2 - 1/2^k-1)^T → (a, 2)^T，说明这个迭代对所有的初始向量都收敛，但极限向量不唯一， 它们构成xy-平面内的一维仿射集{(a, 2)^T: a是任意实数}。

雅可比迭代法和高斯-赛德尔迭代法现在我们可以集中讨论求解线性方程组Ax = b的雅可比迭代法和高斯-赛德尔迭代法了，其中A为n阶可逆方阵。这两种产自德意志的经典迭代法都是基于在矩阵分解A = N - P中对N和P的特殊选取，其迭代格式都可写为
x_k = N^-1Px_k-1 + N^-1b, k = 1, 2, 3, …。
根据如上的谱分析结论，迭代法对所有初始点收敛当且仅当ρ(N^-1P) < 1。
用D表示由A中的对角元素a₁₁, a₂₂, …, a_nn保持原先位置和次序而得到的相应对角矩阵，L表示A中严格位于对角线之下的所有元素a_ij (i > j)按照原先位置而构成的严格下三角矩阵，U表示A中严格位于对角线之上的所有元素a_ij (i < j)按照原先位置而构成的严格上三角矩阵。这里习惯采用的三个字母分别是三个英文单词diagonal（对角）、lower（下）和 upper（上）的大写首字母，一目了然。于是，A = D + L + U = N - P。请注意L和U的所有对角元素为0，此结构的一个推论是Lⁿ = Uⁿ = 0，这将在后面用到。
为了实现上述两种迭代法，必须先假设A的对角元都不为零。对于雅可比方法，N = D和 P = -(L + U)，故迭代矩阵N^-1P = -D^-1(L + U) = I - D^-1A；而对于高斯-赛德尔方法，N = D + L和P = -U，故迭代矩阵N^-1P = -(D + L)^-1U。雅可比迭代法的计算格式是
x_k = (I - D^-1A)x_k-1 + D^-1b, k = 1, 2, 3, …；
高斯-赛德尔迭代法的计算格式是
x_k = -(D + L)^-1Ux_k-1 + (D + L)^-1b, k = 1, 2, 3, …。
我们试图写出上述两个迭代法的分量迭代公式。对i, j = 1, 2, …, n, 记A的第i行、第j列元素为a_ij。为了避免矩阵求逆运算，将雅可比迭代的循环公式写成如下形式
Dx_k = -(L + U)x_k-1 + b, k = 1, 2, 3, …。
如果将第k步迭代向量x_k的n个分量写成x_k,1, x_k,2, …, x_k,n，则对i = 1, 2, …, n，上述等式的第i个分量等式是

两边除以a_ii，我们得到雅可比方法的各分量迭代公式
 至于高斯-赛德尔方法，我们先将其迭代格式改写成等价的“无逆矩阵形式”
(D + L)x_k = -Ux_k-1 + b, k = 1, 2, 3, …
或另一种写法
Dx_k = b - Lx_k - Ux_k-1, k = 1, 2, 3, …。
故对所有的i = 1, 2, …, n，有

这样，高斯-赛德尔方法的各分量迭代公式是

我们比较一下雅可比方法与高斯-赛德尔方法在迭代点分量计算过程中的显著差别。由雅可比迭代法中n个分量的迭代公式可见，在第k - 1个迭代向量所有的分量都计算完毕后，第k个迭代向量的各个分量才一个接着一个地计算出来；换句话说，第k - 1步迭代获得的全部分量都被用于对第k步迭代所有分量的计算。而对于高斯-赛德尔迭代法，每当需要算出第k个迭代向量的第i个分量时，不仅需要已完成计算的第k - 1个迭代向量从第i + 1到第n个分量帮忙，而且还需要正在进行中的本次迭代已得到的第k个迭代向量的第1到第i - 1个分量相助；或言之，高斯-赛德尔法比雅可比法“更急于求成”，命令当前迭代步骤中刚刚收入囊中的迭代点分量取代上一迭代步骤被“擒获”的同下标分量，“披挂上阵”。这说明，高斯和赛德尔比雅可比更会“用兵”，不肯浪费半点士气。

收敛条件自然，无论是雅可比迭代法还是高斯-赛德尔迭代法，对矩阵A缺乏某些额外的假设条件，都不能保证其收敛性。那么，当“A的对角线元素个个不能为零”这个先决条件被满足后，什么是保证它们收敛的充分条件呢？
我们已经知道，它们收敛的充分必要条件分别是ρ(I - D^-1A) < 1和ρ[-(D + L)^-1U] < 1。然而，矩阵的谱半径不是那么好得到的，因为它用到所有特征值。但是矩阵的1-范数和∞-范数的计算可以说是举手之劳。由于矩阵I - D^-1A的所有对角元素均是0，我们马上发现它的1-范数等于，∞-范数等于。因此我们获得雅可比方法收敛的一个充分条件：
如果下列条件

之一满足，则雅可比迭代法对所有的初始列向量都收敛到Ax = b的唯一解p = A^-1b。
上述条件(ii)可以用来定义一类矩阵。一个n阶方阵A被称为严格对角占优的条件是，对所有的i = 1, 2, …, n，严格不等式都成立。由此可见，如果线性方程组Ax = b中的矩阵A是严格对角占优的，那么雅可比迭代法收敛。这时，A的非奇异性自动由严格对角占优的假设保证，因为A = D[I - (I - D^-1A)]是两个非奇异矩阵的乘积，后一因子矩阵I - (I - D^-1A)的非奇异性是由于，从不等式ρ(I - D^-1A) ≤ ||I - D^-1A||_∞ < 1可以推出：1非I - D^-1A的特征值，故而I - (I - D^-1A)没有特征值0。此外还可证明，若严格对角占优矩阵A是实对称（即A^T = A）的并且对角元素均为正数，那么A必定是正定矩阵，即对所有的n维非零列向量x，严格不等式x^TAx > 0都成立。
这里，保证雅可比迭代法收敛性的严格对角占优的要求可以被改成不可约弱对角占优，即上面的所有严格不等式换为一般不等式 “≥”，且至少有一个不等式严格成立，但增加了矩阵不可约的附加条件。如果一个方阵的行和列能以同样的方式重新排列（不再排列也算“重新排列”，因为单位矩阵也是排列矩阵），结果变成一个2×2块上三角矩阵，即矩阵的左下角有至少为1×1的0矩阵，那么我们说这个矩阵是可约的，否则就说它是不可约的。
至于高斯-赛德尔迭代法，我们下面花点笔墨，证明在A为严格对角占优的假设下，它的迭代矩阵-(D + L)^-1U的∞-范数不大于雅可比迭代矩阵I - D^-1A的∞-范数。这个证明主要基于“序 (order) ”的概念。法国布尔巴基数学学派认为纯粹数学主要研究三种结构：关于代数运算的代数结构、关于连续性质的拓扑结构、关于大小关系的有序结构。两个同一尺寸的实矩阵B和C若被写成B ≤ C，指的是B的每一个元素b_ij小于或等于C的对应元素c_ij。对于B，我们记|B|为它的绝对值矩阵，即|B|的每一个元素为B的对应元素b_ij的绝对值|b_ij|。注意，|B|是非负矩阵，而||B||是非负数。另外，我们用字母e表示所有分量均为1的任意维数的列向量。
我们将雅可比方法和高斯-赛德尔方法的迭代矩阵分别记为M₁ = I - D^-1A = -D^-1(L + U) = -D^-1L - D^-1U和M₂ = -(D + L)^-1U = -(I + D^-1L)^-1D^-1U。因为|M₁|的每行所有元素之和都不大于它们中的最大值||M₁||_∞，而一个矩阵乘以特殊列向量e的结果就是将该矩阵每行加起来的和所构成的列向量，故有向量不等式
|M₁| e ≤ ||M₁||_∞ e。
另一方面，由于-D^-1L和-D^-1U分别为严格下三角矩阵和严格上三角矩阵，显然有等式
|M₁| = |-D^-1L - D^-1U| = |D^-1L| + |D^-1U|。
上述的不等式和等式推出
|D^-1U|e = (|M₁| - |D^-1L|)e ≤ [I - |D^-1L| - (1 - ||M₁||_∞)I]e。   (1)
由于-D^-1L和|D^-1L|都是严格下三角矩阵，(-D^-1L)ⁿ 和|D^-1L|ⁿ均为零矩阵，故有
(I + D^-1L)[I - D^-1L + (D^-1L)² - … + (-D^-1L)^n-1] = I，(I - |D^-1L|)[I + |D^-1L| + |D^-1L|² + … + |D^-1L|^n-1] = I。
因此矩阵I + D^-1L和I - |D^-1L|均为可求逆的，且有
|(I + D^-1L)^-1| = |I - D^-1L + (D^-1L)² - … + (-D^-1L)^n-1|≤ I + |D^-1L| + |D^-1L|² + … + |D^-1L|^n-1= (I - |D^-1L|)^-1,
从而上式以及不等式(1)保证
|M₂|e = |-(I + D^-1L)^-1D^-1U|e ≤ |(I + D^-1L)^-1|  |D^-1U|e≤ (I - |D^-1L|)^-1[(I - |D^-1L|) - (1 - ||M₁||_∞)I]e= [I - (1 - ||M₁||_∞)(I - |D^-1L|)^-1]e。
因为由假设||M₁||_∞ < 1以及(I - |D^-1L|)^-1 = I + |D^-1L| + |D^-1L|² + … + |D^-1L|^n-1 ≥ I，我们得到
|M₂|e ≤ ||M₁||_∞e，
从而||M₂||_∞ = max{(|M₂|e)_i: i = 1, 2, …, n} ≤ ||M₁||_∞，故我们证明了下列高斯-赛德尔迭代法收敛的一个简单充分条件：
如果雅可比迭代法的收敛条件(ii)成立，则高斯-赛德尔迭代法对任意的初始列向量都收敛到Ax = b的唯一解p = A^-1b。
当A为严格对角占优矩阵时，高斯-赛德尔迭代矩阵的∞-范数不大于雅可比迭代矩阵的∞-范数，而这个范数越小，收敛速度通常就越快，因而就收敛快慢而言，前一种方法不亚于后一种方法。
我们把高斯-赛德尔方法另一个独特的收敛性命题端上餐桌，作为本文的最后一道数学菜肴：
如果线性方程组Ax = b的系数矩阵A是对称正定的，那么高斯-赛德尔迭代法收敛。
在证明它前，我们指出，对实对称正定矩阵A和所有复向量x ≠ 0，也有x^HAx > 0，其中x^H是x的共轭转置。在此强调之，是因为后面的特征向量一般是复的。在下列的推理中，我们常用到简单事实：对两个可乘的复矩阵B和C，有(BC)^H = C^HB^H；对任意复向量x，二次型x^HBx是个复数；对实矩阵B，有B^H = B^T；对一阶复矩阵[b]，有。
因为A 为对称正定矩阵，U = L^T。令-λ为高斯-赛德尔迭代矩阵M₂ = -(L + D)^-1L^T的一个复特征值，且x为对应的一个复特征向量，则L^Tx = λ(L + D)x，因而
x^HL^Tx = λx^H(L + D)x。  (2)
上式两端各加上x^H(L + D)x，并由于A = L + D + L^T，得0 < x^HAx = (1 + λ)x^H(L + D)x。故λ ≠ -1且等式右端等于其共轭转置，即

由上式及可得

故合并左右两端的同类型给出
(1 - |λ|²) x^H(L + D)^Tx = (1 + λ)x^HDx。
所以，利用等式(2)和上式，有
(1 - |λ|²) x^HAx = (1 - |λ|²) [x^H(L + D)x + x^HL^Tx]= (1 +)x^HDx + (1 - |λ|²) λx^H(L + D)x= (1 + )x^HDx + λ(1 + ) x^HDx= (1 + + λ + |λ|²) x^HDx = |1 + λ|² x^HDx。 因为上述等式的右端是正数，所以左端的1 - |λ|² > 0，即|-λ| = |λ| < 1。由于-λ是M₂的任一特征值，这就证明了ρ(M₂) < 1。迭代矩阵谱半径小于1等价于线性迭代法收敛，命题得证。

尾声德国数学家在两百年前点燃的线性迭代法火炬，百年之后远渡重洋，传递到了科技腾飞、电脑出世的美洲大陆。为了优化迭代效率，从经典的高斯-赛德尔方法出发，在1950年，美国数学家、计算机科学家杨大卫 (David M. Young Jr.，1923-2008，“返朴”将另行介绍他) 和美国物理学家、计算机科学家弗兰克尔 (Stanley Phillips Frankel，1919-1978) 几乎同时引入了一个松弛因子ω进行某种仿射组合，引出了“逐次超松弛迭代法 (method of successive over-relaxation，SOR) ”。前者在哈佛大学数学系的博士学位论文Iterative Methods for Solving Partial Difference Equations of Elliptic Type（《求解椭圆型偏差分方程的迭代方法》）中提出了SOR方法；后者在美国数学会期刊（十年后改名为Mathematics of Computation[《计算数学》]）上发表论文Convergence rates of iterative treatments of partial differential equations（《偏微分方程迭代处理的收敛率》），其中对他设计并称之为“Extrapolated Liebmann method（外推利伯曼法）的SOR方法进行了全面的分析。这两项提出SOR方法的先驱性研究，都与在科学和工程中大量出现的偏微分方程有关，用电子计算机求解这些连续方程离散化后的大型稀疏线性方程组，迭代法是首要之选，这就是SOR方法应运而生的时代背景。
当松弛因子ω取值为1时，SOR方法就还原成高斯-赛德尔方法。我们不细表此法，就此结束本文，因为我们已经通过了解两个古典的迭代法而知悉线性迭代法收敛性准则的基本思想。总而言之，十九世纪三名德国数学家的发明创造，是现代大型、稀疏、结构化线性方程组迭代解研究浪潮的源头。
写于2023年11月19日星期日美国哈蒂斯堡夏日山庄

本文受科普中国·星空计划项目扶持

出品：中国科协科普部

监制：中国科学技术出版社有限公司、北京中科星河文化传媒有限公司

相关阅读

近期推荐

特 别 提 示

1. 进入『返朴』微信公众号底部菜单“精品专栏“，可查阅不同主题系列科普文章。

2. 『返朴』提供按月检索文章功能。关注公众号，回复四位数组成的年份+月份，如“1903”，可获取2019年3月的文章索引，以此类推。

版权说明：欢迎个人转发，任何形式的媒体或机构未经授权，不得转载和摘编。转载授权请在「返朴」微信公众号内联系后台。

找不到《返朴》了？快加星标！！

长按下方图片关注「返朴」，查看更多历史文章

微信实行乱序推送，常点“在看”，可防失联