使用 MATLAB 求三角函数的 Minimax 多项式逼近

Question

我正在尝试使用 MATLAB 中的 remez 交换算法找到正弦和余弦的极小极大多项式近似值。需要精确到 23 位，因为我正在为 IEEE-754 浮点实现正弦和余弦函数。

使用此链接here（请参阅第 8 页至第 15 页），给出了使用 Mathematica 和 Maple 查找多项式的说明，但是，我不确定如何为 MATLAB 推断这些方法。

根据表 3，我需要使用 5 阶或 6 阶多项式来获得大约 23 位（小数点后）的精度。

我计划首先将所有传入的 theta 范围缩小到 -pi/4 到 +pi/4 之间，然后根据需要执行正弦或余弦函数（最终目标是实现 exp(i*x) = cos(x) + i*sin(x)。

我自己也许可以按照本文的说明进行操作，但我不知道如何将 remez 函数用于我的目的。另外，我不明白为什么作者使用方程（6）（第 9 页），我也不明白 k 的方程（第 11 页）是如何确定的（2796201 来自哪里？）以及为什么定义我们希望得到的多项式形式最终变为 sin9x) = x + kx^3 + x^5*P(x^2)。

使用 firpm 函数会更好吗（因为 remez 已弃用）？

谢谢您，我们非常感谢您提供所有帮助和指导，并进行编辑以确保为我的问题提供最佳答案。

Answer 1

我不会费心尝试开发您自己的近似值。更简单的是拿起 Hart 等人的“Computer Approximations”的副本。一个好的大学图书馆应该有它。 23 位大约是 7 个十进制数字，因此只需选择一个近似值，即可获得所需的准确度。您可以选择简单的多项式近似，也可以使用有理多项式，只要您能容忍除法，通常会更好一些。

范围缩小确实有意义，事实上，我在自己的工具中选择了相同的范围 (+/-pi/4)，因为这种范围选择特别容易使用。

编辑：（使用在 Hart 中可以找到的近似值的示例。）

在这里，我将找到 sin(x) 的近似值，其中 x 位于区间 [0,pi/4] 中。我的目标是在该区间内选择绝对准确度至少为 1.e-7 的近似值。当然，如果 x 为负值，我们知道 sin(x) 是一个奇函数，所以这是微不足道的。

我注意到 Hart 中的近似值倾向于采用 sin(alphapix) 的形式，其中 x 位于区间 [0,1] 中。如果我随后选择 alpha = 1/2 的近似值，我将得到一个在所选间隔内有效的近似值。因此，对于区间 [0,pi/4] 的近似值，我们寻找 alpha = 1/4。

接下来，我寻找一个绝对准确度至少为 7 位左右的近似值，我更喜欢使用有理多项式近似值，因为它们的效率更高一些。浏览第 118 页上的表格（我的 Hart 副本来自 1978 年），我发现 alpha = 1/4 的近似值符合要求：索引 3060。

这个近似的形式是

sin(alpha*pi*x) = x*P(x^2)/Q(x^2)

所以现在我切换到给出 SIN 3060 系数的页面。在我的副本中，它位于第 199-200 页。有5个系数，P00、P01、P02、Q00、Q01。 （注意这里使用的有点不标准的科学记数法。）因此，P（分子多项式）有 3 项，而分母 Q 有 2 项。写出来，我明白了：

sin(alpha*pi*x) = (52.81860134812 -4.644800481954*x^3 + 0.0867545069521*x^5)/ ...
    (67.250731777791 + x^2)

现在让我们在 MATLAB 中尝试一下。

x = linspace(0,pi/4,10001);
xt = x*4/pi; % transform to [0,1]
sine = @(x) (52.81860134812*x -4.644800481954*x.^3 + ...
     0.0867545069521*x.^5)./(67.250731777791 + x.^2);

max(abs(sin(x) -sine(xt)))
ans =
   1.6424e-09

plot(x,sin(x)- sine(xt),'-')

注意 y 轴上的 1e-9。

看起来这是在该特定区间内对 sin(x) 进行近似的最合理选择，尽管这提供了大约 29 位的准确度，而不是您要求的 23 位。如果您愿意选择不同的范围缩小区间，有一些选择可能会节省一个术语，可能会以您不需要的几位为代价。

log2(max(abs(sin(x) -sine(xt))))
ans =
      -29.182