【洛谷日报 #53】浅谈一些求近似值的方法

发表于 2018-08-09 更新于 2024-07-02 分类于算法竞赛本文字数： 2.5k 阅读时长 ≈ 9 分钟

~~只是为了拿本子随便写的~~

老文章，可能有很多错误，懒得修了

前言

对于低于 5 次的多项式方程，我们有通用的公式解法求零点的精确值

对于一些特殊高次多项式方程 (例如可以因式分解的或者满足一些特定形式的方程) 的和一些特殊的超越方程，我们也有方法求零点的精确值

但是其余的情况呢？

目前来说我们只能求近似值 QwQ (而且在实际应用中，精确值往往也会被转换成近似值)

下面简要介绍几种方法

二分法

可以说相当常见了

对于在区间 \([l,r]\) 内单调、连续且有 \(f(l)\cdot f(r)<0\) 成立的 \(f(x)\), 做如下操作:

计算 \(m=\frac{l+r}{2}\)
若 \(f(l)f(m)<0\), 则令 \(r=m\), 否则令 \(l=m\)
如果达到预定精度，跳转到 4, 否则跳转到 1
结束

循环次数:

\[ \left\lceil\log_2 {\frac{r-l}{\epsilon}}\right\rceil \]

附程序:

binary_search.cppview raw

double binary_search(double l, double r, double EPS) {
  double mid;
  while (r - l > EPS) {
    mid = l + (r - l) / 2;  // 防溢出
    if (F(mid) * F(l) <= 0) r = mid;
    else l = mid;
  }
  return mid;
}

0.618 法

常用于求单峰函数最值

先证明一下它的最优性 (摘自人教版高中数学选修 4-7 ~~没错真的有这本选修~~)

为了使每次去掉的区间有一定的规律性，我们这样来考虑: 每次舍去的区间占舍去前的区间的比例数相同

下面进一步分析如何按上述两个原则确定合适的试点。如图 2-1 设第 1 试点、第 2 试点分别为 \(x_1\) 和 \(x_2\), \(x_2<x_1\) 且 \(x_1,x_2\) 关于 \([a,b]\) 的中心对称，即 \(x_2-a=b-x_1\)

(图 2-1, 由 GeoGebra 生成)

显然，不论点 \(x_2\)(或点 \(x_1\)) 是好点还是差点，由对称性，舍去的区间长度都等于 \(b-x_1\), 不妨设 \(x_2\) 是好点，\(x_1\) 是差点，于是舍去 \((x_1,b]\). 再在存优范围 \([a,x_1]\) 内安排第 3 次试验，设试点为 \(x_3\), \(x_3\) 与 \(x_2\) 关于 \([a,x_1]\) 的中心对称 (如图 2-2 所示)

(图 2-2, 由 GeoGebra 生成)

点 \(x_3\) 应在点 \(x_2\) 左侧，因为如果点 \(x_3\) 在点 \(x_2\) 的右侧，那么当 \(x_3\) 是好点，\(x_2\) 是差点时，要舍去区间 \([a,x_2]\), 而它的长度与上次舍去的区间 \((x_1,b]\) 的长度相同，违背成比例舍去的原则。于是，不论点 \(x_3\) (或点 \(x_2\)) 是好点还是差点，被舍去的区间长度都等于 \(x_1-x_2\), 按成比例舍去的原则，我们有等式

\[ \frac{b-x_1}{b-a}=\frac{x_1-x_2}{x_1-a}\tag{1} \]

其中，左边是第一次舍去的比例数，右边是第二次舍去的比例数，对式 (1) 变形，得

\[ 1-\frac{b-x_1}{b-a}=1-\frac{x_1-x_2}{x_1-a} \]

即

\[ \frac{x_1-a}{b-a}=\frac{x_2-a}{x_1-a}\tag{2} \]

式 (2) 两边分别是两次舍弁后的存优范围占舍弃前全区间的比例数，设每次舍弃后的存优范围占舍弃前全区间的比例数为 \(t\), 即

\[ \frac{x_1-a}{b-a}=t\tag{3} \]

则由 \(b-x_2=x_1-a\) 可得

\[ \frac{x_2-a}{b-a}=1-t\tag{4} \]

由式 (2) 得

\[ \frac{x_1-a}{b-a}=\frac{\frac{x_2-a}{b-a}}{\frac{x_1-a}{b-a}} \]

把 (3) 与 (4) 代入 (5), 得

\[ t=\frac{1-t}{t} \]

即

\[ t^2+t-1=0 \]

解得 \(t_1=\displaystyle\frac{-1+\sqrt{5}}{2}, t_2=\frac{-1-\sqrt{5}}{2}\), 其中 \(t_1\) 为对本问题有意义的根，这就是黄金分割常数 , 用 \(\varphi\) 表示 (注：原文用 \(\omega\) 表示)

一句话概括就是在缩小区间后可以只计算一个试点坐标，从而保证最优

流程如下

计算 \(m_1=l\varphi+r(1-\varphi)\), \(m_2=l(1-\varphi)+r\varphi\)
若 \(f(l) f(m_1)>0\)
1. 则令 \(l=m_1,m_1=m_2,m_2=l(1-\varphi)+r\varphi\)
2. 否则令 \(r=m_2,m_2=m_1,m_1=l\varphi+r(1-\varphi)\)
如果达到预定精度，跳转到 4, 否则跳转到 2 (注意这里跳转到 2)
结束

附程序:

gold_search.cppview raw

const double PHI = 0.61803399, mPHI = 0.38196601;
double gold_search(double l, double r, double EPS) {
  double mid1 = l + mPHI * (r - l), mid2 = l + PHI * (r - l);
  while (r - l > EPS) {
    if (F(mid1) < F(mid2)) {
      l = mid1;
      mid1 = mid2;
      mid2 = l + PHI * (r - l);
    } else {
      r = mid2;
      mid2 = mid1;
      mid1 = l + mPHI * (r - l);
    }
  }
  return (mid1 + mid2) / 2;
}

读者们可以在洛谷 P3382 中测试一下 (～ o ￣ 3 ￣)～

关于这个还有一个类似方法：斐波那契法。有兴趣的读者可以查阅相关资料~~才不是笔者不想写~~_(:3」∠)_

Taylor 公式

先讲这个是为了为下文 Newton 迭代法二次收敛的证明做铺垫，不想看证明的可以略过 QwQ (不过还是推荐了解一下，挺有趣的)

这里假定函数 \(f(x)\) 在 \(x_0\) 处有任意阶导数

我们可以很容易地求出多项式和类指数函数的近似值，但是像三角函数、对数函数这样的我们又该如何求近似值呢

对了，就是用 Taylor 公式 QwQ

Taylor 公式的想法很简单，就是构造一个多项式函数 \(g(x)=\displaystyle\sum_{k=0}^n{a_kx^k}\), 使得它与函数 \(f(x)\) 在 \(x_0\) 处的原函数值和各阶导数均相等，即

\[ \begin{aligned} f(x_0)&=g(x_0)\\ f^\prime(x_0)&=g^\prime(x_0)\\ f^{\prime\prime}(x_0)&=g^{\prime\prime}(x_0)\\ &...\\ f^{(n)}(x_0)&=g^{(n)}(x_0) \end{aligned} \]

因为

\[ g^{(m)}(x)=\sum_{k=m}^n{\frac{k!}{(k-m)!}a_kx^{k-m}} \]

于是便有

\[ g(x)=\sum_{k=0}^n{\frac{f^{(k)}(x_0)}{k!}(x-x_0)^k} \]

当 \(n\rightarrow\infty\) 时，我们可以认为 \(f(x)=g(x)\)

而当 \(n\) 有一个确定的值时，\(f(x)\) 就可以写成 \(g(x)+R_n(x)\) 了

其中 \(R_n(x)\) 是余项，它有好几种不同的写法，比如 Lagrange 余项

\[ R_k(x)=\frac{f^{(k+1)}(\xi_L)}{(k+1)!}(x-x_0)^{k+1} \]

其中 \(\xi_L\) 在 \(x\) 和 \(x_0\) 之间

当 \(n\rightarrow\infty\) 时，有 (Taylor 级数)

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

特别地，当 \(x_0=0\) 时，有 (Maclaurin 级数)

\[ \sum_{k=0}^\infty{\frac{f^{(k)}(0)}{k!}}x^k \]

另外注意应用 Maclaurin 级数并且 \(x\) 在某个范围之外时，得到的结果可能是发散的 (这个不展开讲，有兴趣的读者可以去学习无穷级数相关知识)

附上 Wikipedia 的动图

对证明感兴趣的读者可以自行查阅相关资料

下面给出几个常见的 Taylor 级数

\[ e^x=\sum_{k=0}^\infty{\frac{x^k}{k!}} \]

\[ \sin x=\sum_{k=0}^\infty{(-1)^k\frac{x^{2k+1}}{(2k+1)!}} \]

\[ \cos x=\sum_{k=0}^\infty{(-1)^k\frac{x^{2k}}{(2k)!}} \]

(有上面三个式子就可以证明欧拉公式之 \(e^{i\theta}=\cos\theta+i\sin\theta\) 了)

\[ \ln{(1+x)}=\sum_{k=1}^\infty{(-1)^{k+1}\frac{x^k}{k}} \]

\[ \frac{1}{1-x}=\sum_{k=0}^\infty{x^k} \]

\[ (1+x)^m=\sum_{k=0}^\infty{\binom{m}{k}x^k} \]

Newton 迭代法

先说说过程

随便确定一个数 \(x_0\)
求在 \(f(x_0)\) 处的切线 \(l:[y-f(x_0)]=f^\prime(x_0)(x-x_0)\)
求切线 \(l\) 的零点 \(x_1\)

稍加计算便得到了

\[ x_1=x_0-\frac{f(x_0)}{f^\prime(x_0)} \]

既然是迭代，那么自然就有

\[ x_{n+1}=x_n-\frac{f(x_n)}{f^\prime(x_n)} \]

其中 \(x_n\) 代表第 \(n\) 次迭代

附上 Wikipedia 的动图

二次收敛证明: (Wikipedia 上的，笔者翻译 QwQ)

根据 Taylor's theorem, 任何二阶导数连续的函数 \(f(x)\) (设 \(\alpha\) 是根) 都可以写成

\[ f(\alpha)=f(x_n)+f^\prime(x_n)(\alpha-x_n)+R_1\tag{1} \]

由 Lagrange form of the Taylor series expansion remainder 得

\[ R_1=\frac{1}{2!}f^{\prime\prime}(\xi_n)(\alpha-x_n)^2 \]

其中 \(\xi_n\) 在 \(x_n\) 和 \(\alpha\) 之间

由于 \(\alpha\) 是根，所以 (1) 式变为

\[ 0=f(\alpha)=f(x_n)+f^\prime(x_n)(\alpha-x_n)+\frac{1}{2}f^{\prime\prime}(\xi_n)(\alpha-x_n)^2\tag{2} \]

(2) 式两边同时除以 \(f^\prime(x_n)\), 整理得

\[ \frac{f(x_n)}{f^\prime(x_n)}+(\alpha-x_n)=\frac{-f^{\prime\prime}(\xi_n)}{2f^\prime(x_n)}(\alpha-x_n)^2\tag{3} \]

由于

\[ x_{n+1}=x_n-\frac{f(x_n)}{f^\prime(x_n)}\tag{4} \]

代入 (3) 式，有

\[ \underbrace{\alpha-x_{n+1}}_{\epsilon_{n+1}}=\frac{-f^{\prime\prime}(\xi_n)}{2f^\prime(x_n)}(\underbrace{\alpha-x_n}_{\epsilon_n})^2 \]

即

\[ \epsilon_{n+1}=\frac{-f^{\prime\prime}(\xi_n)}{2f^\prime(x_n)}\epsilon_n^2\tag{5} \]

两边取绝对值，有

\[ |\epsilon_{n+1}|=\frac{|f^{\prime\prime}(\xi_n)|}{2|f^\prime(x_n)|}\epsilon_n^2\tag{6} \]

(6) 式表明，如果函数满足以下条件，其为二次收敛

对于所有的 \(x\in I\) (\(I\) 为区间 \([\alpha-r,\alpha+r]\),\(r\geq |\alpha-x_0|\) (即 \(x_0\in I\)) ) , 有 \(f^\prime(x)\neq 0\)

对于所有的 \(x\in I\), \(f^{\prime\prime}(x)\) 连续

\(x_0\) 足够接近根

"足够接近" 意为

Taylor 近似足够准确 (可以忽略高阶项)

对于 \(C<\infty\)(原文如此), \(\displaystyle\frac{1}{2}\Big|\frac{f^{\prime\prime}(x_n)}{f^\prime(x_n)}\Big|<C\Big|\frac{f^{\prime\prime}(\alpha)}{f^\prime(\alpha)}\Big|\)

对于 \(n\in\Bbb{N}\), \(\displaystyle C\Big|\frac{f^{\prime\prime}(\alpha)}{f^\prime(\alpha)}\Big|\epsilon_n<1\)

当满足上述条件时，(6) 式可以写为:

\[ |\epsilon_{n+1}|\leq M\epsilon_n^2 \]

其中

\[ M=\sup_{x\in I}{\frac{1}{2}\bigg|\frac{f^{\prime\prime}(x)}{f^\prime(x)}\bigg|} \]

由条件 3 得 \(M|\epsilon_0|<1\)

程序:

newton.cppview raw

double newton(double x0, int n) {
  for (int i = 0; i < n; ++i) x0 -= f(x0) / df(x0);
  return x0;
}

当然，上述各方法的应用范围远不止于此，有兴趣的读者可以自行查阅相关资料 QwQ

赠品

\(\cos x=x\) 的解析解

对于 PhOer 来说，\(\cos x=x\) 这个方程应该是相当熟悉了 QwQ

笔者在这里放上解析解 (近似值 \(x=0.739\)) , 详情见参考文献 ¹ (文献里讲的是 \(t\sin x=x-m\) 的解法，不过笔者太弱了，实在是看不懂 QwQ)

\[ \frac{\pi}{2}\exp\left(\frac{1}{\pi}\int_0^1\frac{\arctan\left({(\pi x+2)\log\left(\frac{\sqrt{1-x^2}+1}{x}\right)x\over x^2\log^2\left(\frac{\sqrt{1-x^2}+1}{x}\right)-\pi x-1}\right)}{x}\mathrm{d}x\right) \]

快速求 \({1\over\sqrt{x}}\)

关于这个有一个相当有名的故事: 一个 Sqrt 函数引发的血案 (这是笔者能找到的最早一篇了 QwQ)

测试程序:

test.cppview raw

#include <bits/stdc++.h>

float InvSqrt(float x, int it) {
  float xhalf = 0.5f * x;
  int i = *(int *)&x;

  i = 0x5f3759df - (i >> 1);
  x = *(float *)&i;
  for (int i = 1; i <= it; ++i) x = x * (1.5f - xhalf * x * x);
  return x;
}
float GetRErr(float test, float real) {
  return fabs((test - real) / real * 100.0);
}
int main() {
  int i = 0;
  while (1) {
    float t;
    scanf("%f", &t);
    if (t < 0) return 0;
    float sys = 1.0 / sqrt(t), invs1 = InvSqrt(t, 1), invs2 = InvSqrt(t, 2),
          invs3 = InvSqrt(t, 3), invs4 = InvSqrt(t, 4), invs5 = InvSqrt(t, 5),
          invs6 = InvSqrt(t, 6);
    printf(
      "Test#%d n=%f\n"
      "System:\n%f\nTest:\n"
      "it#1:%f RelativeErr:%.2f%%\n"
      "it#2:%f RelativeErr:%.2f%%\n"
      "it#3:%f RelativeErr:%.2f%%\n"
      "it#4:%f RelativeErr:%.2f%%\n"
      "it#5:%f RelativeErr:%.2f%%\n"
      "it#6:%f RelativeErr:%.2f%%\n\n",
      ++i,
      t,
      sys,
      invs1,
      GetRErr(invs1, sys),
      invs2,
      GetRErr(invs2, sys),
      invs3,
      GetRErr(invs3, sys),
      invs4,
      GetRErr(invs4, sys),
      invs5,
      GetRErr(invs5, sys),
      invs6,
      GetRErr(invs6, sys));
  }
  return 0;
}

后记

这篇文章偏数学一些，如果不能理解的话请多读几遍 QwQ

其实可写的还有很多，限于篇幅就到此为止了 (现在在后台打个字都要卡两秒 QwQ)

因为笔者是个蒟蒻，所以如果有错误，烦请各位 dalao 不吝赐教

主要参考资料

人教版高中数学选修 4-7 优选法与试验设计初步
Newton's method - Wikipedia
Taylor's theorem - Wikipedia

Siewert C E, Burniston E E. An exact analytical solution of Kepler's equation[J]. Celestial Mechanics, 1972, 6(3):294-304.↩︎