8.4 帕德逼近

简介

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

例子

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

通解

  首先我们要先求麦克劳林级数，这个场景太多了，我就不详细说了，假设求得得麦克劳林级数为:
$$f(x)=c_0+c_{1}x+c_x{x}^2+\cdots +c_{m+n}{x}^{m+n}+O(x^{m+n+1})$$
  然后解方程组求得分母的所有系数，使用LUP分解可以快速求解：
$$\left[\begin{array}{ccc}{c}_n& \dots & c_1\\ ⋮& \ddots & ⋮\\ c_{2n-1}& \cdots & c_n\end{array}\right]\left[\begin{array}{c}{b}_1\\ ⋮\\ b_n\end{array}\right]=\left[\begin{array}{c}-c_{n+1}\\ ⋮\\ -c_{2n}\end{array}\right]$$
  而$b_0=1$，如此，求得了分母了，然后按下式带入，求得分子：
$$\begin{gathered} a_0=c_0 \\ {a}_1=c_0{b}_1+c_1 \\ ⋮ \\ {a}_m=\sum _{i=0}^m{c}_i{b}_{m-i} \end{gathered}$$

python代码

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

# -*- coding:UTF-8 -*-
import copy

import sympy
from sympy.core.expr import Expr

MAX_ERROR = 1e-15


def solve(l, u, values):
    n = len(values)
    y = [0] * n
    for i in range(0, n):
        #
        # 遍历其后的所有行,让values的值减掉
        y[i] = values[i]
        for j in range(i + 1, n):
            values[j] -= l.lines[j][i] * values[i]
    # 再解Ux=y
    solution = [0] * n
    for i in range(n - 1, -1, -1):
        solution[i] = round(y[i] / u.lines[i][i], 2)
        for j in range(0, i):
            y[j] -= u.lines[j][i] * solution[i]
    return solution


class Matrix:
    # 矩阵

    def __init__(self, lines):
        self.__lines = lines

    def __str__(self):
        result: str = ''
        for line in self.__lines:
            result += str(line)
            result += '\n'
        return result

    # 暴露属性
    @property
    def lines(self):
        return self.__lines

    def append(self, matrix):
        # define an array
        rows = len(self.__lines)
        columns = len(self.__lines[0]) + len(matrix.__lines[0])
        unit = [[0 for _ in range(0, columns)] for _ in range(0, rows)]
        for y in range(0, len(unit)):
            self_line = self.__lines[y]
            other_line = matrix.__lines[y]
            # 0~ this this ~other
            for i in range(0, len(self_line)):
                unit[y][i] = self_line[i]

            for i in range(0, len(other_line)):
                unit[y][i + len(self_line)] = other_line[i]
        return Matrix(unit)

    def sub_matrix(self, row_start, row_end, column_start, column_end):
        unit = [[self.__lines[i][j] for j in range(column_start, column_end)] for i in range(row_start, row_end)]
        return Matrix(unit)

    def __sub__(self, other):
        return Matrix([[e - other.__lines[y][x] for x, e in enumerate(line)] for y, line in enumerate(self.__lines)])

    # LUP分解
    def lup_decomposition(self):

        n = len(self.__lines)
        # 初始化L U为同样大小的0矩阵
        l = copy.deepcopy(self.__lines)
        u = copy.deepcopy(self.__lines)
        a = copy.deepcopy(self.__lines)
        p = [i for i in range(0, n)]
        # i代表处理的行和列，因为每次都是从左上到右下不断缩小这个矩阵
        for i in range(0, n):
            # 先求出自己的置换
            # 注意不能用是否等于0判断，而是考虑浮点数精度进行判断
            # 找最大绝对值
            swap_line = i

            # 进行置换
            # 找出不为0的一行，进行置换
            for k in range(i + 1, n):
                if abs(a[k][i]) > abs(a[swap_line][i]):
                    swap_line = k
            if abs(a[swap_line][i]) <= MAX_ERROR:
                raise Exception("矩阵无解")
            # 交换两行
            a[i], a[swap_line] = a[swap_line], a[i]
            l[i], l[swap_line] = l[swap_line], l[i]
            p[i], p[swap_line] = p[swap_line], p[i]

            # 剩下的就是LU分解了
            # 每一行的l进行赋值
            # 这个循环其实是对角线方向
            l[i][i] = 1
            u[i][i] = a[i][i]
            for j in range(i + 1, n):
                l[i][j] = 0
                l[j][i] = a[j][i] / a[i][i]
                u[i][j] = a[i][j]
                u[j][i] = 0
                # 处理A比较麻烦，需要计算舒尔补
                for k in range(i + 1, n):
                    # 舒尔补是 a[j][k] -= a[j][i] * a[i][k] / a[i][i]
                    # 因为前面计算l[j][i] = a[j][i] / a[i][i]
                    # 所以可以直接替换
                    a[j][k] = round(a[j][k] - l[j][i] * a[i][k], 2)
        return Matrix(l), Matrix(u), Matrix([[1 if j == i else 0 for j in range(0, n)] for i in p])

    def solve_by_lup(self, values):
        l, u, p = self.lup_decomposition()
        # lux=pb
        pb = p * Matrix([[num] for num in values])
        pb_values = [i[0] for i in pb.__lines]
        # 因为l是单位下三角矩阵，所以直接代入法
        return solve(l, u, pb_values)

  接下来是我的帕德代码，里面包含了有理函数类：

from com.youngthing.mathalgorithm.matrix import Matrix


class Rational:

    def __init__(self, a, b):
        self.__numerator = list(reversed(a))
        self.__denominator = list(reversed(b))

    def __call__(self, *args, **kwargs):
        return Rational.calc(self.__numerator, args[0]) / Rational.calc(self.__denominator, args[0])

    def __str__(self):
        a_str = Rational.poly_str(self.__numerator)
        b_str = Rational.poly_str(self.__denominator)
        return f"({a_str})/({b_str})"

    @staticmethod
    def calc(polynomial, val):
        result = 0
        for coefficient in polynomial:
            result = result * val + coefficient
        return result

    @staticmethod
    def poly_str(polynomial):
        result = ''
        first = True
        n = len(polynomial)
        for i, coefficient in enumerate(polynomial):
            if coefficient == 0:
                continue
            if first:
                first = False
            else:
                result += '+'
            if coefficient != 1:
                result += str(coefficient) + '*'
            order = n - i - 1
            if order > 0:
                result += 'x'
            else:
                result += '1'
            if order > 1:
                result += '**' + str(order)
        return result



def pade(coefficients, m, n):
    # 第一步解线性方程组，得出分母所有系数
    lines = []
    for i in range(0, n):
        line = list(reversed(coefficients[i + 1:i + n + 1]))
        lines.append(line)
    matrix = Matrix(lines)
    values = [-x for x in coefficients[m + 1:m + n + 1]]
    b = matrix.solve_by_lup(values)
    b.insert(0, 1)
    a = []
    for i in range(0, m + 1):
        a_i = 0
        for j in range(0, i + 1):
            a_i += coefficients[j] * b[i - j]
        a.append(a_i)
    return Rational(a, b)

测试

  我们以sin(x)的$3,3$阶为例子，进行测试。首先计算$\sin (x)$的麦克劳林级数：
$$sin(x)=x-\frac{x^3}{3!}+\frac{x^5}{5!}+O(x^7)$$
  那么测试数据输入单元测试，写出测试代码：

import unittest

from com.youngthing.mathalgorithm.approximate.pade import pade
from com.youngthing.mathalgorithm.permutations.heap import permutations


class MyTestCase(unittest.TestCase):
    def test_sin_3_3(self):
        coefficients = [0, 1, 0, -1 / 6, 0, 1 / 120, 0]
        r = pade(coefficients, 3, 3)
        print(r)

    if __name__ == '__main__':
        unittest.main()

  测试结果：

(-0.10666666666666666*x**3+x)/(0.06*x**2+1)

  用LaTex plot一下，看看结果：

\documentclass[a4paper,UTF8]{article}
\usepackage{ctex}
\usepackage{tikz}
\begin{document}
	\begin{tikzpicture}
		
		% 画坐标系
		\draw (0,-3) [->, thick]-- (0,3) node (yaxis) [above] {y};
		\draw (-4,0) [->, thick]-- (4,0) node (yaxis) [above] {x};
	
		% 画X轴
		\foreach \x in {-4,-3,-2,-1,1,2,3,4}{
			\node at (\x,0) [below] {\x};
			\draw (\x,0) -- (\x, 2pt);
		}
			
		\foreach \y in {-3,-2,-1,0,1,2,3}{
			\node at (0,\y) [left] {\y};
		}
		\foreach \y in {-3,-2,-1,1,2,3}{
			\draw (0,\y) -- (2pt,\y);
		}
		
		\draw [blue,smooth,domain=-4:4] plot function {sin(x)};
		\draw [green,smooth,domain=-4:4] plot function {(-0.10666666666666666*x**3+x)/(0.06*x**2+1)};
	
	
	\end{tikzpicture}
\end{document}

  结果如下，可以看到，比$2/2$阶的帕德逼近精确多了：