8.4 帕德逼近

  帕德逼近Padé approximant是一种对任意函数的有理函数逼近。这个是高中数学内容，经常在高考题中出现，但是呢，很多人忘掉了，以为是高等数学内容。有理函数不消我多说了吧，是分子分母都是多项式的函数。帕德逼近有阶的概念，如果分子是m阶多项式，分母是n阶多项式，那么帕德逼近就是 m / n m/n m/n阶的帕德逼近，符号为 [ m / n ] [m/n] [m/n]或者 R m , n R_{m,n} Rm,n​。当n=0的时候，就是多项似乎逼近了。多项式逼近有多种，常见的有切比雪夫逼近和泰勒级数，帕德逼近在n=0时，是泰勒级数的前m项。   帕德逼近还有个小细节，它的基本形式如下： [ m / n ] = R m , n = p 0 + p 1 x + p 2 x 2 + ⋯ + p m x m 1 + q 1 x + q 2 x 2 + ⋯ + q n x n [m/n]=R_{m,n}=\frac{p_0+p_1x+p_2x^2+\dots+p_mx^m}{1+q_1x+q_2x^2+\dots+q_nx^n} [m/n]=Rm,n​=1+q1​x+q2​x2+⋯+qn​xnp0​+p1​x+p2​x2+⋯+pm​xm​   这个细节就是分母是以1开始的，也就是 q 0 = 1 q_0=1 q0​=1。   与切比雪夫近似值不同的是，帕德逼近没有限定定义域，也就是在目标函数的整个定义域内有唯一的帕德逼近。而切比雪夫近似值是目标函数在不同的定义域内有不同的切比雪夫近似值。帕德逼近是以0为基础的，在 x = 0 x=0 x=0点的误差最小，离0越远，误差越大。   讲了这么多，帕德逼近怎么求呢？切比雪夫近似值的求法是很麻烦的，要用到复杂的雷米兹算法，但是帕德逼近的求法很简单。我拿个例子来说吧。

  以求 f ( x ) = s i n ( x ) f(x)=sin(x) f(x)=sin(x)的 [ 2 / 2 ] [2/2] [2/2]阶帕德逼近为例子，我讲讲帕德逼近的过程。因为分母有2阶，而分子也有2阶，所以总体复杂度有4阶。而 s i n ( x ) sin(x) sin(x)是一个很“无理”的函数，所以要让他“有理”。那么有两种办法，第一种办法是泰勒级数，第二种办法是切比雪夫近似值。因为切比雪夫近似值需要有定义域，所以直接排除，那么只剩下泰勒级数了。那么我们就展开到4阶吧。 s i n ( x ) = x − x 3 3 ! + O ( x 5 ) sin(x) = x-\frac{x^3}{3!}+O(x^5) sin(x)=x−3!x3​+O(x5)   它的 3 / 3 3/3 3/3阶帕德逼近形式为： R 2 , 2 ( x ) = a ( x ) b ( x ) = a 0 + a 1 x + a 2 x 2 1 + b 1 x + b 2 x 2 R_{2,2}(x)=\frac{a(x)}{b(x)}\\ =\frac{a_0+a_1x+a_2x^2}{1+b_1x+b_2x^2} R2,2​(x)=b(x)a(x)​=1+b1​x+b2​x2a0​+a1​x+a2​x2​   要求泰勒级数去掉尾项后减去帕德逼近要约等于0，于是误差只在泰勒级数尾项中： x − x 3 3 ! − a 0 + a 1 x + a 2 x 2 1 + b 1 x + b 2 x 2 ≈ 0 x-\frac{x^3}{3!}-\frac{a_0+a_1x+a_2x^2}{1+b_1x+b_2x^2}\approx0 x−3!x3​−1+b1​x+b2​x2a0​+a1​x+a2​x2​≈0   移动一下，将约等号变成等号： x − x 3 3 ! = a 0 + a 1 x + a 2 x 2 1 + b 1 x + b 2 x 2 ( x − x 3 3 ! ) ( 1 + b 1 x + b 2 x 2 ) = a 0 + a 1 x + a 2 x 2 x-\frac{x^3}{3!}=\frac{a_0+a_1x+a_2x^2}{1+b_1x+b_2x^2}\\ (x-\frac{x^3}{3!})({1+b_1x+b_2x^2})={a_0+a_1x+a_2x^2} x−3!x3​=1+b1​x+b2​x2a0​+a1​x+a2​x2​(x−3!x3​)(1+b1​x+b2​x2)=a0​+a1​x+a2​x2   将左边展开： ( x − x 3 3 ! ) ( 1 + b 1 x + b 2 x 2 ) = x − 1 6 x 3 + b 1 x 2 − 1 6 b 1 x 4 + b 2 x 3 − 1 6 b 2 x 5 + = x + b 1 x 2 + ( b 2 − 1 6 ) x 3 − 1 6 b 1 x 4 − 1 6 b 2 x 5 (x-\frac{x^3}{3!})({1+b_1x+b_2x^2})\\ =x-\frac16x^3+\\ b_1x^2-\frac16b_1x^4+\\ b_2x^3-\frac16b_2x^5+\\ =x+b_1x^2+(b_2-\frac16)x^3-\frac16b_1x^4-\frac16b_2x^5 (x−3!x3​)(1+b1​x+b2​x2)=x−61​x3+b1​x2−61​b1​x4+b2​x3−61​b2​x5+=x+b1​x2+(b2​−61​)x3−61​b1​x4−61​b2​x5   所以有： x + b 1 x 2 + ( b 2 − 1 6 ) x 3 − 1 6 b 1 x 4 − 1 6 b 2 x 5 = a 0 + a 1 x + a 2 x 2 x+b_1x^2+(b_2-\frac16)x^3-\frac16b_1x^4-\frac16b_2x^5\\ = a_0+a_1x+a_2x^2 x+b1​x2+(b2​−61​)x3−61​b1​x4−61​b2​x5=a0​+a1​x+a2​x2   左边减去右边，有： − a 0 + ( 1 − a 1 ) x + ( b 1 − a 2 ) x 2 + ( b 2 − 1 6 ) x 3 − 1 6 b 1 x 4 − 1 6 b 2 x 5 = 0 -a_0+(1-a_1)x+(b_1-a_2)x^2+(b_2-\frac16)x^3-\frac16b_1x^4-\frac16b_2x^5=0 −a0​+(1−a1​)x+(b1​−a2​)x2+(b2​−61​)x3−61​b1​x4−61​b2​x5=0   把最高次的 x 5 x^5 x5项忽略，剩余的每项等于0，可以得到一个方程组： − a 0 = 0 1 − a 1 = 0 b 1 − a 2 = 0 b 2 − 1 6 = 0 − 1 6 b 1 = 0 -a_0=0\\ 1-a_1=0\\ b_1-a_2=0\\ b_2-\frac16=0\\ -\frac16b_1=0\\ −a0​=01−a1​=0b1​−a2​=0b2​−61​=0−61​b1​=0   把上面的方程组解开得： b 1 = 0 b 2 = 1 6 b_1=0\\ b_2=\frac1{6} b1​=0b2​=61​   代入得： a 0 = 0 a 1 = 1 a 2 = 0 a_0=0\\ a_1=1\\ a_2=0 a0​=0a1​=1a2​=0   于是最终 s i n ( x ) sin(x) sin(x)的 2 / 2 2/2 2/2阶帕德逼近为： R 2 , 2 ( x ) = x 1 + 1 6 x 2 R_{2,2}(x)=\frac{x}{1+\frac{1}{6}x^2} R2,2​(x)=1+61​x2x​   在坐标系上，帕德逼近和原函数在越接近于0的地方，误差越小，离0越远，就误差越大，越来越离，如 R 2 , 2 R_{2,2} R2,2​与 s i n ( x ) sin(x) sin(x)的图像如下：

  首先我们要先求麦克劳林级数，这个场景太多了，我就不详细说了，假设求得得麦克劳林级数为: f ( x ) = c 0 + c 1 x + c x x 2 + ⋯ + c m + n x m + n + O ( x m + n + 1 ) f(x)=c_0+c_1x+c_xx^2+\dots+c_{m+n}x^{m+n}+O(x^{m+n+1}) f(x)=c0​+c1​x+cx​x2+⋯+cm+n​xm+n+O(xm+n+1)   然后解方程组求得分母的所有系数，使用LUP分解可以快速求解： [ c n … c 1 ⋮ ⋱ ⋮ c 2 n − 1 ⋯ c n ] [ b 1 ⋮ b n ] = [ − c n + 1 ⋮ − c 2 n ] \begin{bmatrix} c_n & \dots & c_1\\ \vdots & \ddots & \vdots\\ c_{2n-1}&\cdots&c_n\\ \end{bmatrix}\begin{bmatrix} b_1\\ \vdots\\ b_n \end{bmatrix}= \begin{bmatrix} -c_{n+1}\\ \vdots\\ -c_{2n} \end{bmatrix} ⎣ ⎡​cn​⋮c2n−1​​…⋱⋯​c1​⋮cn​​⎦ ⎤​⎣ ⎡​b1​⋮bn​​⎦ ⎤​=⎣ ⎡​−cn+1​⋮−c2n​​⎦ ⎤​   而 b 0 = 1 b_0=1 b0​=1，如此，求得了分母了，然后按下式带入，求得分子： a 0 = c 0 a 1 = c 0 b 1 + c 1 ⋮ a m = ∑ i = 0 m c i b m − i a_0=c_0\\ a_1=c_0b_1+c_1\\ \vdots\\ a_m=\sum_{i=0}^{m}c_ib_{m-i } a0​=c0​a1​=c0​b1​+c1​⋮am​=i=0∑m​ci​bm−i​

  python代码就是根据泰勒级数求帕德逼近的 m + n + 1 m+n+1 m+n+1个未知参数。当然，根据方程求麦克劳林级数（在 x = 0 x=0 x=0点的泰勒级数）的代码我是不会写的，因为场景太多太复杂。对于解方程，我这里复制了一个简单的矩阵类：