拉格朗日插值学习总结

简介

在数值分析中,拉格朗日插值法是以法国18世纪数学家约瑟夫·拉格朗日命名的一种多项式插值方法。如果对实践中的某个物理量进行观测,在若干个不同的地方得到相应的观测值,拉格朗日插值法可以找到一个多项式,其恰好在各个观测的点取到观测到的值。上面这样的多项式就称为拉格朗日(插值)多项式。

题目大意

给出$n$个点,将过这$n$个点的最多$n-1$次的多项式记为$f$$left(x ight)$,求$f$$left(k ight)$的值。

P4781 【模板】拉格朗日插值

拉格朗日插值法

一个最显然的思路就是直接高斯消元求出多项式的系数,但是这样做复杂度巨大$Oleft(n^3 ight)$且根据算法实现不同往往会存在精度问题。

而拉格朗日插值法可以在$Oleft(n^2 ight)$的复杂度内完美解决上述问题。

如图所示,将每一个点$left(x_i,y_i ight)$在$x$轴上的投影$left(x_i,0 ight)$记为$H_i$。对于每一个$i$,我们选择一个点集${P_i}cup{H_j|1le ile n,j e i}$,作过这$n$个点的至多$n-1$次的线$g_ileft(x ight)$。图中$fleft(x ight)$用黑线表示,$g_ileft(x ight)$用彩色线表示。

这样,我们得到了$n$个$g_ileft(x ight)left(1le ile n ight)$,它们都在各自对应的$x_i$处的值$y_i$为,并且在其它$x_jleft(j e i ight)$处值为$0$。所以很容易构造出$g_ileft(x ight)$的表达式:

$$g_ileft(x ight)=y_iprodlimits_{j e i}frac{x-x_j}{x_i-x_j}$$

令$$ell_ileft(x ight)=prodlimits_{j e i}frac{x-x_j}{x_i-x_j}$$

每个$ell_ileft(x ight)$就是拉格朗日基本多项式(或称插值基函数

拉格朗日基本多项式$ell_ileft(x ight)$的特点是在$x_i$上取值为$1$,在其他的点$x_j,j e i$上取值为$0$。

最后,我们所求的$fleft(x ight)=sumlimits_{i=1}^ng_ileft(x ight)$,即各$gleft(x ight)$之和。因为对于每一个$i$,都只有一条$g_i$经过$P_i$,其余$g_j$都经过$H_i$,故它们相加后在$x_i$的取值仍为$y_i$,即最后的和函数总是过所有$P_i$的。

公式整理即得拉格朗日插值法:

$$fleft(x ight)=sumlimits_{i=1}^ny_iprodlimits_{j e i}frac{x-x_j}{x_i-x_j}$$

对于本题来说,只用求出$fleft(k ight)$的值,直接带入公式求即可。

复杂度:$Oleft(n^2 ight)$

 Talk is cheap,show you code.

 1 int n, k;
 2 int a, b, ans;
 3 int x[N ], y[N ];
 4 inline int ksm (int a, int b)
 5 {
 6     int res = 1;
 7     while(b)
 8     {
 9         if(b & 1) res = res * a % P ;
10         a = a * a % P ;
11         b >>= 1;
12     }
13     return res;
14 }
15 inline int inv(int x)
16 {
17     return ksm (x, P - 2);
18 }
19 signed main()
20 {
21     read(n); read(k);
22     for(R int i = 1; i <= n; i++) read(x[i]), read(y[i]);
23     for(R int i = 1; i <= n; i++)
24     {
25         a = y[i] % P ;
26         b = 1;
27         for(R int j = 1; j <= n; j++)
28         {
29             if(i == j) continue;
30             a = a * (k - x[j]) % P ;
31             b = b * (x[i] - x[j]) % P ;
32         }
33         ans = (ans + a * inv(b) % P ) % P ;
34     }
35     writeln((ans + P ) % P );
36     return 0;
37 }

(% P + P) % P 熟练得令人心疼(懂的都懂

优点与缺点

拉格朗日插值法的公式结构整齐紧凑,在理论分析中十分方便,然而在计算中,当插值点增加或减少一个时,所对应的基本多项式就需要全部重新计算,于是整个公式都会变化,非常繁琐。这时可以用重心拉格朗日插值法或牛顿插值法来代替。此外,当插值点比较多的时候,拉格朗日插值多项式的次数可能会很高,因此具有数值不稳定的特点,也就是说尽管在已知的几个点取到给定的数值,但在附近却会和“实际上”的值之间有很大的偏差(如下图)。这类现象也被称为龙格现象,解决的办法是分段用较低次数的插值多项式。

重心拉格朗日插值法

$$fleft(x ight)=sumlimits_{i=1}^n y_i prodlimits_{j e i}frac{x-x_j}{x_i-x_j}$$

$$=sumlimits_{i=1}^n y_i frac{prodlimits_{j e i}left(x-x_j ight)}{prodlimits_{j e i}left(x_i-x_j ight)}$$

$$=sumlimits_{i=1}^n y_i frac{prodlimits_{i=1}^n left(x-x_i ight)}{prodlimits_{j e i}left(x_i-x_j ight)*left(x-x_i ight)}$$

$$=prodlimits_{i=1}^n left(x-x_i ight)sumlimits_{i=1}^n frac{y_i}{prodlimits_{j e i}left(x_i-x_j ight)*left(x-x_i ight)}$$

令$g=prodlimits_{i=1}^n left(x-x_i ight)$

令$wleft(i ight)=prodlimits_{j e i} x_i-x_j$

那么柿子变成$$fleft(x ight)=gsumlimits_{i=1}^n frac{y_i}{left(x-x_i ight)wleft(i ight)}$$

对于增加的插值点,我们只用$Oleft(n ight)$更新所有的$wleft(i ight)$。

对于询问求值,我们先$Oleft(n ight)$求出$g$,然后再$Oleft(n ight)$套公式即可。

记得先预处理逆元。不然白送一个$log_2 n$

Talk is cheap,show you code.

 1 int n, opt, x0, y0;
 2 int x[N ], y[N ], w [N ], num ;
 3 inline int ksm (int a, int b)
 4 {
 5     int res = 1;
 6     while(b)
 7     {
 8         if(b & 1) res = res * a % P ;
 9         a = a * a % P ;
10         b >>= 1;
11     }
12     return res;
13 }
14 inline int inv(int x)
15 {
16     return ksm (x, P - 2); 
17 }
18 inline void insert(int x0, int y0)
19 {
20     num ++;
21     x[num ] = x0;
22     y[num ] = y0;
23     w [num] = 1;
24     for(R int i = 1; i <= num - 1; i++)
25     {
26         w [i] = w [i] * (x[i] - x0) % P ;
27         w [num ] = w [num ] * (x0 - x[i]) % P ;
28     }
29 }
30 inline int work(int x0)
31 {
32     int g = 1, ans = 0;
33     for(R int i = 1; i <= num ; i++)
34     {
35         if(x[i] == x0) return y[i];
36         g = g * (x0 - x[i]) % P ;
37     }
38     for(R int i = 1; i <= num ; i++)
39         ans = (ans + y[i] * inv((x0 - x[i]) * w [i])) % P ;//我来白送
40     return (g * ans % P + P ) % P ;    
41 }
42 signed main()
43 {
44     read(n);
45     for(R int i = 1; i <= n; i++)
46     {
47         read(opt);
48         if(opt == 1)
49         {
50             read(x0); read(y0);
51             x0 %= P ; y0 %= P ;
52             insert(x0, y0);
53         }
54         else
55         {
56             read(x0);
57             writeln(work(x0));
58         }
59     }
60     return 0;
61 }

参考资料

拉格朗日插值

拉格朗日插值学习小结

拉格朗日插值法

原文地址:https://www.cnblogs.com/e-e-thinker/p/13394015.html