用MATLAB找到三角函数的Minimax多项式逼近

Question

我试图在MATLAB中使用remez交换算法找到正弦和余弦的minimax多项式近似.因为我正在实现IEEE-754浮点的正弦和余弦函数,所以需要精度到23位.

在这里使用此链接(参见第8页到第15页),给出了使用Mathematica和Maple查找多项式的指令,但是,我不确定如何为MATLAB推断这些方法.

根据表3,我需要使用5阶或6阶多项式来获得~23位(小数点后)的精度.

我计划首先将所有输入θ的范围缩小到-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等人.一个好的大学图书馆应该有它.23位大约是7位十进制数字,因此只需选择一个近似值即可获得所需的精度.您可以选择简单的多项式近似,或使用有理多项式,只要您能够容忍除法,通常会更好一些.

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

编辑:(使用Hart中可以找到的近似值的一个例子.)

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

我注意到Hart中的近似值倾向于sin(alpha pi x),其中x位于区间[0,1]中.如果我然后选择alpha = 1/2的近似值,我会得到一个在所选区间内有效的近似值.因此,对于区间[0,pi/4]的近似值,我们寻找alpha = 1/4.

接下来,我寻找一个近似值,表示绝对精度至少为7位左右,我更倾向于使用有理多项式近似,因为它们往往更有效率.扫描表格(第118页)(哈特的副本是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),'-')

注意1e-9连接到y轴.

看起来这是在特定区间内对sin(x)进行近似的最合理选择,尽管这给出了大约29位的精度而不是你要求的23位.如果您愿意选择不同的范围缩小间隔,则可以选择一些可以节省一个术语的选项,可能需要花费您不需要的几个位.

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