文章目录
简介
帕德逼近 Padé approximant_是一种对任意函数的有理函数逼近。这个是高中数学内容，经常在高考题中出现，但是呢，很多人忘掉了，以为是高等数学内容。有理函数不消我多说了吧，是分子分母都是多项式的函数。帕德逼近有阶的概念，如果分子是m阶多项式，分母是n阶多项式，那么帕德逼近就是m / n m/n m /n阶的帕德逼近，符号为[ m / n ] [m/n][m /n ]或者R m , n R{m,n}R m ,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 ]=R m ,n =1 +q 1 x +q 2 x 2 +⋯+q n x n p 0 +p 1 x +p 2 x 2 +⋯+p m x m
这个细节就是分母是以1开始的，也就是q 0 = 1 q_0=1 q 0 =1。
与切比雪夫近似值不同的是，帕德逼近没有限定定义域，也就是在目标函数的整个定义域内有唯一的帕德逼近。而切比雪夫近似值是目标函数在不同的定义域内有不同的切比雪夫近似值。帕德逼近是以0为基础的，在x = 0 x=0 x =0点的误差最小，离0越远，误差越大。
讲了这么多，帕德逼近怎么求呢？切比雪夫近似值的求法是很麻烦的，要用到复杂的雷米兹算法，但是帕德逼近的求法很简单。我拿个例子来说吧。
例子
以求f ( x ) = s i n ( x ) f(x)=sin(x)f (x )=s in (x )的[ 2 / 2 ] [2/2][2/2 ]阶帕德逼近为例子，我讲讲帕德逼近的过程。因为分母有2阶，而分子也有2阶，所以总体复杂度有4阶。而s i n ( x ) sin(x)s in (x )是一个很”无理”的函数，所以要让他”有理”。那么有两种办法，第一种办法是泰勒级数，第二种办法是切比雪夫近似值。因为切比雪夫近似值需要有定义域，所以直接排除，那么只剩下泰勒级数了。那么我们就展开到4阶吧。
s i n ( x ) = x − x 3 3 ! + O ( x 5 ) sin(x) = x-\frac{x^3}{3!}+O(x^5)s in (x )=x −3 !x 3 +O (x 5 )
它的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}R 2 ,2 (x )=b (x )a (x )=1 +b 1 x +b 2 x 2 a 0 +a 1 x +a 2 x 2
要求泰勒级数去掉尾项后减去帕德逼近要约等于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 !x 3 −1 +b 1 x +b 2 x 2 a 0 +a 1 x +a 2 x 2 ≈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 !x 3 =1 +b 1 x +b 2 x 2 a 0 +a 1 x +a 2 x 2 (x −3 !x 3 )(1 +b 1 x +b 2 x 2 )=a 0 +a 1 x +a 2 x 2
将左边展开：
( 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 !x 3 )(1 +b 1 x +b 2 x 2 )=x −6 1 x 3 +b 1 x 2 −6 1 b 1 x 4 +b 2 x 3 −6 1 b 2 x 5 +=x +b 1 x 2 +(b 2 −6 1 )x 3 −6 1 b 1 x 4 −6 1 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 = 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 +b 1 x 2 +(b 2 −6 1 )x 3 −6 1 b 1 x 4 −6 1 b 2 x 5 =a 0 +a 1 x +a 2 x 2
左边减去右边，有：
− 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 −a 0 +(1 −a 1 )x +(b 1 −a 2 )x 2 +(b 2 −6 1 )x 3 −6 1 b 1 x 4 −6 1 b 2 x 5 =0
把最高次的x 5 x^5 x 5项忽略，剩余的每项等于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\−a 0 =0 1 −a 1 =0 b 1 −a 2 =0 b 2 −6 1 =0 −6 1 b 1 =0
把上面的方程组解开得：
b 1 = 0 b 2 = 1 6 b_1=0\ b_2=\frac1{6}b 1 =0 b 2 =6 1
代入得：
a 0 = 0 a 1 = 1 a 2 = 0 a_0=0\ a_1=1\ a_2=0 a 0 =0 a 1 =1 a 2 =0
于是最终s i n ( x ) sin(x)s in (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}R 2 ,2 (x )=1 +6 1 x 2 x
在坐标系上，帕德逼近和原函数在越接近于0的地方，误差越小，离0越远，就误差越大，越来越离，如R 2 , 2 R_{2,2}R 2 ,2 与s i n ( x ) sin(x)s in (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 )=c 0 +c 1 x +c x x 2 +⋯+c m +n x m +n +O (x m +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}⎣⎡c n ⋮c 2 n −1 …⋱⋯c 1 ⋮c n ⎦⎤⎣⎡b 1 ⋮b n ⎦⎤=⎣⎡−c n +1 ⋮−c 2 n ⎦⎤
而b 0 = 1 b_0=1 b 0 =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 }a 0 =c 0 a 1 =c 0 b 1 +c 1 ⋮a m =i =0 ∑m c i b m −i
python代码
python代码就是根据泰勒级数求帕德逼近的m + n + 1 m+n+1 m +n +1个未知参数。当然，根据方程求麦克劳林级数（在x = 0 x=0 x =0点的泰勒级数）的代码我是不会写的，因为场景太多太复杂。对于解方程，我这里复制了一个简单的矩阵类：
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):
y[i] = values[i]
for j in range(i + 1, n):
values[j] -= l.lines[j][i] * values[i]
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):
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]
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)])
def lup_decomposition(self):
n = len(self.__lines)
l = copy.deepcopy(self.__lines)
u = copy.deepcopy(self.__lines)
a = copy.deepcopy(self.__lines)
p = [i for i in range(0, n)]
for i in range(0, n):
swap_line = i
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]
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
for k in range(i + 1, n):
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()
pb = p * Matrix([[num] for num in values])
pb_values = [i[0] for i in pb.__lines]
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 3,3 3 ,3阶为例子，进行测试。首先计算sin ( x ) \sin(x)sin (x )的麦克劳林级数：
s i n ( x ) = x − x 3 3 ! + x 5 5 ! + O ( x 7 ) sin(x)=x-\frac{x^3}{3!}+\frac{x^5}{5!}+O(x^7)s in (x )=x −3 !x 3 +5 !x 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]
\draw (-4,0) [->, thick]
% 画X轴
\foreach \x in {-4,-3,-2,-1,1,2,3,4}{
\node at (\x,0) [below] {\x};
\draw (\x,0)
}
\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)
}
\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 2/2 2/2阶的帕德逼近精确多了：
Original: https://blog.csdn.net/m0_66201040/article/details/125866121
Author: 醒过来摸鱼
Title: 8.4 帕德逼近
原创文章受到原创版权保护。转载请注明出处：https://www.johngo689.com/623451/
转载文章受原作者版权保护。转载请注明原作者出处！