使用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代码来获得这样的拟合.对于几种情况,GNU Scientific Library的多项式拟合代码(如此链接)结果正确.但是对于上面的数据(以及我的数据集中的其他几个案例),我得到的结果是有缺陷的.

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

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)的样子:

在这样的图像中,(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 的错误，该错误相对于y该 的实际值约为 300% x。

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

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