拉格郎日插值法的计算精度完全丢失

涩兔子

Lagrange Interpolating拉格郎日插值法仅仅用线性逼近非线性，但是，我的方法让计算精度完全丢失了，请大家帮帮看看

SVN地址：
https://sirch.svn.sourceforge.ne ... -analysis/lagrange/

1. 
   
2. #include <stdio.h>
   
3. #include <stdlib.h>
   
4. #include <string.h>
   
5. 
   
6. #define DEBUG 1
   
7. 
   
8. static float m_p = 0.5635;
   
9. static char *m_xi = NULL;
   
10. static char *m_yi = NULL;
    
11. static float m_arr[] = {};
    
12. static unsigned int m_n = 0;
    
13. 
    
14. static void m_lagrange_interpolating_usage();
    
15. static int m_lagrange_interpolating_set_n_arr();
    
16. static float m_lagrange_interpolating();
    
17. 
    
18. static void
    
19. m_lagrange_interpolating_usage()
    
20. {
    
21.   printf("Lagrange Interpolating Usage:\n");
    
22.   printf("bin/li -p XXX -x XXX -y XXX\n");
    
23.   printf("-p 0.56350\n");
    
24.   printf("-x 0.56160,0.56280,0.56401,0.56521\n");
    
25.   printf("-y 0.82741,0.82659,0.82577,0.82495\n");
    
26. }
    
27. 
    
28. static int
    
29. m_lagrange_interpolating_set_n_arr()
    
30. {
    
31.   int xi_n = 0;
    
32.   int yi_n = 0;
    
33.   int n = 0;
    
34.   char *token = NULL;
    
35.   const char *delim = ",";
    
36. 
    
37.   token = strtok(m_xi, delim);
    
38.   while (token)
    
39.   {
    
40.     xi_n++;
    
41.     n++;
    
42.     m_arr[n] = (float) atof(token);
    
43.     token = strtok(NULL, delim);
    
44.   }
    
45. 
    
46.   token = strtok(m_yi, delim);
    
47.   while (token)
    
48.   {
    
49.     yi_n++;
    
50.     n++;
    
51.     m_arr[n] = (float) atof(token);
    
52.     token = strtok(NULL, delim);
    
53.   }
    
54.   
    
55.   if (xi_n != yi_n)
    
56.   {
    
57.     return 0;
    
58.   }
    
59. 
    
60.   m_n = xi_n = yi_n;
    
61. 
    
62.   return 1;
    
63. }
    
64. 
    
65. static float
    
66. m_lagrange_interpolating()
    
67. {
    
68.   float li = 0.00000;
    
69.   float l = 0.00000;
    
70.   int i;
    
71.   int j;
    
72. 
    
73. #if DEBUG
    
74.   printf("DEBUG at %d: debug mode\n", __LINE__);
    
75.   unsigned int m_n = 4;
    
76.   float m_p = 0.5635;
    
77.   float m_x[] = {0.56160, 0.56280, 0.56401, 0.56521};
    
78.   float m_y[] = {0.82741, 0.82659, 0.82577, 0.82495};
    
79. 
    
80.   for (i = 0; i < m_n; i++)
    
81.   {
    
82.     for (j = 0; j < m_n; j++)
    
83.     {
    
84.       if (j != i)
    
85.       {
    
86.         li = li * (m_p - m_x[j]) / (m_x[i] - m_x[j]);
    
87.       }
    
88.     }
    
89.     l = l + li * m_y[i];
    
90.   }
    
91. #else
    
92.   if (!m_lagrange_interpolating_set_n_arr())
    
93.   {
    
94.     printf("Error: xi and yi counter is diff\n");
    
95.    
    
96.     return 0;
    
97.   }
    
98. 
    
99.   for (i = 0; i < m_n; i++)
    
100.   {
     
101.     for (j = 0; j < m_n; j++)
     
102.     {
     
103.       if (j != i)
     
104.       {
     
105.         li = li * (m_p - m_arr[j + 1]) / (m_arr[i + 1] - m_arr[j + 1]);
     
106.       }
     
107.     }
     
108.     l = l + li * m_arr[i + m_n + 1];
     
109.   }
     
110. #endif
     
111. 
     
112.   return l;
     
113. }
     
114. 
     
115. int
     
116. main(int argc, char **argv)
     
117. {
     
118.   int c;
     
119.   extern char *optarg;
     
120.   float res;
     
121. 
     
122.   if (argc < 2)
     
123.   {
     
124.     m_lagrange_interpolating_usage();
     
125.    
     
126.     return -1;
     
127.   }
     
128. 
     
129.   while ((c = getopt(argc, argv, "p:x:y:")) != -1)
     
130.   {
     
131.     switch (c)
     
132.     {
     
133.       case 'p':
     
134.               m_p = (float) atof(optarg);
     
135.       case 'x':
     
136.         m_xi = optarg;
     
137.         break;
     
138.       case 'y':
     
139.         m_yi = optarg;
     
140.         break;
     
141.     }
     
142.   }
     
143. 
     
144.   res = m_lagrange_interpolating();
     
145.   printf("Lagrange interpolating with %.5f result is %.5f\n", m_p, res);
     
146.   
     
147.   return 0;
     
148. }
     

复制代码

涩兔子

啊啊啊，天，公式中的Li的初始值搞错了

1. 
   
2.    for (i = 0; i < m_n; i++)
   
3.   {
   
4.     li = 1.0000000;
   
5.     for (j = 0; j < m_n; j++)
   
6.     {
   
7.       if (j != i)
   
8.       {
   
9.         li = li * (m_p - m_arr[j + 1]) / (m_arr[i + 1] - m_arr[j + 1]);
   
10.       }
    
11.     }
    
12.     l = l + li * m_arr[i + m_n + 1];
    
13.   }
    

复制代码

https://sirch.svn.sourceforge.ne ... ange-interpolating/