+++
date = '2024-11-30T15:26:27+08:00'
draft = true
title = '拉格朗日插值和数值微积分'
+++
初次发布于我的个人文档。（每次都是个人文档优先发布哦）
本文想简要介绍和推导一下拉格朗日插值和数值积分方法。
什么是插值？
所谓的插值就是已知几个离散点的信息视图求一个满足这些信息的函数的过程。
如拉格朗日插值和牛顿插值就是已知f(x)在若干点的函数值希望找一个多项式函数穿过这些点。
而艾尔米特插值则更进一步要求函数在各个点的各阶导数值都等于指定的值。
本文只介绍最基础的拉格朗日插值和牛顿插值。
插值的唯一性
这两种插值方法都是已知一系列点$(x_i,y_i)$，找一个多项式p(x)它穿过了这些点。
其实这个问题你应该很容易就能找到第一个方法，待定系数法嘛。
我们假设$p(x)=a_0+a_1x+a_2x+...+a_nx^n$
这里有n+1个系数，所以我们需要给定n+1个不同的点，或者说n+1个点需要用n次插值多项式插值。
如果我们将这n+1个点代入就能得到
$a_0+a_1x_i+a_2x_i+...+a_nx_i^n=y_i$
显然，这是一个线性方程组，并且还有更令人兴奋的，
它的系数矩阵的行列式就是范德蒙德行列式。
从而，其系数矩阵的行列式是$x_i-x_j$的乘积，又由于这n+1个点是不同的，所以它的系数矩阵的行列式不为0。
那么根据克莱姆法则，这个线性方程组就有唯一解。
所以插值多项式是唯一的。
拉格朗日插值
那么怎么求插值多项式呢？
你是可以硬解刚刚的方程组啦，但是这有点痛苦。
拉格朗日的方法是，这样的：
你不是要求$p(x_i)=y_i$吗？
如果我能找到一系列多项式$l_k(x)$在$x_k$处取1，其他地方（指$x_1,x_2,x_3,...,x_n$但不包含$x_k$)取0,
那么p(x)不就是$y_i l_k(x)$吗？
这就是拉格朗日插值法的思想。
其中$l_k(x)$被称为插值基函数。
那么我们怎么找到插值基函数呢？
其实很简单，$l_k(x)$在其他点取0从而其他点都是插值基函数的零点，所以插值基函数$l_k(x)$有因子
$$\prod_{i \neq k}(x-x_i)$$
那么我们怎么保证$l_k(x)$在$x_k$处取1呢？
很简单啊，把$x_k$代入刚刚的可能的因子，把代入的结果除掉不就行了？
也就是说$l_k(x)=\frac{\prod_{i \neq k}(x-x_i)}{\prod_{i \neq k}(x_k-x_i)}$
这就是拉格朗日插值基函数了。
则拉格朗日插值多项式就是
$$p(x)=\sum_{k=1}^n y_k l_k(x)$$
这就是拉格朗日插值公式了。
把拉格朗日插值进一步优化，你就能得到埃特金插值和牛顿插值。
这个就不多介绍了，比较复杂而且插值公式是唯一的，你用任何方法得到的插值公式都是一样的。而且市面上大把的直接计算插值的软件和库，他们直接内置了其他的插值方法。
拉格朗日插值余项
下面简单给个拉格朗日插值余弦的结论，其定理的证明网上导出都是，我也没必要再赘述了。
关于拉格朗日插值的误差，有如下的拉格朗日插值余项定理，
设[a,b]上有插值节点$x_0,x_1,...,x_n$，f(x)在[a,b]上有连续的直到n+1阶导数，且希望$f(x_i)=y_i$，
那么就有当$x \in [a,b]$时，有
余项$R_n(x)=f(x)-p_n(x)=\frac{f{(n+1)}(\xi)}{(n+1)!}\prod_{k=0}n(x-x_k)$
其中$\xi \in [a,b]$
这就是拉格朗日插值的余项了。由于插值多项式是唯一的，所以其他的插值方法的余项也是这个。
数值积分方法
下面我们来介绍一下怎么用插值法来得到数值方法计算函数积分。
对于大部分函数，我们其实都是很难求其积分的，甚至很多函数例如$e{x2}$这样的函数压根就没有初等原函数。所以我们需要寻找数值方法来计算他们的积分。
一种可能的思路就是用刚刚的插值多项式来近似替代这个函数。而且我们是已知了差值余项的，余项的积分就是我们数值积分方法的误差。
而这些都是多项式啊，都很好计算的。
这种求积分的数值方法得到的积分公式我们都叫做差值型的求积公式。
前面我说这只是数值求积的一种思路，其最本真的思路其实是取函数的部分函数值，用他们的线性组合来近似积分。
也就是最一般化的求积公式是$\int_a^b f(x)dx=\sum_{k=0}^n A_k y_k$。
只需要知道$A_k$就知道怎么求积了。
而对于插值型求积公式，我们很容易就能得到$A_k$。
插值型求积公式就是用f(x)的插值多项式近似f(x)，从而有
$\int_a^bf(x)dx \approx \int_a^bp_n(x)$
再把拉格朗日插值公式代入得
$\int_a^bf(x)dx \approx \int_a^b \sum_{k=1}^n y_k l_k(x)dx$
$=\sum_{k=1}^n y_k\int_a^bl_k(x)dx$
而我们最一般的求积公式是$\int_a^b f(x)dx=\sum_{k=0}^n A_k y_k$，从而
$A_k=\int_a^b l_k(x)dx$
牛顿-科特斯公式
如果我们进一步要求，插值点的间距相等（步长为h），我们可以得到一个特殊的结论。
我们竟然将区间长度提出来了，
$\int_a^b f(x)dx=(b-a)\sum_{k=0}^n C_k y_k$
并且，$C_k$竟然还是常数。
我们称$C_k$为科特斯系数，可以直接查科特斯系数表得到，与具体的函数f(x)无关，只和你插值多项式的次数n有关!
这就是牛顿-科特斯公式。
其中n=1得到的求积公式被称为梯形公式，n=2的是辛普森公式。
对牛顿-科特斯公式的玩法有很多，例如复化得到复化梯形公式和复化辛普森公式，还有递推化得到龙贝格公式等等。
还有变步长使得精度尽可能高的高斯公式等等玩法。
数值微分方法
有了积分当基石我们就可以用数值方法求解微分方程。
我们以标准的微分方程初值问题为例介绍几个简单的方法。
对微分方程
$$\frac{dy}{dx}=f(x,y)$$
我们两边同时积分得到
$$\int_{x_n}{x_{n+1}}\frac{dy}{dx}=\int_{x_n}{x_{n+1}}f(x,y)dx$$
也就是
$$y_{n+1}-y_n=\int_{x_n}^{x_{n+1}}f(x,y)dx$$
从而
$$y_{n+1}=y_n+\int_{x_n}^{x_{n+1}}f(x,y)dx$$
而如果我们可以用数值方法求出这个积分，不就解出了这个微分方程？
例如，我们直接用积分区间左端点得到的矩形的面积来近似积分（这个叫做左矩形求积公式）
得到
$$y_{n+1}=y_n+\int_{x_n}^{x_{n+1}}f(x,y)=y_n+(x_{n+1}-x_n)f(x_n,y_n)$$
这就是欧拉格式。(注意，数值微分里我们称为格式而不是公式，这里没有打错字)
你可能会发现我这里给的欧拉格式和网上一般的不一样，其实你只要设$h=x_{n+1}-x_n$（当我设了h就默认等步长了）就能得到一般的欧拉格式了。
$$y_{n+1}=y_n+hf(x_n,y_n)$$
同样地，我们可以用积分区间右端点得到的矩形面积（右矩形公式）近似积分得到
$$y_{n+1}=y_n+\int_{x_n}^{x_{n+1}}f(x,y)=y_n+(x_{n+1}-x_n)f(x_{n+1},y_{n+1})$$
类似地，我们设$h=x_{n+1}-x_n$就有，
$$y_{n+1}=y_n+hf(x_{n+1},y_{n+1})$$
这就是隐式欧拉格式了。
那如果我用区间中点的矩形，还能得到两步欧拉格式：
对微分方程
$$\frac{dy}{dx}=f(x,y)$$
我们两边同时积分得到
$$\int_{x_n-1}{x_{n+1}}\frac{dy}{dx}=\int_{x_n-1}{x_{n+1}}f(x,y)dx$$
注意，这次是从$x_{n-1}$积到$x_n$，并且我们要求步长为h，则$x_n$是积分区间中点。
有$y_{n+1}-y_{n-1}=2hf(x_n,y_n)$
从而，$y_{n+1}=y_{n-1}+2hf(x_n,y_n)$
这就是两步欧拉格式。
还可以用刚刚说的梯形公式来求积分得到梯形格式，或者引入预报-矫正系统得到改进的欧拉格式。
以及另外一个思路可以得到一整套的龙格-库塔方法和亚当姆斯方法。
