Numerical Approximation Chapter 5 Notes

Posted 2023-04-08 WenDavid的数学随笔

考虑在某个范围\\([a,b]\\)上的\\(n+1\\)个点值所确定的插值\\(n\\)次多项式\\(p\\)满足

\\[p(x_i)=f(x_i) \\]

其中\\(f\\in \\mathscr L[a,b]\\) , \\(p\\)是\\(\\mathscr P_n\\)中满足条件的polynomial. 考虑到按照条件可以列出一个\\(n+1\\)元方程组，未知数表示对应的系数，这里多项式系数空间的形态取决于对应方程组的形态，考虑到\\(x_i\\) distinct可以证明\\(\\beginbmatrix1&x_i&x_i^2&\\cdots&x_i^n\\endbmatrix\\)之间linear independent，所以\\(n\\)次多项式对应的解只有一个，这是trivial的。这样的插值多项式中\\(x^n\\)的系数被定义为divided difference of order \\(n\\)，记为\\(f[x_0,x_1,\\dots,x_n]\\).注意这只是一个系数，不是一个函数.

一个\\(n\\)次多项式使用\\(n+1\\)个点值，这是正确的。实际上我们还可以想到比如Discrete Fourier Transform里面的\\(n\\)个单位圆上的点值，在复平面上面延拓的多项式插值拥有更好的性质.

其形式其实就是Lagrange Interpolation的最高项. 考虑Lagrange Interpolation

\\[p=\\sum\\limits_k=0^nf(x_k)\\prod\\limits_j\\neq k\\fracx-x_jx_k-x_j \\]

那么最高项就是\\(\\sum\\limits\\)里面每一项把\\(x\\)乘一起，显然

\\[f[x_0,x_1,\\dots,x_n]=\\sum\\limits_k=0^n\\fracf(x_k)\\prod\\limits_j\\neq k(x_k-x_j) \\]

这个divided difference是关于\\(f(x_i)\\)的linear function，但是实际中直接计算容易出现巨大误差(e.g. 考虑\\(n\\)比较大，同时\\(x_j\\)的误差会通过倒数快速放大)，后面有一个使用递推的法子。

插值的中值定理

Theorem 5.1 插值的高次中值定理误差项
考虑\\(f\\in\\mathscr L^(n)[a,b]\\)，也就是在\\([a,b]\\)上\\(n\\)次可导，\\(x_i\\)是\\([a,b]\\)上面选的一堆点，则存在\\(\\zeta\\) in the smallest interval that contains points \\(x_i\\) such that

\\[f[x_0,x_1,\\dots,x_n]=\\fracf^(n)(\\zeta)n! \\]

• \\(f\\)在\\([a,b]\\)上\\(n\\)次可导，考虑一个Taylor展开到\\(n\\)次，意思就是Taylor展开的\\(x^n\\)系数.
• 因为\\(x^n\\)到头了（指Taylor展开到\\(x^n\\)后面的是高阶无穷小），所以直接把\\(f\\)求导即可.
  Proof. 考虑\\(e(x)=f(x)-p(x)\\),\\(x\\in [a,b]\\). 因为\\(p(x)\\)是interpolation polynomial，显然\\(e(x)\\in\\mathscr L^(n)[a,b]\\). 考虑一堆\\(x_i\\)把\\([a,b]\\)分成若干段，这些段的边界上考虑interpolation condition，\\(e(x_i)=0\\)，由罗尔定理（用\\(n\\)次）在这些每段里面做一下，每段都存在\\(\\zeta\\)使得\\(e^(n)(\\zeta)=0\\). 也就是说

\\[p^(n)(\\zeta)=f^(n)(\\zeta) \\]

考虑\\(p\\)这个\\(n\\)次多项式求\\(n\\)阶导之后是\\(f[x_0,x_1,\\dots,x_n]n!\\).得证.

这个定理有啥用？

This theorem is an important part of the standard method of checking tabulated values of a function for errors. 也就是说，验证一组数据是否符合某个模型。
考虑这样一种情况，一个\\(f\\in\\mathscr L^(n)[a,b]\\)，提供了\\(m\\)个点值来做一个插值。假如\\(n<<m\\),同时\\(m\\)在\\([a,b]\\)里面比较密集。实际上这就是用来验证的法子，\\(f\\)是给定模型的某个奇怪的函数，然后可能就是有\\(m\\)组数据的点值，我们需要验证这些点值是否符合\\(f\\)。假设\\(a\\leq x_0<x_1<\\dots<x_m\\leq b\\)，这里验证的方法是考虑一个滑动的窗口，在这里面取一段连续的\\(n\\)个点，例如\\(x_j,x_j+1,\\dots,x_j+n\\)这样，因为\\(n<<m\\)，\\([x_j,x_j+n]\\)其实很短，所以\\(f^(n)\\)几乎没怎么变，也就是说\\(f[x_j,x_j+1,\\dots,x_j+n]\\)随着\\(j=0,1,\\dots,m-n\\)其实变化得非常慢。在后面我们会提到计算\\(f[x_j,x_j+1,\\dots,x_j+n]\\)的方法，假如计算出来之后发现某些地方的点不太光滑，那么就表明有一些数据不符合\\(f\\)，至少相当于是在\\(n\\)阶这个级别上的增长趋势不太符合。似乎我们可以考虑一步步验证\\(n=0,1,\\dots,n\\)这样.

Newton\'s interpolation method

考虑我们有\\(n+1\\)个点值插值出的\\(n\\)次的多项式，这个时候来多一个点值会有什么变化呢？

\\[p_n=\\sum\\limits_k=0^nf(x_k)\\prod\\limits_j\\neq k\\fracx-x_jx_k-x_j \\]

Theorem 5.2 插值多项式的递推关系

\\[p_n+1(x)=p_n(x)+\\left\\\\prod\\limits_j=0^n(x-x_j)\\right\\f[x_0,x_1,\\dots,x_n+1] \\]

Proof. 有了结果证明其实并不困难，考虑右边往上加的时候，我们考虑\\(f[x_0,x_1,\\dots,x_n+1]=\\sum\\limits_k=0^n+1\\fracf(x_k)\\prod\\limits_j\\neq k(x_k-x_j)\\)和\\(p_n(x)\\)中\\(k\\)一样的项加在一起（\\(f[x_0,x_1,\\dots,x_n+1]\\)多一个)，形如

\\[\\beginaligned &\\fracf(x_k)\\prod\\limits_j\\neq k^n+1(x_k-x_j)\\prod_j=0^n(x-x_j)+f(x_k)\\prod_j\\neq k^n\\fracx-x_jx_k-x_j\\\\ \\endaligned \\]

这两边有以下不同：

• \\((x-x_j)\\)乘积，左边项多一个\\(j=k\\)的
• \\((x_k-x_j)\\)乘积，左边项多一个\\(j=n+1\\)的
  考虑\\(\\fracx-x_k(x_k-x_n+1)+1=\\fracx-x_n+1x_k-x_n+1\\)，所以其实相当于多乘了一项\\(n+1\\)对应的，非常巧妙！
  额外多的\\(\\left\\\\prod\\limits_j=0^n(x-x_j)\\right\\\\fracf(x_n+1)\\prod\\limits_j\\neq n+1(x_n+1-x_j)\\)实际上就是插值里面多的一项。

Divided differences的递归求解

上面提到的

• 比较精确地求解插值多项式系数的问题（使用递推）
• 验证数据是否符合某个模型
  都需要用到divided differences的计算。我们可以使用递归的法子来算。
  Theorem 5.3 divided differences递归求解

\\[f[x_j,\\dots,x_j+k+1]=\\fracf[x_j+1,\\dots,x_j+k+1]-f[x_j,\\dots,x_j+k]x_j+k+1-x_j \\]

note. 一个长度\\(k+1\\)的区间被分到两个\\(k\\)的区间求解，实际上有点像组合数，可以顺便思考一下能不能分到三个\\(k-1\\)的区间，以此类推（或者两个\\(k-1\\)的区间？）.
proof.1 直接代入公式，会发现相减之后除去公共的，剩下的就是\\(j\\)和\\(j+k+1\\)这两项在直接对决，所以要除掉额外的项，总之是可以推，不过另一个证明方法比较优雅.
proof.2 考虑\\(p_k+1(x)\\)的形式关于关于\\(f(x_j),\\dots,f(x_j+k+1)\\)总是唯一的，所以实际上可以考虑用\\(f[x_j+1,\\dots,x_j+k+1]\\)和\\(f[x_j,\\dots,x_j+k]\\)凑一个interpolation polynomial出来，只要次数一样，满足interpolation condition那就是一样的.分别将\\(f[x_j+1,\\dots,x_j+k+1]\\)和\\(f[x_j,\\dots,x_j+k]\\)对应的Lagrange Interpolation polynomial记为\\(q_k\\)和\\(p_k\\)，那么

\\[p_k+1(x)=\\fracx-x_jx_j+k+1-x_jq_k(x)+\\fracx-x_j+k+1x_j-x_j+k+1p_k(x) \\]

总之是基于Lagrange Interpolation的思想可以构造以上式子. 只需要关注往上乘\\(x\\)的最高次的变动就能得出原来的式子.

Hermite Interpolation

哈密顿插值？厄米特插值？

现在我们的插值不满足于\\(p(x_i)=f(x_i)\\)，我们要求其高阶导数也要一致，即\\(p^(j)(x_i)=f^(j)(x_i)\\).在不同的位置，我们需要拟合的导数阶数也不同，比如靠近零点附近点比较密集，导数阶数可以少一点，远处点比较稀疏，但是希望能拟合长远的趋势，需要拟合高阶的趋势，等等.
考虑点\\(x_i\\)需要拟合\\(j=0,\\dots,l_i\\)阶的导数，考虑到解方程，\\(p\\)的系数个数应该等于number of data，所以

\\[n+1=\\sum\\limits_i (l_i+1) \\]

Theorem 5.4 哈密顿插值的唯一性

• only one polynomial in \\(\\mathscr P_n\\) that satisfies the conditions where the \\(n\\) is defined above.

proof. 考虑\\(f_1-f_2\\)显然最多是\\(n\\)次多项式，然后这个多项式在\\(x_i\\)处up to \\(l_i\\)阶导数都是\\(0\\)，多项式有子式\\((x-x_i)^l_i+1\\).全乘一块发现有\\(x^n+1\\)，所以\\(f_1-f_2=0\\).

实际上这个用\\((x-x_i)^k\\)的思路已经可以比较轻松地构造导数了.导到某阶的时候，让剩下的都有\\((x-x_i)\\)的子式，只剩对应的导数的常数就可以了. 不过事实上Newton\'s Interpolation method仍然可以用在这个上面，不过需要一定的修改来恰好塞进高阶导数的约束. 考虑比如我们对\\(x_i\\)有\\(l_i\\)阶导数的约束，我们就把\\(x_i\\)重复写\\(l_i+1\\)次，然后做Newton\'s Interpolation method. 这个时候形如\\(x_0,x_0,,\\dots,x_0,x_1,x_1,\\dots,x_1,\\dots,x_m,x_m,\\dots,x_m\\). 这个时候会出现一个问题，我们Theorem 5.3考虑的是distinct points，而

\\[f[x_j,\\dots,x_j+k+1]=\\fracf[x_j+1,\\dots,x_j+k+1]-f[x_j,\\dots,x_j+k]x_j+k+1-x_j \\]

在\\(x_j+k+1=x_j\\)的时候不是well defined，这怎么办呢？我们使用Theorem 5.1，

\\[f[x_0,x_1,\\dots,x_n]=\\fracf^(n)(\\zeta)n! \\]

因为这里相当于\\(x_0=x_1=\\dots=x_n\\)的情况，所以这里的\\(f[x_j,x_j+1,\\dots,x_j+k+1]\\)可以直接写出为\\(\\fracf^(k+1)(x_i)(k+1)!\\).刚好可以把对应的导数的约束塞进去！以上这个法子被称为extension of Newton\'s method.
Theorem 5.5 Extension of Newton\'s method的正确性

• 用extension of Newton\'s method可以完成Hermite interpolation的任务.
  proof. 在\\(x_i\\)附近做Taylor展开到\\(l_i\\)其实就可以说明了.

本文来自博客园，作者：WenDavid552，转载请注明原文链接：https://www.cnblogs.com/wdmath/p/17297833.html

ARC-100 C - Linear Approximation

题面在这里！

 

    可以看成点集{a[i]-i}和b之间距离的和，于是找到中位数就可以直接算了2333.

 

#include<bits/stdc++.h>
#define ll long long
using namespace std;
const int N=200005;

int a[N],n,num;
ll ans=0;

int main(){
	scanf("%d",&n);
	for(int i=1;i<=n;i++) scanf("%d",a+i),a[i]-=i;
	sort(a+1,a+n+1),num=a[(n+1)>>1];
	for(int i=1;i<=n;i++) ans+=(ll)abs(num-a[i]);
	cout<<ans<<endl;
	return 0;
}