c - 使用 GNU Scientific C 的多项式拟合

Question

我需要从 n+1 个数据点得到一个 n 次函数。通过使用以下 gnuplot 脚本，我得到了正确的拟合:

f(x) = a + b*x + c*x**2 + d*x**3 + e*x**4 + f*x**5 + g*x**6 + h*x**7 + i*x**8 + j*x**9 + k*x**10 + l*x**11 + m*x**12

# Initial values for parameters
a = 0.1
b = 0.1
c = 0.1
d = 0.1
e = 0.1
f = 0.1
g = 0.1
h = 0.1
i = 0.1
j = 0.1
k = 0.1
l = 0.1
m = 0.1


# Fit f to the following data by modifying the variables a, b, c
fit f(x) '-' via a, b, c, d, e, f, g, h, i, j, k, l, m
   4.877263 45.036000
   4.794907 44.421000
   4.703827 43.808000
   4.618065 43.251000
   4.530520 42.634000
   4.443111 42.070000
   4.357077 41.485000
   4.274298 40.913000
   4.188404 40.335000
   4.109381 39.795000
   4.027594 39.201000
   3.946413 38.650000
   3.874360 38.085000
e

拟合后得到如下系数:

a               = -781956        
b               = -2.52463e+06   
c               = 2.75682e+06    
d               = -553791        
e               = 693880         
f               = -1.51285e+06   
g               = 1.21157e+06    
h               = -522243        
i               = 138121         
j               = -23268.8       
k               = 2450.79        
l               = -147.834       
m               = 3.91268 

然后，通过将数据和 f(x) 一起绘制，似乎给定的系数是正确的:

但是，我需要使用 C 代码来获得这样的拟合。在某些情况下，多项式拟合 ( as in this link ) 的 GNU 科学图书馆代码结果正确。但是对于上述数据(以及我数据集中的其他几个案例)，我得到的结果是有缺陷的。

例如，以下代码(使用与上述示例相同的数据):

void testOfPolynomialFit(){
   double x[13] = {4.877263, 4.794907, 4.703827, 4.618065, 4.530520, 4.443111, 4.357077, 4.274298, 4.188404, 4.109381, 4.027594, 3.946413, 3.874360};
   double y[13] = {45.036000, 44.421000, 43.808000, 43.251000, 42.634000, 42.070000, 41.485000, 40.913000, 40.335000, 39.795000, 39.201000, 38.650000, 38.085000};
   double coefficients[13];

   polynomialfit(13, 13, x, y, coefficients);

   int i, n = 13;

   for (i = 0; i < n; i++)
   {
    printf("%lf\t", coefficients[i]); 
   }
   printf("\n"); 

}

结果:

-6817581083.803348      12796304366.105989      -9942834843.404181      3892080279.353104             
 -630964566.517794        -75914607.005088         49505072.518952        -5062100.000931 
   -1426228.491628           514259.312320           -70903.844354            4852.824607
       -136.738756

对应于以下形式的函数:

c(x)=-6837615134.799868+12834646330.586414*x**1-9973474377.668280*x**2+3904659818.834625*x**3-633282611.288889*x**4-76066283.747942*x**5+49670960.939126*x**6-5091123.449217*x**7-1426628.818192*x**8+515175.778491*x**9-71055.177018*x**10+4863.969973*x**11-137.065848*x**12

可以在这里检查 c(x) 的样子:

在这样的图像中，a(x) 和 b(x) 是使用“多项式拟合”拟合的函数，仅针对少数几个点(4 和 7)。

那么，关于我在这里做错了什么有什么提示吗？还有一些其他的 C 代码可以提供合适的配件吗？

Answer 1

这两种解决方法之间的一个主要区别是，当您使用 gnuplot 时，您正在执行 fit 以在同一程序中设置系数并从那里绘制函数，而使用 GSL您正在将数字从一个程序复制到另一个程序。

如果您使用 printf("%lf", ...) 的输出作为您的第二个 gnuplot 程序的输入，您会失去很多准确性，因为 printf 对数字的舍入方式比任何一个程序的任何内部操作都要多。因为这是一个数值不稳定的问题，所以一点点四舍五入会造成很大的伤害。

当 x 为 4.877263 时，x**12 大约为 181181603.850932 所以如果您的 m 偏离 0.000001(printf 的默认舍入级别)引入了 181.181603850932 的错误，相对于 x 的实际 y 值大约 300% 的相对错误。

尝试 %.60lf，看看它是否变得更好。

如果其中一个程序在内部使用 long double 而另一个程序没有，那么无论您做什么，您都可能无法获得良好的匹配。