Blogs

Chebyshev近似值以及如何帮助您省钱,赢得朋友和影响人们

杰森萨赫斯2012年9月30日20评论

好吧......也许这是一个伸展。我不认为我可以推荐任何东西来帮助你赢得朋友。不是我的强项。

但我将尝试说服你为什么要了解Chebyshev近似,这是一种用于计算如何尽可能接近计算数学函数的方法,具有最小的设计工作量和CPU力量。让我们探索两个用例:

  • Amy has a low-power 8-bit microcontroller and needs to compute \( \sqrt{x} \) (square root) to get an accurate estimate of RMS power dissipated in a current sense resistor. She's finding that the math library implementation for sqrt() takes a lot of CPU time, and she's short on free program memory space, and is wondering if there is another way to compute square roots that takes less time and program memory.
  • 比尔正在设计二氧化硫(所以2)探测器,具有输出依赖的电压的传感器2 气体浓度,以及将该电压转换为读取的微控制器2 浓度所以可以显示它。他对CPU时间限制并不大关,但不知何故,他必须弄清楚如何将气体传感器电压转换为浓度,因为传感器相当非线性。

在艾米的情况下,她有一个已知的功能似乎难以计算。在法案的情况下,他有一个未知的功能,他必须弄清楚如何模拟。这些都是函数评估的情况,具有两个非常不同的优先级。两者都在嵌入式系统上,这将我带来了嵌入式系统上数学计算的中心假设(Chomcoes?它需要一个更好的名字,我不想称之为Sachs假设):

本文以PDF格式提供,便于打印
需要精确浮点计算的唯一嵌入式系统是桌面计算器。

让那个沉入一分钟。你的电脑是一件好事。它可以以毫秒为单位计算PI至数百个小数位,或者在喷射发动机刀片周围模拟气流模式。计算机的起源一直以数学计算始终集中,在第二次世界大战结束时以及通过固态晶体管的早期发展巨大的发展。如果您在1940年代和1950年期间需要计算机,那就实现了一些精确的数学计算。拥有和运营的计算机非常昂贵;如果您需要快速估计一些数学,则存在幻灯片规则。这也是整个书籍被公布的含有工程功能表的时代,因此如果您需要将Bessel函数计算为5个小数位,可以通过在表格中查找它来完成它。

桌面计算机使用“全机精度”的计算策略:对于给定的位表示,SIN(x)只有一个正确的答案,并且它是一个比特模式,它是从精确的数学中产生最小的错误计算。有人已经写了并测试了这样做的数学软件图书馆,从那里出来,它只是一个函数调用。

另一方面,嵌入式系统往往包含廉价的处理器,这些处理器正在使用模数转换器测量进行。您可以为IEEE-754的完整53位进行数学,但如果您的ADC读数或其前面的信号调节,则只需准确到0.1%(10位),那么就会被数学错误赋予数学错误信号调理误差。

当您有一个能够容忍数学错误的系统时,并且具有计算资源限制的成本驱动,突然的数学计算变化的优先级。现在我们必须在想要计算某些东西时决定我们真正需要的东西。

让我们第一次取消艾米的问题: 我们需要计算平方根.

该陈述只是概述所需的内容。有两个主要决定需要做出前线:

  • 什么 范围 是否需要计算x(\ sqrt {x} \)?
  • 什么样的 准确性 是否需要\(\ sqrt {x})? (它如何随x而变化?)

一般全机精密数学图书馆对这些问题的简单答案:需要容忍任何输入,对于每个输入,也有一个正确的输出,即使正确的输出是错误(例如\(\ sqrt {x} \)对于不使用复数的系统中的负输入)。

在嵌入式系统中,准确性和范围要求是判断呼叫。

准确性是一件容易理解的事情:有算法占据计算表达式的快捷方式,并且可以容忍更加不准确,可以拍摄更多的快捷方式。

范围要求有点微妙。显然,您必须支持对算法的预期输入范围,也许您想要留下一些工程边缘。但是为什么计算\(\ sqrt {x})在更大的范围内更昂贵?如果我要告诉你,我更关心计算的范围比准确性,你会相信我吗?

让我们使用一个真实的例子。艾米认为她需要\(\ sqrt {u})计算在范围\(u \中[0.2,5.0] \中)。她的第一个本能就是去寻找她的微积分教科书并查找泰勒系列:

如何不评估功能,第1部分:泰勒系列

这是泰勒系列for \(\ sqrt {1 + x}):   (方程式提供 维基百科)

但是有一点脚注说这只适用于-1和1之间的x,这意味着它不会用于计算\(\ sqrt {2} \)或任何2或更大的方形根源。对于\(\ sqrt {5} \),她必须使用一个技巧,即\(\ sqrt {5} \)= \(2 \ times \ sqrt {5/4} \)。因此,她决定限制泰勒系列的实际输入范围,以覆盖\(\ sqrt {1 + x} \)的范围\(x \中[-0.8,+ 0.25] \),如果\(u = 1 + x>1.25 \)那么艾米将做除了第4个技巧:(我们会在这里使用Python;假设amy在c)中相同的东西

def sqrt1(x):
# returns sqrt(1+x) for x between -0.8 and +0.25
   return 1 + x/2 - x*x/8 + x*x*x/16
def sqrtapprox(u):
# valid for x between 0.2 to 5.0
   if (u >= 1.25):
      return sqrt1(u/4-1)*2
   else:
      return sqrt1(u-1)

这足够好吗?嗯,让我们这么想。我们希望了解函数与其近似之间的错误。让我们的图表错误。这是一个快速的Python函数,它使用numpy和matplotlib绘制结果:

import numpy as np
import matplotlib.pyplot as plt
def showplots(f,approxlist,a,b):
    x = np.linspace(a,b,1000)
    plt.figure(1)
    plt.subplot(211)
    plt.plot(x,f(x))
    vfuncs = [np.vectorize(approx) for approx in approxlist]
    for vf in vfuncs:
        plt.plot(x,vf(x))
    plt.xlim(a,b)
    plt.ylabel('f(x) and approximations fa(x)')
    plt.subplot(212)
    for vf in vfuncs:
        plt.plot(x,f(x)-vf(x))
    plt.xlim(a,b)
    plt.ylabel('error = f(x)-fa(x)')
    plt.xlabel('x')
    plt.show()

The results of showplots(np.sqrt,[sqrtapprox],0.2,5.0) will show us the error between the full-machine-precision sqrt() 和 our approximation:

哎呀。最坏情况的误差为0.038。让我们在Taylor系列中尝试更多条款。这意味着我们必须弄清楚系数如何工作。一点铅笔和纸代数显示从二次术语开始的系数是 $$\begin{eqnarray} a_2&=& - \ frac {1} {4 \ times 2} = - \ frac {1} {8},\ cr a_3&=&+ \ frac {3 \ times 1} {6 \ times 4 \ times 2} = \ frac {1} {16},\ cr a_4&=&\& - \ \ frac {5 \ times 3 \ times 1} {8 \ times 6 \ times 4 \ times 2} = \ frac {5} {128},\ cr a_5&=&+ \ frac {7 \ times 5 \ times 3 \ times 1} {10 \ times 8 \ times 6 \ times 4 \ times 2} = \ frac {7} {256},\ cr A_6&=& - \& - \ frac {9 \ times 7 \ times 5 \ times 3 \ times 1} {12 \ times 10 \ times 8 \ times 6 \ times 4 \ times 2} = - \ frac {21} {1024} \end{eqnarray}$$ 等等,所以我们再试一次:

def sqrt1a(x):
# returns sqrt(1+x) for x between -0.8 and +0.25
   return 1 + x/2 - x*x/8 + x*x*x/16 - 5.0/128*x*x*x*x
def sqrtapproxa(u):
# valid for x between 0.2 to 5.0
   if (u >= 1.25):
      return sqrt1a(u/4-1)*2
   else:
      return sqrt1a(u-1)
def sqrt1b(x):
# returns sqrt(1+x) for x between -0.8 and +0.25
   return 1 + x/2 - x*x/8 + x*x*x/16 - 5.0/128*x*x*x*x + 7.0/256*x*x*x*x*x
def sqrtapproxb(u):
# valid for x between 0.2 to 5.0
   if (u >= 1.25):
      return sqrt1b(u/4-1)*2
   else:
      return sqrt1b(u-1)

以下是Sqrtappoxa(第4位Taylor系列)的结果:

以下是Sqrtappoxb的结果(第5位Taylor系列):

哎呀。我们添加了两种术语,我们只会将错误降至0.015。

以下是您需要了解泰勒系列进行功能评估。 (或者相反,阅读这一点,无论你不理解什么,都相信我。)

为了在点\(x = x_0 \)周围的区域中的函数f的函数f的系数,评估\(f(x_0)\)和\(f(x)\)的所有衍生物在\(x = x_0 \)。您不必在别的地方查看功能;以\(x = x_0 \)的行为以某种方式让您知道\(f(x)\)是别的地方。对于某些功能(如\(\ sqrt {1 + x}))收敛区域是有限的;对于\(\ SIN X \)或\(e ^ x \)收敛区域到处都是。

泰勒系列在大多数日常功能上工作,因为这些功能有一个名为“Analyticity”的财产:它们的所有衍生品都存在,因此是连续和平滑的,并且通过数学的魔力这意味着只在1点看分析功能告诉您足够的信息来在附近的一些地区计算它。

右边\(x = x_0 \),泰勒系列很棒。相应泰勒序列的确切功能与前N​​项之间的误差将非常小,左右非常小,x = x_0 \)。当您远离\(x = x_0 \)时,错误开始变大。通常,您使用的术语越多,近似越突出\(x = x_0 \),但是当近似在其范围边缘开始做得不好时,它会更快地分流。

总结一下: 泰勒级近似不公平分布近似误差。 误差在\(x = x_0 \)处小,但其误差在其范围的边缘处趋于更大。

现在,我将向您展示\(\ sqrt {x} \)\(x \中的[0.2,5.0] \)中的两个尝试,但在解释我的想出之前,我将等待他们。

def sqrt1c(x):
# valid for x between 0.2 to 1.25
   coeffs = [0.20569678, -0.94612133, 1.79170646, -1.89014568, 1.66083189, 0.17814197]
   y = 0
   for c in coeffs:
      y = y * x + c
   return y
def sqrtapproxc(u):
# valid for u between 0.2 to 5.0
   if (u >= 1.25):
      return sqrt1c(u/4)*2
   else:
      return sqrt1c(u)
def sqrtapproxd(x):
# valid for x between 0.2 to 5.0
   coeffs = [0.00117581,-0.01915684,0.12329254,-0.41444219,1.04368339, 0.26700714]
   y = 0
   for c in coeffs:
      y = y * x + c
   return y

近似“C”拆分\([0.2,1.25)\)之间的范围,并如前所述\([1.25,5] \);近似“d”在整个范围内操作([0.2,5] \)。两者都是5度多项式。这是近似“C”的性能:

您会发现一个约0.00042的峰值误差。这比近似“B”的错误小35倍(使用具有相同范围分割的第5阶多项式的Taylor系列版本)。

这是近似值“d”的错误:

这是约0.011的最坏情况误差,比近似值“b”小约30%。我们得到了更好的最坏情况错误,而无需分开范围,具有相同程度的多项式!

从这项练习中要注意的几件事:

  • 良好的函数近似具有具有特征“Wiggle”的错误图。 理想情况下,它们形成了所谓的最小函数,其中正负峰值所有排队。这在整个输入参数范围内传播了错误,并且可以是比同一学位的泰勒串多项式更好的级别的级。在实践中,难以实现理想的最小误差图,但是有一个简单的方法靠近它,我将在一下解释。
  • 范围减少可以显着提高功能评估的准确性。 近似“C”和“D”两者使用5度多项式,但近似“D”使用整个输入范围(25:1比率为最大输入到最小输入),而近似“C”仅使用了输入范围的一部分(最大值为6.25:1比率为最小输入),“C”的比例更高25倍。

如何不评估函数,第2部分:查找表

好吧,让我们扔掉泰勒系列。下一个方法是查找表。这个想法很简单:你有一个针对你的函数的预定值数组,并使用那些作为表来评估函数的表。让我们在\([0.2,5.0] \)范围内使用带有32个条目的查找表来查看\(\ sqrt {x} \)的解决方案:

import math
def maketable0(func, a, b, n):
    table = func([(x+0.5)/n*(b-a)+a for x in range(0,n)])
    def tablefunc(x):
        if x < a or x > b:
            raise "out of range error"
        i = math.floor((x-a)/(b-a)*n)
        if i == n:
            i = n-1
        return table[i]
    return tablefunc
sqrtapproxe = maketable0(np.sqrt, 0.2, 5.0, 32

使用查找表涉及以下过程,这是相当微不足道的,但值得讨论:

  • 在设计时间(或初始化时间),将函数的输入范围划分为n个间隔并评估这些间隔内的功能。 (通常间隔是等于宽度,并且功能评估位于每个间隔的中点。)
  • 在运行时:
    • 处理超出范围的情况(抛出错误或返回对应于MIN / MAX输入的值)
    • 计算函数中的间隔的索引
    • 查找该索引的表中的值

以下是近似值“e”的错误看起来像:

请注意楼梯级结果和锯齿错误绘图。这对于普通表的方法是典型的。最大误差大致与间隔内函数的斜率大致成比例(对于x的小值,误差较大,其中\(\ sqrt {x} \)较陡峭),并且与间隔大小成比例。如果我们想将误差减少10倍,我们需要10倍的查找表条目。

如果我们想要改善情况,我们可以使用的下一个改进是线性插值:

  • 在设计时间(或初始化时间),将函数的输入范围划分为n个间隔,并评估这些间隔的N + 1端点处的功能。 (通常间隔是相等的宽度)
  • 在运行时:
    • 处理超出范围的情况(抛出错误或返回对应于MIN / MAX输入的值)
    • 计算函数中的index k
    • 查找该索引和下一个索引的表中的值
    • perform a linear interpolation: y = table[k]*(1-u) + table[k+1]*u, where u is the number between 0 and 1 representing the position of the input within the interval (0 = beginning of interval, 1 = end)
import math
def maketable1(func, a, b, n):
    table = func([float(x)/n*(b-a)+a for x in range(0,n+1)])
    def tablefunc(x):
        if x < a or x > b:
            raise "out of range error"
        ii = (x-a)/(b-a)*n
        i = math.floor(ii)
        if i == n:
            i = n-1
        u = ii-i
        return table[i]*(1-u) + table[i+1]*u
    return tablefunc
sqrtapproxf = maketable1(np.sqrt, 0.2, 5.0, 32)

这是近似“f”的错误图:

更好 - 请注意包含抛物线部分的错误曲线。这对于插值表查找是典型的。最大误差大致与间隔内函数的曲率(第2衍生物)成比例,并且与间隔大小的平方成比例。如果我们想将误差减少10倍,我们只需要3倍的查找表条目。

这是表插值可能比多项式评估更合适的罕见情况之一:这里的最大误差是大约1/2,如近似“D”所用的5阶多项式。

有关执行时间的比较:表明没有插值的表明大致相当于二次(第二阶)函数,表插值大致相当于四分之一(第4阶)多项式。这只是一个估计;它真的取决于你是否正在进行浮点或定点计算,是表格是否具有2的电源的长度(它应该是,因此您可以使用按位操作来加快数学),在处理器上执行表查找需要多长时间。

在我看来,只有几种情况,表查找是合适的:

  • 如果精度要求非常低,则可能就是不插入的表查找。
  • 如果函数非常“nonpolynomial”,那么与可比执行时间的多项式评估相比,内插表查找可能会产生卓越的准确性。 (这似乎是\(\ sqrt {x} \)\([0.2,5.0] \)范围内的情况。我们很快就定义了“nonpolynomial”;因为现在只是将它视为一个函数很多Wiggles或曲折。)
  • 如果功能非常难以在整个范围内与多项式近似,但需要非常高的精度,则一个选项是将该范围分成4或8个子内部,并使用表查找以获得每个子宫的多项式系数。

表查找是 不是 适用于您需要如此大的桌子大小,您必须担心内存要求。

另请注意,近似值“f”的错误不相当分布:在输入范围的低端处发生最大误差,在其余的输入范围内更小。您可以使用技巧来使错误分布更公平,即使用变量宽度间隔,但在运行时仍然使用有效计算时没有一般过程。

Chebyshev多项式

使用多项式来评估函数的关键是不是想到由\(1,x,x ^ 2,x ^ 3 \)等的线性组合组成的多项式,而是作为线性组合 Chebyshev多项式 \(t_n(x)\)。正如Wikipedia页面中更详细讨论的那样,这些具有以下特殊属性:

  • 它们在范围内是最小的函数\([-1,1]),例如所有最小值/最大/极值/极值位于±1(见下图)
  • 它们在具有权重的范围\([-1,1] \)上正交\(\ frac {1} {\ sqrt {1-x ^ 2}})
  • \(t_0(x)= 1,t_1(x)= x,t_ {n + 1}(x)= 2xt_n(x) - t_ {n-1}(x)\)
  • \(t_n(x)= \ cos(n \ cos ^ { - 1} x)\)

图:chebyshev函数\(t_n(x)\),由维基百科提供

假设我们在Range \(x \中有一个给定的函数\(f(x)\))。我们可以表达\(f(x)= \ sum c_k t_k(u)\)其中\(u = \ frac {2x-ab} {ba} \)和\(x = \ frac {ba} {2} u + \ FRAC {A + B} {2} \),它是将间隔\(x \中的[a,b] \中的x \)映射到\(u \在[-1,1] \中)的转换。然后,问题变得有问题,用于找出可以使用的系数\(C_K \),这可以使用该系数\(c_k \) Chebyshev节点:

这些更容易理解图形;它们仅仅是半圆形的相等间隔圆弧的X轴上的投影:

您将注意到,节点间隔在0附近0,更加压缩至±1。

为了通过第一n Chebyshev多项式(k = 0至n-1)的线性组合来近似函数,系数\(c_k \)简单地等于产品的平均值的\(a(k)\)倍(t_k(u)f(x)\)在n Chebyshev节点上评估,其中\(k = 0 \)和所有其他k的\(a = 2 \)的\(a = 1 \)。

让我们用一个例子更清楚地说明这个过程。

假设\(f(x)= \ frac {1} {3} x ^ 3 + 2x ^ 2 + x - 10 \)在范围内\([-1,3] \),with \(u = \ frac {x-1} {2},x = 2u + 1 \):

让我们使用\(n = 5 \)。 N = 5的Chebyshev节点是\(u = -0.951057,-0.58785,0,+0.587785,+0.951057 \)。这对应于\(x = -0.902113,-0.175571,1,2.175771,2.902113 \),并且在这些节点评估的功能\(f(x)\)是\(y = -9.51921,0.10.11572,-6.66667, 5.07419,17.89408 \)。

对于\(c_0 \),我们计算Y = \的平均值((-9.51921-10.11572-6.6667 + 5.07419 + 17.89408)/ 5 = -0.66667 \)。

对于\(c_1 \),我们计算2×u×y = \的平均值(2 \ times(-9.51921 \ times-0.951057 + -10.11572 \ times-0.587785 \)\(+ - 6.66667 \ times 0 + 5.07419 \时代0.587785 + 17.89408 \ times 0.951057)= 14 \)

for \(c_2 \)我们计算2×平均值\(t_2(u)\ times y =(2u ^ 2-1)\ times y \:\ longrightarrow \:c_2 = 6 \)

对于\(c_3 \),我们计算2×平均值\(t_3(u)\ times y =(4u ^ 3-3u)\ times y \:\ longrightarrow \:c_3 = 0.66667 \)

for \(c_4 \)我们计算2×平均值\(t_4(u)\ times y =(8u ^ 4-8u ^ 2 + 1)\ times y \:\ longrightarrow \:c_4 = 0 \)

这里结论是\(f(x)= - \ frac {2} {3} t_0(u)+ 14t_1(u)+ 6t_2(u)+ \ frac {2} {3} t_3(u)\) ;如果您通过代数进行咬嚼,您将看到结果与\(f(x)= \ frac {1} {3} x ^ 3 + 2x ^ 2 + x-10 \)相同。

那么这里的重点是什么......?好吧,如果你 已经 具有多项式的系数,在使用Chebyshev多项式中没有任何点。 (如果输入范围与\([1,1])有显着不同,则从x到u的线性变换,其中\(u \在[-1,1]中)将为您提供更好的多项式数值稳定性,Chebyshev多项式是做这种转变的一种方法。更稍后的是,如果你想用较小程度的多项式近似多项式,Chebyshev多项式将给你最好的近似值。

例如,让我们使用二次函数来在有问题的输入范围内近似\(f(x)\)。我们需要做的就是截断Chebyshev系数:让我们计算$$ \ begin {eqnarray} f_2(x)&& - \ frac {2} {3} t_0(u)+ 14t_1(u)+ 6t_2(u) \ Cr. &=& - \ frac {2} {3} + 14u + 6 \ times(2u ^ 2-1)\ cr &=& - \ frac {20} {3} + 14u + 12u ^ 2 \ cr &=& - \ frac {20} {3} + 7(x-1)+ 3(x-​​1)^ 2 \ cr &=&3x ^ 2 + x - \ frac {32} {3}。 \end{eqnarray} $$

瞧!我们具有靠近原始立方体功能的二次函数,并且错误的幅度只是Chebyshev多项式的系数,我们删除了= \(\ FRAC {2} {3} {3})。

您将注意到这两个多项式(\(\ frac {1} {3} x ^ 3 + 2x ^ 2 + x - 10 \)和\(3x ^ 2 + x - \ frac {32} {3} \) )以正常的功率系列形式表达,彼此没有明显的关系,而Chebyshev系数除外相同,除了系数\(c_3 \),我们设置为零以获得二次函数。

在我的观点中,这是在其Chebyshev系数方面表达函数的主要实用程序:系数\(C_0 \)告诉您在输入范围内的函数的平均值; \(c_1 \)告诉您函数的“线性”,\(c_2 \)告诉您函数的“二次性”,\(c_3 \)告诉您函数的“立方体”等。

这是一个更典型和更不凡的例子。假设我们希望在Range \([1,2] \)上近似\(\ log_2 x \)。

Chebyshev系数的分解高达6度产生以下内容:

$$\begin{eqnarray} c_0 &=& +0.54311 \cr c_1 &=& +0.49505\cr c_2 &=& -0.042469\cr C_3&=&+0.0048577 \ cr C_4&&&-6.2508 \ times 10 ^ { - 4} \ cr C_5&=&+8.5757 \ times 10 ^ { - 5} \ cr c_6&&&-1.1996 \ times 10 ^ { - 5} \ cr \end{eqnarray} $$



您将注意到,每个系数仅是先前系数的大小的一小部分。这是一个指标,其中所讨论的函数(在指定的输入范围内)对多项式近似值非常好。如果Chebyshev系数不会明显减小第五系数,则我将在感兴趣范围内分类为“非突变”。

\(\ log_2 x \)的6度Chebyshev近似值的最大误差仅为[1,2] \中)仅为2.2×10-6,这大约是下一个Chebyshev系数的值。

如果我们将系列截断到4度Chebyshev近似值,我们会获得大约等于\(C_5 = 8.6×10 ^ { - 5} \)的最大错误:

整洁,呵呵?

基本功能的一些进一步的例子

以下是一个常用函数的列表,以及前6个Chebyshev系数:

\(f(x)\) 范围 c0 c1 c2 c3 c4 c5
\(\ sin \ pi x \) [-0.5,0.5] 0 1.1336 0 -0.13807 0 0.0045584
\(\ sin \ pi x \) [-0.25,0.25] 0 0.72638 0 -0.01942 0 0.00015225
\(\ cos \ pi x \) [-0.5,0.5] 0.472 0 -0.4994 0 0.027985 0
\(\ cos \ pi x \) [-0.25,0.25] 0.85163 0 -0.14644 0 0.0019214 0
\(\ sqrt {x} \) [1,4] 1.542 0.49296 -0.040488 0.0066968 -0.0013836 0.00030211
\(\ log_2 x \) [1,2] 0.54311 0.49505 -0.042469 0.0048576 -0.00062481 8.3994E-05.
\(e ^ x \) [0,1] 1.7534 0.85039 0.10521 0.0087221 0.00054344 2.7075E-05.
\(\ frac {2} {\ pi} \ tan ^ { - 1} x \) [-1,1] 0 0.5274 0 -0.030213 0 0.0034855
\(\ frac {1} {1 + e ^ { - x}} \) [-1,1] 0.5 0.23557 0 -0.0046202 0 0.00011249
\(\ frac {1} {1 + e ^ { - x}} \) [-3,3] 0.5 0.50547 0 -0.061348 0 0.01109
\(\ frac {1} {1 + x ^ 2} \) [-1,1] 0.70707 0 -0.24242 0 0.040404 0
\(\ frac {1} {1 + x ^ 2} \) [-3,3] 0.30404 0 -0.29876 0 0.12222 0

其他的建议

Chebyshev功能的很大参考是 数值食谱 通过印刷机,Teukolsky,Vetterling和Flannery,其覆盖Chebyshev近似详细。

评估Chebyshev功能时有几件事可以注意到:

  • 最好直接计算函数而不是尝试将Chebyshev近似转换为标准多项式形式。 (即,如果要计算,如果要计算\(3t_0(x)+ 5t_1(x) - 0.1t_2(x)\),请勿转换为多项式,而是直接计算T值。)数字配方推荐克伦斯彻底的复发公式,但我只是使用以下简单算法:
u = (2*x-a-b)/(b-a)
Tprev = 1
T = u
y = c[0]
for i = 1:n
    y=y+T*c[i]
   Tnext = 2*u*T - Tprev
   Tprev = T
   T = Tnext
  • 请记住从X转换为U - 这会保持算法的数值稳定性,因为多项式倾向于在评估-1和1之间的参数的良好调节,并且对于具有非常大或非常小的争论的参数的病变。

以下是用于计算和应用Chebyshev系数的Python类:

import math
import numpy as np
def chebspace(npts):
    t = (np.array(range(0,npts)) + 0.5) / npts
    return -np.cos(t*math.pi)
def chebmat(u, N):
    T = np.column_stack((np.ones(len(u)), u))
    for n in range(2,N+1):
        Tnext = 2*u*T[:,n-1] - T[:,n-2]
        T = np.column_stack((T,Tnext))
    return T
class Cheby(object):
    def __init__(self, a, b, *coeffs):
        self.c = (a+b)/2.0
        self.m = (b-a)/2.0
        self.coeffs = np.array(coeffs, ndmin=1)
    def rangestart(self):
        return self.c-self.m
    def rangeend(self):
        return self.c+self.m
    def range(self):
        return (self.rangestart(), self.rangeend())
    def degree(self):
        return len(self.coeffs)-1
    def truncate(self, n):
        return Cheby(self.rangestart(), self.rangeend(), *self.coeffs[0:n+1])
    def asTaylor(self, x0=0, m0=1.0):
        n = self.degree()+1
        Tprev = np.zeros(n)
        T     = np.zeros(n)
        Tprev[0] = 1
        T[1]  = 1
        # evaluate y = Chebyshev functions as polynomials in u
        y     = self.coeffs[0] * Tprev
        for co in self.coeffs[1:]:
            y = y + T*co
            xT = np.roll(T, 1)
            xT[0] = 0
            Tnext = 2*xT - Tprev
            Tprev = T
            T = Tnext
        # now evaluate y2 = polynomials in x
        P     = np.zeros(n)
        y2    = np.zeros(n)
        P[0]  = 1
        k0 = -self.c/self.m
        k1 = 1.0/self.m
        k0 = k0 + k1*x0
        k1 = k1/m0
        for yi in y:
            y2 = y2 + P*yi
            Pnext = np.roll(P, 1)*k1
            Pnext[0] = 0
            P = Pnext + k0*P
        return y2
    def __call__(self, x):
        xa = np.array(x, copy=False, ndmin=1)
        u = np.array((xa-self.c)/self.m)
        Tprev = np.ones(len(u))
        y = self.coeffs[0] * Tprev
        if self.degree() > 0:
            y = y + u*self.coeffs[1]
            T = u
        for n in range(2,self.degree()+1):
            Tnext = 2*u*T - Tprev
            Tprev = T
            T = Tnext
            y = y + T*self.coeffs[n]
        return y
    def __repr__(self):
        return "Cheby%s" % (self.range()+tuple(c for c in self.coeffs)).__repr__()
    @staticmethod
    def fit(func, a, b, degree):
        N = degree+1
        u = chebspace(N)
        x = (u*(b-a) + (b+a))/2.0
        y = func(x)
        T = chebmat(u, N=degree)
        c = 2.0/N * np.dot(y,T)
        c[0] = c[0]/2
        return Cheby(a,b,*c)

您可以将其放在文件Chebby.py中,然后用5度多项式在0和\(\ frac {\ pi} {2} \)之间的\(\ sin x \)之间使用它在角度0,\(\ frac {\ pi} {6} \),\(\ frac {\ pi} {4} \),\(\ frac {\ pi} {3} {3} \)) :

>>> import cheby
>>> import numpy as np
>>> import math
>>> c = cheby.Cheby.fit(np.sin,0,math.pi/2,5)
>>> c
Cheby(0.0, 1.5707963267948966, 0.60219470125550711, 0.51362516668030367, -0.10354634422944738, -0.013732035086651754, 0.001358650338492214, 0.00010765948465629727)
>>> c(math.pi*np.array([0, 1.0/6, 1.0/4, 1.0/3]))
array([  6.21628624e-06,   5.00003074e-01,   7.07099696e-01,
         8.66028717e-01])

经验数据的功能近似

让我们回到我们对账单问题的讨论开始。他具有与电压测量的物理量(气体浓度)相关的非线性传感器;他面临的问题是如何反转该关系,以便他可以估计原始物理量:\(y = f(x)\),其中x是测量,y是物理量。

Chebyshev近似在这里工作好,但我们有问题。我们不能只评估Chebyshev节点的\(f(x)\)以获取Chebyshev系数,因为我们不知道某些\(f(x)\)实际上是什么。相反,我们必须进行一些参考测量来确定那些具有已知物理量的系数(例如,单个情况下已知的气体浓度),并制表这些X和Y值,然后使用它们来估计系数。我们在选择测量时有两种选择:

  • 尝试规划测量,使得值X在Chebyshev节点附近通过预期范围。这将使我们能够使用较少的测量。
  • 制作大量测量,因此整个输入范围都很好

在某些情况下,它可以快速且易于进行任意参考。 (例如,如果您尝试校准麦克风,则可以使用电子设备自动更改参考声音信号的幅度和频率),第二个选项是可行的。在其他情况下,与气体浓度问题一样,制备大量测量可能是不方便或昂贵的,并且您必须仔细选择参考值。

我发现从经验数据估算Chebyshev系数的最佳方式是使用 加权最小二乘:

  • 将测量x和物理输入量y收集到向量x和y中
  • 对于测量X的向量,计算u的x,并计算Chebyshev PolyNomials \(t_0(u),t_1(u),t_2(u),\ ldots t_n(u)\)的值并将其保存为矩阵T的载体。每个多项式将用作基本函数。
  • 计算加权值\(w(x)\)以补偿值x的浓度。 (例如,如果大量x值近\(x = 0 \)但很少有x(x = 1 \),您可能必须选择\(w(x)\)近\(x = 0 \)近\(x = 1 \)。)将这些值表达为对角线矩阵W.
  • 使用加权最小二乘来计算C系数\(C_0,C_1,C_2,\ LDOTS C_N \)的矢量C来解决C:
    $$(t ^ twt)c = t ^ twy。$$ 为此,您需要某种矩阵库。在Matlab中,这将被计算为 C = (T'*diag(W)*T)\(T'*diag(W)*Y)。在python用numpy,你可以跳过计算t矩阵,只需使用 numpy.polynomial.chebyshev.chebfit 在X和Y测量上的功能;加权矢量是可选输入。

我不会使用大于5的程度。如果您发现您仍然没有与第5度多项式获得良好的近似,则可能使用错误的方法;该函数太尖锐或摇摆或尖叫。更高的多项式可能具有数值准确性的问题。

在任何情况下,经验数据的函数近似值比已知功能的近似略微棘手。 总是 确保通过绘制原始数据以及您提出的功能来仔细检查近似值。

包起来

下次需要转向函数近似时,请给乳酪近似拍摄!它不仅可能是近似多项式的最佳和最简单的方法,但它还将让您知道多项式近似于Chebyshev系数的行为。

确保明智地选择功能所需的准确性和输入范围。输入范围太高,输入范围太高,需要更高的多项式和更长的执行时间,并且高度多项式通常是数值不均衡的。

更有效的函数评估将意味着您可以使用较少的CPU资源以获得所需的准确性。所以这 能够 节省你的钱并影响你的同事。

快乐近似!


©2012 Jason M. Sachs,保留所有权利。


[]
评论 markrages.2012年9月30日
好东西。这些matplotlib剧会提醒我约翰亨特。 http://numfocus.org/johnhunter/
[]
评论 rickyroko2013年7月24日
伟大的帖子!这真的很方便,谢谢你的比较。我将制作一个小的变化是清楚地将复发关系中的生成功能定义为t_n(x):= cos(n cos-1(x))而不是cos(cos-1nx)。
[]
评论 JMS_NH.2015年8月25日
固定的!
[]
评论 JMS_NH.2015年7月26日
yikes - 这是一个错误,而不仅仅是一个建议。谢谢,我有机会修复。 (+将本文中的方程升级到Mathjax。)
[]
评论 jkvasan.2012年10月9日
伟大的帖子!你赢了我作为朋友,你的帖子肯定会影响我。
[]
评论 尼米什特2015年8月23日
Matlab中的Chebfun包是一个非常彻底的扩张这个想法
[]
评论 Macsdev.2015年11月16日
什么是一篇优秀的文章!谢谢。可能是你可以编写一个新的文章,详细介绍了帖子的这一部分:)
“我发现从经验数据估算Chebyshev系数的最佳方式是使用加权最小二乘”。
[]
评论 默博2013年7月28日
Chebyshev大约与atan的cordic(a / b):)
[]
评论 Davidf12.2015年7月21日
您的语法突出显示已破坏
[]
评论 lemzwerg.2015年8月25日
在句子后的Python代码中

以下是用于计算和应用Chebyshev系数的Python类:

第一个评论行后,一切都是灰色的
[]
评论 Calkuta.2015年8月27日
我遇到同样的问题,看起来在评论行之后都是灰色的任何东西。酷文!
[]
评论 JMS_NH.2015年7月21日
请添加更多详细信息 - 似乎对我来说很好。
[]
评论 jbwhitmore.2015年9月21日
语法突出显示也为我打破了(使用Chrome)。这是我所看到的屏幕截图: //www.dropbox.com/s/qtgnyvimiqicvvt/Screenshot%202015-09-22%2011.15.09.png?dl = 0

[]
评论 adam892015年8月23日
在您的定义tn(x)中的略带拼写。你写了:

tn(x)= cos(cos-1 nx)

当然,这简化了tn(x)= nx。我相信你的意思是:

tn(x)= cos(n * cos-1 x)
[]
评论 JMS_NH.2015年8月25日
固定的!
[]
评论 JMS_NH.2015年8月23日
YEP,Rickyroko注意到它也+它在我的Todo列表中以及转换为Mathjax。谢谢!
[]
评论 Samir12.2016年4月23日
感谢提供这篇好文章。当你计算C(0)对于AMY的问题时,我认为它应该是2 *平均值。不是吗?提前致谢。
[]
评论 Sopie62.2017年1月7日

我以为我明白了这篇文章 - 但我错过了一些东西。

对我来说,酸测试是试图了解地球上你有这些价值 sqrt1c(x)sqrtappoxd(x)。我尝试计算Chebyshev多项式,使用Excel电子表格进行逐步逐步,并与您的价值观出来,这些值看起来像你的东西。

我错过了什么?它是否与如何使用简单评估Chebyshev多项式进行评估。 

y = y * x + c

......当您的代码进一步下页面时更复杂?

[]
评论 atomq.2019年1月4日

非常感谢,伟大的东西Alltogratr。潜入故事后,我对W(x)对角线矩阵的系数有问题 加权最小二乘 方法。是关于如何将它们算出的过程是什么程序?我的猜测是系数的值与测量点值有关,它们的平均值与它们跨越范围或极端范围的密度?提前致谢!

[]
评论 JMS_NH.2019年9月29日

抱歉迟到的回复......我不是100%确定是否有一种来计算w(x)的系统。低密度=重量高,密度高=低重量。这可能意味着您希望将权重的累积和累积和k的线性函数是k的线性函数。 Blah Blah Blah,练习读者等。

我在10-12年前在以前的职业生涯中做了这样的事情,但我不记得细节。祝你好运!

要发布回复评论,请单击连接到每个注释的“回复”按钮。发布新的评论(不是回复评论),请在评论的顶部查看“写评论”选项卡。

注册将允许您参加所有相关网站的论坛,并为您提供所有PDF下载。

注册

我同意 使用条款隐私政策.

尝试我们偶尔但流行的时事通讯。非常容易取消订阅。
或登录