浅谈拉格朗日插值

浅谈拉格朗日插值

数值分析中,拉格朗日插值法是以法国十八世纪数学家约瑟夫·拉格朗日命名的一种多项式插值方法。许多实际问题中都用函数来表示某种内在联系或规律,而不少函数都只能通过实验和观测来了解。拉格朗日插值法可以找到一个多项式,其恰好在各个观测的点取到观测到的值。这样的多项式称为拉格朗日(插值)多项式。

——百度百科

通俗地说,拉格朗日插值法可以找出一个恰好经过直角坐标系内n个给定点的函数

众所周知,n个点能唯一确定的多项式最高次数是n1

这个可以两点确定一次函数,三点确定二次函数来推出,或者我们由方程组有唯一解的充要条件也能得出

知道了这个之后我们就容易想到直接用高斯消元来搞出多项式的系数,但是复杂度消耗太大,是O(n3)

而拉格朗日插值法就是一种一般可以在O(n2)的复杂度下求出多项式的方法(不过通常用来求一个点值,以下的讲述一般以此为参考)


一般的拉格朗日插值法

我们设要求的多项式为f(x),点的坐标为(xi,yi),我们要找多项式在x0处的取值

先上公式,由拉格朗日插值法:

f(x0)=i=1nyiijx0xjxixj

看起来不是很好理解,其实很简单,我们把原来的某个给定点xk代入以下有:

f(xk)=i=1nyiijxkxjxixj

容易发现,当ki时,后面的xkxjxixj的分子总有一项是0,此时ijxkxjxixj=0

k=i时,后面的xkxjxixj上下完全相同,此时ijxkxjxixj=1

即对于f(xk)来说,这个多项式的确给出了对应的yk的值

不难发现这个方法对所有点都适用,因此它是正确的

从上面的式子可以看出每次计算要枚举两次,因此复杂度很简单,就是O(n2)


x取值连续时的插值法

因为很多时候我们做题都是先发现某个函数是多少次的多项式,然后自己随意取一些值代入插值

这样的话为了省事横坐标的取值完全可以从1开始连续取,那么我们把上面的式子中的xi换成i就有:

f(x0)=i=1nyiijx0jij

考虑怎么快速求ijkjij,我们分别考虑:

分子的话容易发现就是t=1nx0tx0i

分母比较复杂,ij的累乘可以分成两个阶乘部分,因此推导一下就是(1)nii!(ni)!

这样我们一般就可以O(nlogn)来算了(log主要是有求逆元的过程,当然你预处理O(n)求也没有问题)


重心拉格朗日插值法

我们考虑到朴素的拉格朗日插值每次多加入一个点时就要整个重新算过,很浪费时间

那么能不能把重复算的一些东西利用起来?

我们对于

f(x0)=i=1nyiijx0xjxixj

把分子提取出来,设为g=i=1nx0xi,则此时:

f(x0)=gi=1nijyi(x0xi)(xixj)

设一个ti=yiijxixj,则:

f(x0)=gi=1ntix0xi

因此每次多加入一个点只需要重新O(n)算它的ti就好了


拉格朗日插值的应用以及常用解题思路

一个经典例子:自然数幂和,即求i=1nik

之前也提到过用第二类斯特林数做的方法,但那种方法是O(k2)的,不够优秀

但是现在我们观察这个式子,如果不看求和的话ik就是一个k次多项式

那么前缀和(其实就是差分)之后,次数要+1,即此时答案为一个关于nk+1多项式

那我们直接代入k个值(取值连续,反正自己定)之后插值算即可,复杂度O(klogk)

那么很多具体题目怎么办呢,一般就是先推出某个式子,然后证明它是关于某个自变量的多少次函数(一定要判断出次数!),然后自己选一些点(一般是连续的)代入插值即可

其实做了一些题目之后就很套路了


拉格朗日插值的模板

模板看Luogu P4781 【模板】拉格朗日插值 ,就是用一般的O(n2)来做就好了

CODE

#include<cstdio>
#define RI register int
#define CI const int&
using namespace std;
const int N=2005,mod=998244353;
int n,x[N],y[N],k;
inline int sub(CI x,CI y)
{
	int t=x-y; return t<0?t+mod:t;
}
inline int inv(int x,int p=mod-2,int mul=1)
{
	for (;p;p>>=1,x=1LL*x*x%mod) if (p&1) mul=1LL*mul*x%mod; return mul;
}
inline int Lagrange(CI n,int *x,int *y,CI k)
{
	int ret=0; for (RI i=1;i<=n;++i)
	{
		int s1=1,s2=1; for (RI j=1;j<=n;++j) if (i!=j)
		s1=1LL*s1*sub(k,x[j])%mod,s2=1LL*s2*sub(x[i],x[j])%mod;
		(ret+=1LL*y[i]*s1%mod*inv(s2)%mod)%=mod;
	}
	return ret;
}
int main()
{
	scanf("%d%d",&n,&k); for (RI i=1;i<=n;++i) scanf("%d%d",&x[i],&y[i]);
	return printf("%d",Lagrange(n,x,y,k)),0;
}

一些入门例题

posted @   空気力学の詩  阅读(5019)  评论(4编辑  收藏  举报
编辑推荐:
· 如何编写易于单元测试的代码
· 10年+ .NET Coder 心语,封装的思维:从隐藏、稳定开始理解其本质意义
· .NET Core 中如何实现缓存的预热?
· 从 HTTP 原因短语缺失研究 HTTP/2 和 HTTP/3 的设计差异
· AI与.NET技术实操系列:向量存储与相似性搜索在 .NET 中的实现
阅读排行:
· 地球OL攻略 —— 某应届生求职总结
· 周边上新:园子的第一款马克杯温暖上架
· Open-Sora 2.0 重磅开源!
· 提示词工程——AI应用必不可少的技术
· .NET周刊【3月第1期 2025-03-02】
点击右上角即可分享
微信分享提示