使用泰勒公式的马华算法

mokongbro · 2025-08-12 16:59:41 · 算法·理论

来源（）

积分优化

根据麦克劳林公式（假设 f 是一个解析函数）：

直接逐项积分得到：

所以，由 \displaystyle\int_{a}^{b}f(x)\mathrm{d}x=\int_{a}^{c}f(x)\mathrm{d}x+\int_{c}^{b}f(x)\mathrm{d}x 我们有：

我们发现这个算法竟然远优于定积分的定义式：

对于每一份输入造就的唯一 f(x) 来说，我们只需预处理出 f^{(n)}(0) 以及 n! 的数组，

就可以处理几个求导运算增长慢于 \displaystyle\int_{0}^{\infty}e^{-x}x^n\mathrm{d}x 的函数的积分，这时最多计算几百项就可以逼近积分值，但是 b^{n+1},a^{n+1} 会非常大（如果您先前的 a 和 b 都是大约 \color{green}27 这种数的话）,可能会把 \text{long double} 给爆掉。

被吓傻了（好像我罪不至此）。

可以尝试去优化一下阶乘和幂次，使得他们互相抵消，并且误差在可以接受的范围内，如记录数组

由于我们可以使求和的顺序使得程序沿 i 递增的方向求和，所以我们大致不必记录 2 个数组这么夸张，只需要在更新值的时候用

an=an*a/(i+1)
bn=bn*b/(i+1)
ans+=df0[i]*(bn-an)

那么它们不仅不会过于巨大，反而会在一个模糊的边界之后迅速变小，我们在察觉到这一点后就可以准备 \text{return} 我们的答案。

微分预备

微分预备，就是要求一个函数 f 在原点处的 n 阶导数。

众所周知，导数的定义是

而多阶导数则是递归定义的，这就是：

这就启发我们使用一个函数来封装它的这么多阶导数，由于我们只关心在原点的导数，所以它确实可以实现记忆化。

例如，我们定义一个 \text{Deri} 函数，让它直接返回 f 的 n 阶导函数，按照上面的定义，我们可以写出如下片段：

DX=int(3e10)
def Deri(f,n:int,dx=None):
    if dx==None:
        dx=DX
    else:
        dx=int(dx)
    if n<=0:
        return f
    return (lambda x:(Deri(f,n-1,dx)(x+1/dx)-Deri(f,n-1,dx)(x))*dx)

我说过，这个竟然适合加记忆化，我个人喜欢这么写。

from functools import wraps
import sys
from decimal import Decimal as dec
sys。setrecursionlimit(0x7fffffff)
def mem(func):
    cache={}
    @wraps(func)
    def wrapper(*args,**kwargs):
        key=args,tuple(sorted(kwargs。items()))
        if key not in cache:
            cache[key]=func(*args,**kwargs)
        return cache[key]
    return wrapper
DX=int(3e10)
@mem
def Deri(f,n:int,dx=None):
    if dx==None:
        dx=DX
    else:
        dx=int(dx)
    if n<=0:
        return f
    return (lambda x:(Deri(f,n-1,dx)(x+1/dx)-Deri(f,n-1,dx)(x))*dx)

优化 \text{DX} 和 \text{dx} 的值，便可以控制这个函数的精度，这里我犯唐了，把这个增量写成倒数了；非常令人卑微与困惑的是，这里的误差会逐层累加，但是在 0 左右，我们可以尝试控制 \text{DX} 的值，在爆精度之间苦苦挣扎。

小谈一手，之前我们通过分部积分法证明了：

但是这里的导数需要考虑一个变量 x ,但是 x+\text{dx} 的精度就没有 \text{dx} 本身那么高，所以这个只是 O.E.L 在数学部分的一个废案。

数值积分

所以，我们经过了上一个部分的准备，我们可以求出 f^{(n)}(0) 了，我们现在可以考虑以下式子：

更近一步地，我们可以得到：

所以，对于一个存在 (m+1) 阶导数的函数，满足：

我们代码的主要内容就是求出前面的求和式，而截断后续的积分运算。

采用（将这个积分代入上上式并展开）：

所以我们就可以用如下代码（假设我们已经处理了数组 \text{df0[M]} 存储 f 在原点处的各阶导数）。

#define ll long long
#define ld long double
ld *df0;
ll DX=3e10;
ld d_int0(auto f,ld x,ll dx=DX){
    ld ans=0,xn=1;
    ll i=0;
    while(df0[i]*xn>1.0l/dx){
        xn=xn*x/(i+1);
        ans+=df0[i++]*xn;
    }return ans;
}

那我这个传入的 f 有什么用啊。

求反函数

整个马华的初衷就是求出一个 x 使得一个方程 I(x)=0，只不过这里 \displaystyle I(x)=\int_{a}^{x}f(t)\mathrm{d}t-s 包含一个积分。

由于我们的积分运算太耗时间了，所以枚举+二分的想法太慢，由于我们可以简单求出一个函数 f 在一个点 x=x_0 的切线 y=f(x_0)+f'(x_0)(x-x_0),

我们可以求出这个切线的零点 \displaystyle x=x_0-\frac{f(x_0)}{f'(x_0)}。

做一个序列

我们来讨论，如果 \displaystyle \lim_{n\to\infty}x_n 存在，那么

考虑当 \displaystyle \lim_{n\to\infty}x_n 存在时，有 A\in\mathbb{R} 满足

\forall\varepsilon>0\exists N\in\mathbb{N}_+\forall n\in\mathbb{N}\wedge n>N\left(|x_n-A|<\varepsilon\right)

找出任意两个数 n,m>N 则有 |x_n-A|<\varepsilon\wedge|x_m-A|<\varepsilon,即

x_n,x_m\in\left(A-\varepsilon,A+\varepsilon\right)

所以

|x_n-x_m|<2\varepsilon

因为 2>0,而 \varepsilon\ne0,所以这个式子实际上等价于

\forall x\in\mathbb{R}_+(|x_n-x_m|<x)

我们知道，说一个数小于集合中所有的数，那么他一定小于等于这个集合的下确界（下确界的定义），于是就可以把这个式子写成

|x_n-x_m|\le\inf\mathbb{R}_+

而绝对值具有非负性，于是我们证明了：

\lim_{n\to\infty,m\to\infty}x_n-x_m=0

取 m=n+1 代入到递推式中，我们得到了

\lim_{n\to\infty}\frac{f(x_n)}{f'(x_n)}=0

现在可以大致分为 2 种情况

1. 我们找到了函数的零点，\displaystyle f\left(\lim_{n\to\infty}x_n\right)=0；

2. 这个可爱的导数太巨大了，\displaystyle \lim_{n\to\infty}f'(x_n)=\infty。

   但如果真的是这样，那么函数在 A 处不可导，并且

   \lim_{\Delta x\to0}\frac{f(A+\Delta x)-f(A)}{\Delta x}=\infty

   说明在 A 的某个小邻域内有至少 f(x_1)-f(x_2)=\mathrm{O}(1),那么在这个区间内密集地取无穷个点，我们得到了 f(A) 的数量级应该至少是 \mathrm{O}(\infty)：

   \lim_{x\to A}f(x)=\infty

   而且他是比 f'(x) 高阶的无穷大，那既然是高阶的无穷大，应该有 \displaystyle\lim_{n\to\infty}\frac{f(x_n)}{f'(x_n)}=\infty 与假设矛盾。

所以 \displaystyle f\left(\lim_{n\to\infty}x_n\right)=0。

那么有的同学就要问了，既然 x_{n+1}=\displaystyle x_n-\frac{f(x_n)}{f'(x_n)} 这种迭代可以，那么是不是我用 x_{n+1}=x_n-f(x_n) 也可以呢？

然后你就会发现：

如此循环三十年，直到大厦崩塌；

一万个无穷大啊，在这数组中奔跑。

没错，对于数列极限存在的情况下，后者的证明更简单，但是后者的极限存在的概率有点碰运气。

比如，你要求一个 f(x)=x^2 的零点，从 x_0=-1 开始跑数列，结果他猛地就朝着负无穷的远方飞去了。

那我们所说的这个方法呢？是否对于一个存在零点的解析函数 f 我们的数列一定存在一个极限 A 呢？

这个方法不一定能够在整个数轴上使用，你可能会遇到数值不存在、除以零错误等令人头疼的意外终止。

除开这里，他一定在某个循环之后陷入一个漩涡中，以逼近零点（否则落入到终止）。

1. 假设 f'(x)>0\wedge f''(x)<0。

   当 x_n 终于跑到了 A 的左侧去，因为 f'(x)>0 这时的 f(x_n)<0,判断出来 \displaystyle \frac{f(x_n)}{f'(x_n)}<0,

   那么 x_{n+1}>x_n,在这一步为递增，因为 f'(x_n)>0,所以 f(x_{n+1})>f(x_n)。

   因为 f''(x_n)<0,所以 f'(x_{n+1})<f'(x_n),

   在 x>x_n 的时候，对于切线 g(x)=f(x_n)+f'(x_n)(x-x_n) 有 g(x)>f(x)。

   于是，x_{n+1}<A。

   这就是说，如果在本假设下，出现了 x_n<A 那么 x_n<x_{n+1}<A。

   单调递增有上界序列必有极限。

2. 假设 f'(x)>0\wedge f''(x)>0。

   重复上例的证明，只不过需要换方向。

   单调递减有下界数列必有极限。

3. 假设 f'(x)\ne0\wedge f''(x)=0。

   神秘一次函数。

剩余的情况可以用取相反数变换变为上述的第 1 或第 2 种情况。

单调有界数列必有极限。

这个方法叫做牛顿迭代法。当我们递推数列的时候遇见了神秘 \color{red}\text{OVERFLOW} 的时候，要么就是函数根本没有零点，要么就是函数的零点在远方，由于我们十分自大，所以这些情况一律视为没有零点。

mh函数

在牛顿迭代法中，令 \displaystyle I(x)=\int_{a}^{x}f(t)\mathrm{d}t-s 得到我们的数列：

于是就有以下的写法。

def mh(f,a,s,dx=None):
    if dx==None:
        dx=DX
    else:
        dx=int(dx)
    xx=a
    i=int(0)
    try:
        while f(xx)==0:
            xx=xx+1
        tmp=(d_int(f,a,xx,dx)-s)/f(xx)
        while abs(tmp)>1/dx:
            xx-=tmp
            tmp=(d_int(f,a,xx,dx)-s)/f(xx)
            i+=1
    except ZeroDivisionError as e:
        print("--------------------------")
        print(f"There's An ERROR:{e} ,so it's {xx} as answer give you")
        print("--------------------------")
        return xx
    except Exception as e:
        print("--------------------------")
        print(f"There's An ERROR:{e} ,\nso it's a piece of shit for you as {xx}")
        print("--------------------------")
        return xx
    return xx

选择 a 作为数列初始值是因为我们想要求的是距离 a 最近的解。

所以完整的 \text{py} 代码是这样的。

from functools import wraps
import sys
from decimal import Decimal as dec
sys。setrecursionlimit(0x7fffffff)
def mem(func):
    cache={}
    @wraps(func)
    def wrapper(*args,**kwargs):
        key=args,tuple(sorted(kwargs。items()))
        if key not in cache:
            cache[key]=func(*args,**kwargs)
        return cache[key]
    return wrapper
DX=int(3e10)
@mem
def Deri(f,n:int,dx=None):
    if dx==None:
        dx=DX
    else:
        dx=int(dx)
    if n<=0:
        return f
    return (lambda x:(Deri(f,n-1,dx)(x+1/dx)-Deri(f,n-1,dx)(x))*dx)
def d_int0(f,x,dx=None):
    if x==0:
        return 0
    if dx==None:
        dx=DX
    else:
        dx=int(dx)
    ans=0
    xn=1
    i=int(0)
    while abs(2*Deri(f,i,dx)(0)*xn)>1/dx or i<5:
        ans+=Deri(f,i,dx)(0)*xn
        xn=xn*x/(i+1)
        i+=1
    return ans
def d_int(f,a,b,dx=None):
    if dx==None:
        dx=DX
    else:
        dx=int(dx)
    return d_int0(f,b,dx)-d_int0(f,a,dx)
def mh(f,a,s,dx=None):
    if dx==None:
        dx=DX
    else:
        dx=int(dx)
    xx=a
    i=int(0)
    try:
        while f(xx)==0:
            xx=xx+1
        tmp=(d_int(f,a,xx,dx)-s)/f(xx)
        while abs(tmp)>1/dx:
            xx-=tmp
            tmp=(d_int(f,a,xx,dx)-s)/f(xx)
            i+=1
    except ZeroDivisionError as e:
        print("--------------------------")
        print(f"There's An ERROR:{e} ,so it's {xx} as answer give you")
        print("--------------------------")
        return xx
    except Exception as e:
        print("--------------------------")
        print(f"There's An ERROR:{e} ,\nso it's a piece of shit for you as {xx}")
        print("--------------------------")
        return xx
    return xx

其他细节

1. 切线

   我们知道一根直线的表达式就是

   y=kx+b

   切线寻找在点 x=x_0 处 y'=f'\wedge y=f 的直线 y=kx+b。

   代入导数的公式我们算出来这里 y'=k,所以 k=f'(x_0)。

   于是有方程

   f(x_0)=f'(x_0)x_0+b b=f(x_0)-f'(x_0)

   将它们代回直线表达式得到

   y=f'(x_0)x+f(x_0)-f'(x_0)

   整理得

   y=f(x_0)+f'(x_0)(x-x_0)

   这就是切线方程。

2. 分部积分法

   分部积分法用于处理函数乘积的积分。

   公式：

   \int f\mathrm{d}g=fg-\int g\mathrm{d}f

   利用不定积分的定义，对两边同时求导

   fg'=(f'g+fg')-fg'

   这是恒成立的。

   这得自于函数相乘的求导法则

   (fg)'=f'g+fg'

   利用定积分的定义，我们发现定积分的微分就是函数本身乘以自变量的微分。

   \int_{a}^{b}f(x)\mathrm{d}x=\lim_{\max\{\Delta x_i\}\to0}\sum f(x_i)\Delta x_i,(a\le x_0<\cdots\le b)

   所以

   \int_{a}^{x+\Delta x}f(t)\mathrm{d}t-\int_{a}^{x}f(t)\mathrm{d}t=\int_{x}^{x+\Delta x}f(t)\mathrm{d}t\to f(x)\Delta x

   因此分布积分法对于定积分适用。

3. 带有积分余项的泰勒公式

   利用定积分的上述性质，我们可以得出

   \int_{x_0}^{x}f'(t)\mathrm{d}t=f(x)-f(x_0)

   等号左侧的 \mathrm{d}t 可以看成分部积分法中的 \mathrm{d}g 来进行分部积分操作。

   因为

   \mathrm{d}(x-t)=-\mathrm{d}t

   从而为了导出泰勒公式的形态，我们将这个积分重新写作

   -\int_{x_0}^{x}f'(t)\mathrm{d}(x-t)

   套用分部积分法的公式，我们得到他等于

   -\left[f'(t)(x-t)\right]_{x_0}^{x}+\int_{x_0}^{x}(x-t)\mathrm{d}(f'(t))

   我们先来处理左侧的方括号。

   被减数：

   f'(x)(x-x)=0

   减数：

   f'(x_0)(x-x_0)

   用被减数减去减数再相反数得到

   f(x)-f(x_0)=f'(x_0)(x-x_0)+\int_{x_0}^{x}f''(t)(x-t)\mathrm{d}t

   于是这成为我们数学归纳法的起始项。

   于是，求通用的积分

   \int_{x_0}^{x}f^{(n+1)}(t)(x-t)^n\mathrm{d}t= -\frac{1}{n+1}\int_{x_0}^xf^{(n+1)}(t)\mathrm{d}\left((x-t)^{n+1}\right)= -\frac{1}{n+1}\left(\left[f^{(n+1)}(x-t)^{n+1}\right]_{x_0}^x-\int_{x_0}^{x}(x-t)^{n+1}\mathrm{d}\left(f^{(n+1)}(t)\right)\right)= \frac{\displaystyle f^{(n+1)}(x_0)}{n+1}(x-x_0)^{n+1}+\frac{1}{n+1}\int_{x_0}^{x}f^{(n+2)}(x-t)^{n+1}\mathrm{d}t

   我们又有 \displaystyle \prod_{i=1}^{n}i:=n! 作为 n 的阶乘，直接就有

   f(x)=\sum_{m=0}^{n}\frac{f^{(m)}(x_0)}{m!}(x-x_0)^m+\frac{1}{n!}\int_{x_0}^{x}f^{(n+1)}(t)(x-t)^{n}\mathrm{d}t

   另外，如果函数 f 满足

   f(x)=\sum_{n=0}^{\infty}\frac{f^{(n)}(x_0)}{n!}(x-x_0)^n

   那么我们就称 f 为解析函数。

   泰勒公式在 x_0=0 处的特例叫做麦克劳林公式。

4. 如此循环三十年，直到大厦崩塌；

   一万个无穷大啊，在这数组中奔跑。

   原歌词：（万能青年旅店《杀死那个石家庄人》）

   如此生活三十年，直到大厦崩塌；

   一万匹脱缰的马，在他脑海中奔跑。

   对不起万能青年旅店，对不起董二千、机耕、小耕、冯江、以及石家庄全体人民。