数值计算——牛顿插值法

置顶 zsyz_ZZY

已于 2023-04-14 19:58:48 修改

阅读量1.8k

点赞数 2

分类专栏：数论文章标签：抽象代数

于 2023-04-06 22:00:39 首次发布

本文链接：https://blog.csdn.net/zsyz_ZZY/article/details/129999636

版权

数论专栏收录该内容

4 篇文章 0 订阅

订阅专栏

一般牛顿插值

①考虑两个点的情况， $f_1(x)=f(x_0)+k_1(x-x_0)$ ，带入 $x_1,f(x_1))$ ，得
$k_1=\frac{f(x_1)-f(x_0)}{x_1-x_0}\\ \therefore f_1(x)=f(x_0)+\frac{f(x_1)-f(x_0)}{x_1-x_0}(x-x_0)$
②考虑三个点的情况， $f_2(x)=f_1(x)+k_2(x-x_0)(x-x_1)$ ，带入 $x_2,f(x_2))$ ，得
$\begin{aligned}k_2&=\frac{f(x_2)-f_1(x_2)}{(x_2-x_0)(x_2-x_1)}\\ &=\frac{[f(x_2)-f(x_0)]-[\frac{f(x_1)-f(x_0)}{x_1-x_0}(x_2-x_0)]}{(x_2-x_0)(x_2-x_1)}\\ &=\frac{[f(x_2)-f(x_0)](x_1-x_0)-[f(x_1)-f(x_0)](x_2-x_0)}{(x_2-x_0)(x_2-x_1)(x_1-x_0)}\\ &=\frac{[f(x_2)-f(x_1)](x_1-x_0)-[f(x_1)-f(x_0)](x_2-x_1)}{(x_2-x_0)(x_2-x_1)(x_1-x_0)}\\ &=\frac{\frac{f(x_2)-f(x_1)}{x_2-x_1}-\frac{f(x_1)-f(x_0)}{x_1-x_0}}{x_2-x_0} \end{aligned}\\ \therefore f_2(x)=f(x_0)+\frac{f(x_1)-f(x_0)}{x_1-x_0}(x-x_0)+\frac{\frac{f(x_2)-f(x_1)}{x_2-x_1}-\frac{f(x_1)-f(x_0)}{x_1-x_0}}{x_2-x_0}(x-x_0)(x-x_1)$

我们定义
$\begin{aligned}&f[x_0,x_1]=\frac{f(x_1)-f(x_0)}{x_1-x_0}\\ &f[x_0,x_1,x_2]=\frac{f[x_1,x_2]-f[x_0,x_1]}{x_2-x_0}\\ &...\\ &f[x_0,x_1,...,x_k]=\frac{f[x_1,x_2,...,x_k]-f[x_0,x_1,...,x_{k-1}]}{x_k-x_0} \end{aligned}$

其中 $f[x_0,x_1]$ 为 $f (x)$ 在 $x_0,f(x_0)),(x_1,f(x_1))$ 的一阶差商，
$f[x_0,x_1,x_2]$ 为 $f (x)$ 在 $x_0,f(x_0)),(x_1,f(x_1)),(x_2,f(x_2))$ 的二阶差商,
…
$f[x_0,x_1,...,x_k]$ 为 $f (x)$ 在 $x_0,f(x_0)),(x_1,f(x_1)),...,(x_k,f(x_k))$ 的 $k$ 阶差商。

归纳一下，得
$\therefore f[x_0,x_1,...,x_k]=\sum_{i=0}^{k}\frac{f(x_i)}{\prod_{j=0,j\not =i}^{k}(x_i-x_j)}\\ \therefore f(x)=f(x_0)+f[x_0,x_1](x-x_0)+f[x_0,x_1,x_2](x-x_0)(x-x_1)+...+f[x_0,x_1,...,x_k](x-x_0)(x-x_1)...(x-x_{k-1})$
code

#include<cstdio>
#include<cstring>
#include<algorithm>
#define LL long long
#define mod 998244353
using namespace std;
	int n;
	LL k;
	struct node{LL x,y;} a[2010];
	LL times[2010][2010],delta[2010];
LL dg(LL x,LL k)
{
	if(!k) return 1;
	LL op=dg(x,k>>1);
	if(k&1) return op*op%mod*x%mod; else return op*op%mod;
}
LL inv(LL x)
{
	return dg(x,mod-2);
}
void init()
{
	for(int i=0;i<n;i++)
	{
		if(i) times[i][0]=(a[i].x-a[0].x+mod)%mod; else times[i][0]=1;
		for(int j=1;j<n;j++)
			if(i!=j) times[i][j]=times[i][j-1]*((a[i].x-a[j].x+mod)%mod)%mod; else times[i][j]=times[i][j-1];
	}
	for(int i=1;i<n;i++)
		for(int j=0;j<=i;j++)
			delta[i]=(delta[i]+a[j].y*inv(times[j][i])%mod)%mod;
}
LL newton(LL k)
{
	LL ans=a[0].y,tmp=1;
	for(int i=1;i<n;i++)
	{
		tmp=(tmp*((k-a[i-1].x+mod)%mod))%mod;
		ans=(ans+delta[i]*tmp%mod)%mod;
	}
	return ans;
}
int main()
{
	scanf("%d %lld",&n,&k);
	for(int i=0;i<n;i++)
		scanf("%lld %lld",&a[i].x,&a[i].y);
	init(); 
	printf("%lld",newton(k));
}

等间距牛顿插值

规定取点的间距为固定值 $\Delta x$ ，即满足 $x_1=x_0+\Delta x,x_2=x_1+\Delta x,...,x_k=x_{k-1}+\Delta x$ 。
我们定义
$\begin{aligned}&\Delta f(x_0)=f(x_0+\Delta x)-f(x_0)\\ &\Delta^2 f(x_0)=\Delta f(x_0+\Delta x)-\Delta f(x_0)\\ &...\\ &\Delta^kf(x_0)=\Delta^{k-1}f(x_0+\Delta x)-\Delta ^{k-1}f(x_0) \end{aligned}$

其中 $\Delta f(x_0)$ 为 $f (x)$ 在 $x_0,f(x_0))$ 的一阶向前差分，
$\Delta^2 f(x_0)$ 为 $f (x)$ 在 $x_0,f(x_0))$ 的二阶向前差分，
…
$\Delta^k f(x_0)$ 为 $f (x)$ 在 $x_0,f(x_0))$ 的 $k$ 阶向前差分。
$\begin{aligned} \therefore &f_1(x)=f(x_0)+\frac{\Delta f(x_0)}{\Delta x}(x-x_0)\\ &f_2(x)=f(x_0)+\frac{\Delta f(x_0)}{\Delta x}(x-x_0)+\frac{\Delta^2f(x_0)}{2\Delta^2 x}(x-x_0)(x-x_1)\\ &...\\ &f_k(x)=f(x_0)+\frac{\Delta f(x_0)}{\Delta x}(x-x_0)+\frac{\Delta^2f(x_0)}{2\Delta^2 x}(x-x_0)(x-x_1)+...+\frac{\Delta^kf(x_0)}{k!\Delta^k x}(x-x_0)(x-x_1)...(x-x_{k-1}) \end{aligned}$

证明同差商。
若 $\Delta x\rightarrow0$ ，因此 $\Delta x_0=\Delta x_1=...=\Delta x_k$ ，则有
$\begin{aligned} &f_1(x)=f(x_0)+f'(x_0)(x-x_0)\\ &f_2(x)=f(x_0)+f'(x_0)(x-x_0)+f''(x_0)(x-x_0)^2\\ &...\\ &f_k(x)=f(x_0)+f'(x_0)(x-x_0)+\frac{f''(x_0)}{2!}(x-x_0)^2+...+\frac{f^{(k)}(x_0)}{k!}(x-x_0)^k \end{aligned}$
正是泰勒展开从 $1$ 到 $k$ 阶的展开式。
而在此意义下，牛顿差值的误差 $R_k(x)=\frac{f^{(k+1)}(\xi)}{(k+1)!}(x-x_0)^{k+1}(x_0<\xi<x)$ 也正是泰勒展开的拉格朗日余项。
也就是说，严格意义上， $f(x)=f_k(x)+R_k(x)$ ，但是在一般计算时不会考虑 $R_k(x)$ ，因此我们一般认为 $f(x)=f_k(x)$ 。

牛顿向前&向后插值

无论是向前插值还是向后插值，都是等间距的插值。

向前插值

其实上面讲的等间距插值就是用的向前插值。
令 $h=\Delta x,x=x_0+hn$ ，有 $x-x_0=hn,x-x_1=h(n-1),...,x-x_{k-1}=h(n-k+1)$ ，上面的式子变为
$\begin{aligned} &f_1(x)=f(x_0)+\frac{\Delta f(x_0)}{h}hn\\ &f_2(x)=f(x_0)+\frac{\Delta f(x_0)}{h}hn+\frac{\Delta^2f(x_0)}{2h^2}hn\cdot h(n-1)\\ &...\\ &f_k(x)=f(x_0)+\frac{\Delta f(x_0)}{h}hn+\frac{\Delta^2f(x_0)}{2h^2}hn\cdot h(n-1)+...+\frac{\Delta^kf(x_0)}{k!h^k}hn\cdot h(n-1)\cdot...\cdot h(n-k+1) \end{aligned}$

化简一下，得
$\begin{aligned} &f_1(x)=f(x_0)+\frac{n}{1}\Delta f(x_0)\\ &f_2(x)=f(x_0)+\frac{n}{1}\Delta f(x_0)+\frac{n(n-1)}{2}\Delta^2f(x_0)\\ &...\\ &f_k(x)=f(x_0)+\frac{n}{1}\Delta f(x_0)+\frac{n(n-1)}{2}\Delta^2f(x_0)+...+\frac{n(n-1)...(n-k+1)}{k!}\Delta^kf(x_0)\end{aligned}$

这就是前插公式。
其余项
$R_k(x)=\frac{n(n-1)...(n-k)}{(k+1)!}h^{k+1}f^{k+1}(\xi),\xi∈(x_0,x_k)$

向后插值

同理，我们定义 $\nabla f(x_0)$ 为 $f (x)$ 在 $x_0,f(x_0))$ 的一阶向后差分，
$\nabla^2 f(x_0)$ 为 $f (x)$ 在 $x_0,f(x_0))$ 的二阶向后差分，
…
$\nabla^k f(x_0)$ 为 $f (x)$ 在 $x_0,f(x_0))$ 的 $k$ 阶向后差分。
$\begin{aligned}&\nabla f(x_0)=f(x)-f(x_0-\Delta x)\\ &\nabla^2 f(x_0)=\nabla f(x_0)-\nabla f(x_0-\Delta x)\\ &...\\ &\nabla^kf(x_0)=\nabla^{k-1}f(x_0)-\nabla^{k-1}f(x_0-\Delta x) \end{aligned}$

同理，令 $h=\Delta x,x=x_0+hn$ ，有 $x-x_0=hn,x-x_1=h(n-1),...,x-x_{k-1}=h(n-k+1)$ 。
后插公式如下：
$\begin{aligned} &f_1(x)=f(x_0)+\frac{n}{1}\Delta f(x_0)\\ &f_2(x)=f(x_0)+\frac{n}{1}\Delta f(x_0)+\frac{n(n+1)}{2}\Delta^2f(x_0)\\ &...\\ &f_k(x)=f(x_0)+\frac{n}{1}\Delta f(x_0)+\frac{n(n+1)}{2}\Delta^2f(x_0)+...+\frac{n(n+1)...(n+k-1)}{k!}\Delta^kf(x_0)\end{aligned}$