精华内容
下载资源
问答
  • 该源代码很简单的实现了拉格朗日求未知点的问题
  • 全域多项式插值指的是在整个插值区域内形成个多项式函数作为插值函数。关于多项式插值的基本知识,见“计算基本理论”。在单项式基插值和牛顿插值形成的表达式中,求该表达式在某一点处的值使用的Horner嵌套算法啊...

    全域多项式插值指的是在整个插值区域内形成一个多项式函数作为插值函数。关于多项式插值的基本知识,见“计算基本理论”。

    在单项式基插值和牛顿插值形成的表达式中,求该表达式在某一点处的值使用的Horner嵌套算法啊,见"Horner嵌套算法"。

    1. 单项式(Monomial)基插值

    1)插值函数基  单项式基插值采用的函数基是最简单的单项式:$$\phi_j(t)=t^{j-1}, j=1,2,...n;\quad f(t)=p_{n-1}(t)=x_1+x_2t+x_3t^2+...x_nt^{n-1}=\sum\limits_{j=1}^nx_jt^{j-1}$$  所要求解的系数即为单项式系数 $x_1,x_2,...x_n$ ,在这里仍然采用1,2,...n的下标记号而不采用和单项式指数对应的0,1,2,...,n-1的下标仅仅是出于和前后讨论一致的需要。

    2)叠加系数

    单项式基插值采用单项式函数基,若有m个离散数据点需要插值,设使用n项单项式基底:

    $$x_1+t_1x_2+t_1^2x_3+...+t_1^{n-1}x_n=y_1\\ x_1+t_2x_2+t_2^2x_3+...+t_2^{n-1}x_n=y_2\\ ......   ......   ......   ......   ......   ......\\ x_1+t_mx_2+t_m^2x_3+...+t_m^{n-1}x_n=y_m$$  系数矩阵为一 $m\times n$ 的矩阵($m\leq n$),范德蒙(Vandermonde)矩阵:

    $$\begin{bmatrix}1&t_1&t_1^2&...&t_1^{n-1}\\1&t_2&t_2^2&...&t_2^{n-1}\\...&...&...&...&...\\1&t_n&t_n^2&...&t_n^{n-1}\end{bmatrix} \begin{bmatrix}x_1\\x_2\\...\\x_n\end{bmatrix}=\begin{bmatrix}y_1\\y_2\\...\\y_n\end{bmatrix}$$  根据计算基本理论中的讨论,多项式插值的函数基一定线性无关,且只要离散数据点两两不同,所构成的矩阵行也一定线性无关,这保证了矩阵一定行满秩。此时当且仅当m=n时,叠加系数有且仅有一组解。因此,n=插值基底的个数=离散数据点的个数=多项式次数+1。

    3)问题条件和算法复杂度

    生成范德蒙矩阵的复杂度大致在 $O(n^2)$ 量级;由于范德蒙矩阵并没有什么好的性质,既没有稀疏性,也没有对称性,因此只能使用标准的稠密矩阵分解(如LU)来解决,复杂度在 $O(n^3)$ 量级。因此,问题的复杂度在 $O(n^3)$ 量级。

    范德蒙矩阵存在的问题是,当矩阵的维数越来越高的时候,解线性方程组时问题将越来越病态,条件数越来越大;从另一个角度来说,单项式基底本身趋势相近,次数增大时将越来越趋于平行(见下图)。这将造成随着数据点的增加,确定的叠加系数的不确定度越来越大,因此虽然单项式基很简单,进行插值时却往往不用这一方法。如果仍然采用单项式基底,有时也会进行两种可以改善以上问题的操作:平移(shifting)和缩放(scaling),即将 $((t-c)/d)^{j-1}$ 作为基底。常见的平移和缩放将所有数据点通过线性变换转移到区间[-1,1]中,即:$c=(t_1+t_n)/2,d=(t_n-t_1)/2$ 。

    20181006114053046299.png

    4)算法实现

    使用MATLAB实现单项式插值代码如下:

    function [ polyfunc, vecx, condition ] = MonoPI( vect, vecy, shift, scale )

    % 计算单项式型插值多项式系数

    % 输入四个参数:插值点行向量vect,插值点函数值行向量vecy,平移shift,压限scale;

    % 输出两个参数:插值多项式各项系数行向量vecx,矩阵条件数condition;

    % 设置缺省值:若只输入两个参数,则不平移不缩放

    if nargin==2

    shift = 0; scale = 1;

    end

    % 求解系数

    vecsize = size(vect, 2);

    basis = (vect - shift * ones(1, vecsize))/scale; % 确定基底在各个数据点的取值向量basis

    Mat = vander(basis); condition = cond(Mat); % 用vander命令生成basis的范德蒙矩阵并求条件数

    [L, U] = lu(Mat); vecx = (U\(L\vecy.‘)).‘; vecx = fliplr(vecx); % 标准lu分解解矩阵方程

    % 生成句柄函数polyfunc

    syms t;

    monomial = (t - shift)/scale; vecsize = size(vecx, 2); funeval = vecx(vecsize);

    for i = vecsize:-1:2 % 生成函数的算法采用Horner算法提高效率

    funeval = vecx(i - 1) + monomial*funeval;

    end

    polyfunc = matlabFunction(funeval, ‘Vars‘, t);

    end

    比如对于点:$(-2,-27),(0,-1),(1,0)$ 它具有唯一的二次插值多项式:$p_2(t)=-1+5t-4t^2$ 。调用以上代码:

    % 命令行输入

    [polyfunc, vecx, condition] = MonoPI(vect, vecy)

    % 命令行输出

    polyfunc =

    包含以下值的 function_handle:

    @(t)-t.*(t.*4.0-5.0)-1.0

    vecx =

    -1 5 -4

    condition =

    6.0809

    和预期完全一致。

    2. 拉格朗日(Lagrange)插值

    1)插值函数基

    拉格朗日插值采用的是一种设计巧妙的多项式基,每个基底都是n-1次多项式,而每个基底函数当且仅当在第i个数据点处取1,在其余数据点均为零。这个多项式基是这样设计的:

    $$l_j(t)=\frac{(t-t_1)(t-t_2)...(t-t_{j-1})(t-t_{j+1})...(t-t_n)}{(t_j-t_1)(t_j-t_2)...(t_j-t_{j-1})(t_j-t_{j+1})...(t_j-t_n)}=\frac{\prod\limits_{k=1,k\neq j}^n(t-t_k)}{\prod\limits_{k=1, k\neq j}^n(t_j-t_k)}$$  因此就有:

    $$l_j(t_i)=\delta_{ij}, i,j=1,2,...n $$  其中,$\delta$ 为克罗内克(Kronecker)记号,当两个下标相等时为1,否则为零;也可以将 $\delta_{ij}$ 理解为一个二阶张量,即单位矩阵。只要将各个$t_i$ 带入定义式,上式是很容易验证的。这意味着拉格朗日插值的叠加系数的求解将会产生很好的性质,即:

    2)叠加系数

    需要求解的插值函数即:$f(t)=\sum\limits_{k=1}^nx_kl_k(t)$ ,而又已知:

    $$l_1(t_1)x_1+l_2(t_1)x_2+...+l_n(t_1)x_n=y_1$$

    $$l_1(t_2)x_1+l_2(t_2)x_2+...+l_n(t_2)x_n=y_2$$

    $$... ... ... ... ...$$

    $$l_1(t_n)x_1+l_2(t_n)x_2+...+l_n(t_n)x_n=y_n$$  写成矩阵形式就是:

    $$\begin{bmatrix}l_1(t_1)&l_2(t_1)&...&l_n(t_1)\\l_1(t_2)&l_2(t_2)&...&l_n(t_2)\\...&...&...&...\\l_1(t_n)&l_2(t_n)&...&l_n(t_n)\end{bmatrix} \begin{bmatrix}x_1\\x_2\\...\\x_n\end{bmatrix}=\begin{bmatrix}1&0&..&0\\0&1&..&0\\..&..&..&..\\0&0&..&1\end{bmatrix} \begin{bmatrix}x_1\\x_2\\...\\x_n\end{bmatrix}=I\boldsymbol{x}=\begin{bmatrix}x_1\\x_2\\...\\x_n\end{bmatrix}=\begin{bmatrix}y_1\\y_2\\...\\y_n\end{bmatrix}$$

    其矩阵即单位矩阵,因此直接得出 $x_i=y_i$ ,$f(t)=p_{n-1}(t)=y_1l_1(t)+y_2l_2(t)+...+y_nl_n(t)=\sum\limits_{k=1}^ny_kl_k(t)$

    3)问题条件和算法复杂度

    拉格朗日插值生成的系数矩阵为单位矩阵,完全不存在条件病态的问题,只需要将各个数据点的取值作为系数即可。同样地,求解系数也将不存在任何复杂度。

    但是,作为代价的是,函数求值开销较大。Horner嵌套算法可以用于单项式和牛顿插值表达式的求值,将总运算量大致控制在n次浮点数加法和n次浮点数乘法,但该算法不适用于拉格朗日插值的表达式。拉格朗日插值的求值的复杂度至少也要3n次浮点乘(除)法和2n次浮点加法以上,这还是在所有的系数(将插值系数和基底的分母合并以后的系数)都已经处理完成之后,而处理系数本身可能就需要 $n^2$ 级别的复杂度。此外,拉格朗日插值表达式也不利于求微分和积分;和牛顿插值相比,当新增数据点时不得不将所有的基底都改写,很不方便。总而言之,拉格朗日插值属于非常容易构造的一种插值,很适用于在理论上讨论某些问题,但在数值计算上仍具有很多劣势。

    4)算法实现

    实现拉格朗日多项式插值一种途径的MATLAB代码如下。此处的输出为多项式函数句柄。这些函数(句柄)只需要在函数名后面加括号变量即可调用,即polyfun(3)这样的形式。

    function [ polyfun ] = LagrangePI( vect, vecy )

    % 生成Lagrange插值多项式

    % 输入两个参数:插值点行向量vect,函数行向量vecy;输出一个参数:插值多项式函数句柄polyfun

    vecsize = size(vect, 2);

    syms t f term;

    f = 0;

    for i = 1:vecsize

    term = vecy(i);

    for j = 1:vecsize

    if (j ~= i)

    term = term*(t - vect(j))/(vect(i) - vect(j));

    end

    end

    f = f + term;

    end

    polyfun = matlabFunction(f);

    end

    但是,由于多项式形式的函数表达式带入后为符号型变量,这意味着每一项的系数都经历了单独计算,每一项的分子也需要单独计算,这将使得拉格朗日插值表达式的函数求值(function evaluation)的复杂度达到 $O(n^2)$ 量级;如果想要使得每次求值能够控制在 $O(n)$ 量级,就必须实现计算出除了含有未知量的函数基分子以外的全部系数,同时在求值时也需要一些技巧。按照如下的书写方法可以实现这一目的:

    function [ coefficients ] = newLagrangePI( vect, vecy )

    % 生成Lagrange插值多项式的系数(计算分母)

    % 输入两个参数:插值点行向量vect,函数行向量vecy;

    % 输出一个参数:插值多项式的系数行向量coefficients;

    vecsize = size(vect, 2);

    coefficients = zeros(1, vecsize);

    for i = 1:vecsize

    tmp = vecy(i); % 取得基底函数对应的系数y_i

    for j = 1:vecsize % 将其除以函数基底的分母

    if (j~=i)

    tmp = tmp/(vect(i) - vect(j));

    end

    end

    coefficients(i) = tmp;

    end

    end

    除了求系数的函数还需要一个特别的求值函数:

    function [ funeval ] = evaLagrangePI( coefficients, vect, vecy, t )

    % Lagrange插值多项式估值

    % 输入四个参数:Lagrange插值的完整系数行向量coefficients,插值点行向量vect,函数行向量vecy,求值点t;

    % 输出一个参数:函数在t处取值funeval

    vecsize = size(vect, 2);

    [found, index] = ismember(t, vect);

    if found % 如果求值点是原数据点,则直接返回原始信息中数据点的函数值

    funeval = vecy(index);

    else % 否则,先计算全部(t-t_i)的乘积

    funeval = 0; product = 1;

    for i = 1:vecsize

    product = product*(t - vect(i));

    end

    for i = 1:vecsize % 然后,计算每一项的值,乘以该项的系数并且除以该项分子不存在的那项(t-t_i)

    funeval = funeval + coefficients(i)*product/(t - vect(i));

    end

    end

    end

    同样是对于三点 $(-2,-27),(0,-1),(1,0)$ ,调用Lagrange插值方法:

    vect = [-2, 0, 1]; vecy = [-27, -1, 0];

    % 命令行输入

    coefficients = newLagrangePI(vect, vecy)

    % 命令行输出

    coefficients =

    -4.5000 0.5000 0

    % 命令行输入

    val = evaLagrangePI(coefficients, vect, vecy, -2)

    % 命令行输出

    val =

    -27

    % 命令行输入

    val = evaLagrangePI(coefficients, vect, vecy, 0.5)

    % 命令行输出

    val =

    0.5000

    所有的输出均和实际的多项式插值 $f(t)=p_2(t)=-1+5t-4t^2$ 吻合。

    3. 牛顿(Newton)插值

    1)插值函数基底

    单项式基底非常简洁,缺点是求解方程组所用的是稠密的范德蒙矩阵,可能非常病态,复杂度也很高;拉格朗日基底比较精巧复杂,因为求解的系数矩阵是单位矩阵,求解很简单很准确,缺点是生成表达式和函数求值复杂度很高。牛顿插值方法在二者之间提供了一个折衷选项:基底不如拉格朗日的函数基那么复杂,而求解又比单项式基底大大简化,这来源于牛顿插值选取的基底:$$\pi_j(t)=\prod\limits_{k=1}^{j-1}(t-t_k)=(t-t_1)(t-t_2)...(t-t_{j-1}), j=1,...,n$$  相对于拉格朗日基底的特殊性($l_j(t_i)=\delta_{ij}$),牛顿插值基底具有一个弱一点的性质:$$\pi_j(t_i)=0,\forall i

    2)叠加系数

    $$\pi_1(t_1)x_1+\pi_2(t_1)x_2+...+\pi_n(t_1)x_n=y_1$$

    $$\pi_1(t_2)x_1+\pi_2(t_2)x_2+...+\pi_n(t_2)x_n=y_2$$

    $$............$$

    $$\pi_1(t_n)x_1+\pi_2(t_n)x_2+...+\pi_n(t_n)x_n=y_n$$  写成矩阵形式:

    $$\begin{bmatrix}\pi_1(t_1)&\pi_2(t_1)&...&\pi_n(t_1)\\ \pi_1(t_2)&\pi_2(t_2)&...&\pi_n(t_2)\\...&...&...&...\\ \pi_1(t_n)&\pi_2(t_n)&...&\pi_n(t_n)\end{bmatrix} \begin{bmatrix}x_1\\x_2\\...\\x_n\end{bmatrix}=\begin{bmatrix}1&0&...&0\\1&t_2-t_1&...&0\\...&...&...&...\\1&t_n-t_1&...&(t_n-t_1)..(t_n-t_{n-1})\end{bmatrix}=\begin{bmatrix}y_1\\y_2\\...\\y_n\end{bmatrix}$$  也就是说,牛顿插值的系数求解矩阵为一个下三角矩阵。

    3)算法性质和算法复杂度

    我们知道对于下三角矩阵,利用向前代入法可以比较便利地解出,其时间复杂度为 $O(n^2)$ 。再来看生成这个下三角矩阵,复杂度也是 $O(n^2)$ 的运算量。因此求解插值系数的总复杂度即 $O(n^2)$ 量级,比稠密矩阵的求解少一个量级。当然,求解牛顿插值系数不一定需要显示地生成矩阵,然后采用矩阵分解的标准套路;牛顿插值有好几种生成系数的方法可供选择,包括差商方法等,采用递归或者迭代都可以获得良好的效果,还能避免上溢出。

    此外,牛顿插值的表达式在求值时适用Horner嵌套算法(太棒了!),这将把求值的复杂度控制在 $O(n)$ 的量级内,如果带上系数比拉格朗日插值表达式的求值要更高效。

    牛顿插值方法有如下优越的性质:

    3.1 当增加数据点时,可以仅仅通过添加一项函数基和增加一个系数,生成新的牛顿插值多项式。这其实是可以理解的,因为当新增点 $(t_{n+1},y_{n+1})$ 时,下三角系数矩阵所有的元素都没有发生任何变化,仅仅需要新增一列和一行即可,而在该矩阵右侧新增的一列全为零。这意味着所有的系数 $x_1,x_2,...x_n$ 不仅满足原线性方程组,也因此必定是新线性方程组解的部分。而基底部分也只是新增了一个基,所以新的插值多项式仅仅是老的插值多项式加一项而已,即 $f(t)^*=p_n(t)=p_{n-1}(t)+x_{n+1}\pi_{n+1}(t)$ 。对于新的这一项 $x_{n+1}\pi_{n+1}(t)$ 它的系数究竟如何取,只需要将新增的数据点带入即可求得:$$f(t_{n+1})^*=p_{n-1}(t_{n+1})+x_{n+1}\pi_{n+1}(t_{n+1})=y_{n+1}\quad \Rightarrow \quad x_{n+1}=\frac{y_{n+1}-p_{n-1}(t_{n+1})}{\pi_{n+1}(t_{n+1})}$$  生成新系数的复杂度大致需要一次函数求值和一次基底求值,大致为 $O(n)$ 复杂度。如果迭代地使用这一公式,就可以生成全部牛顿插值多项式系数,复杂度 $O(n^2)$ ,和矩阵解法也大致是持平的。

    3.2 差商法也可以实现生成牛顿插值多项式的系数。其中,差商 $f[]$ 的定义为:

    $$f[t_1, t_2,...t_k]=\frac{f[t_2, t_3, ... t_{k}-f[t_1, t_2,...t_{k-1}]}{t_k-t_1}$$  而牛顿多项式的系数则取自:$$x_j=f[t_1, t_2... t_j]$$  对于这个可以有证明:

    $$f[t_1]=y_1, x_1=y_1=f[t_1];\quad f[t_1, t_2]=\frac{f[t_2]-f[t_1]}{t_2-t_1}=\frac{y_2-y_1}{t_2-t_1},x_2=\frac{y_2-y_1}{t_2-t_1}=f[t_1, t_2]

    $$

    若$x_j=f[t_1, t_2, ...t_j]=\frac{f[t_2, t_3,...t_j]-f[t_1, t_2,...t_{j-1}]}{t_j-t_1}$ 对于任意 $j\leq k-1$ 成立,当 $j=k$ 时:

    ?考虑对于点 $(t_1, y_1), (t_2,y_2),...(t_{k-1},y_{k-1})$ 的 Newton 插值多项式 $p_1(t)$ ;对于点 $(t_2, y_2),(t_3, y_3),...$$(t_k, y_k)$ 的 Newton 插值多项式 $p_2(t)$ ,并且已知分差插值系数对任意 $j\leq k-1$ 均成立,因而有:

    $$

    p_1(t)=\sum\limits_{j=1}^{k-1}a_j\pi_j(t), \quad p_2(t)=\sum\limits_{j=2}^{k}b_j\pi_j(t),\qquad a_j=f[t_1,...t_{j}],b_j=f[t_2,...t_j]

    $$

    由 $p_1(t)$ 过点 $(t_1, y_1)$ 到 $(t_{k-1},y_{k-1})$ ,$p_2(t)$ 过点 $(t_2,y_2)$ 到 $(t_k,y_k)$ ,构造插值多项式:

    $$

    p(t)=\frac{t_k-t}{t_k-t_1}p_1(t)+\frac{t-t_1}{t_k-t_1}p_2(t)

    $$

    就有该多项式通过点 $(t_1, y_1)$ 到 $(t_k,y_k)$ ,因此即为所求的 Newton 插值多项式。带入 $p_1(t),p_2(t)$ 表达式,并比较等式两端最高次项系数即得:

    $$

    p(t)=\sum\limits_{j=1}^kx_j\pi_j(t)=\frac{t_k-t}{t_k-t_1}\sum\limits_{j=1}^{k-1}a_j\pi_j(t)+\frac{t-t_1}{t_k-t_1}\sum\limits_{j=2}^{k}b_j\pi_j‘(t)\\

    x_k=\frac{-1}{t_k-t_1}a_{k-1}+\frac{1}{t_k-t_1}b_k=\frac{f[t_2,...t_k]-f[t_1,...t_{k-1}]}{t_k-t_1}=f[t_1, ...t_k]\qquad \square

    $$

    这个证明我摘录自奥斯陆大学数学系的 Michael S. Floater 在 Newton Interpolation 讲义里面写的证明。

    4)算法实现

    根据3.1,可以通过新增节点的方法迭代地生成插值系数。利用这种思路的实现代码如下:

    function [ vecx_new, vect_new ] = newNPI( vect, vecx, newPoint )

    %Newton插值算法新增节点函数;

    % 输入三个参数:原插值点行向量vect,原插值系数行向量vecx,新增节点newPoint;

    % 输入两个参数:新系数行向量vecx_new,新插值点行向量vect_new;

    vecsize = size(vecx, 2);

    vecx_new = zeros(1, vecsize + 1); vecx_new(1:vecsize) = vecx;

    vect_new = zeros(1, vecsize + 1); vect_new(1:vecsize) = vect; vect_new(vecsize + 1) = newPoint(1);

    p_new = HornerNPI(vect, vecx, newPoint(1)); w_new = 1;

    for i = 1:vecsize

    w_new = w_new * (newPoint(1) - vect(i));

    end

    vecx_new(vecsize + 1) = (newPoint(2) - p_new) / w_new;

    end

    新增节点函数newNPI可以单独使用;同时也可以反复调用生成牛顿插值系数,如下:

    function [ polyfun, vecx ] = newNewtonPI( cvect, cvecy )

    % 使用新增节点函数逐渐更新产生Newton插值多项式系数;

    % 输入两个参数:插值点行向量cvect,插值系数行向量cvecx;

    % 输出两个参数:多项式函数句柄polyfun,系数行向量vecx;

    % 迭代生成系数行向量

    vect = cvect(1); vecx = cvecy(1);

    vecsize = size(cvect, 2);

    for i=2:vecsize

    [vecx, vect] = newNPI(vect, vecx, [cvect(i), cvecy(i)]);

    end

    % 采用Horner嵌套算法生成多项式函数句柄

    syms f t; f = vecx(vecsize);

    for i = vecsize-1:-1:1

    f = vecx(i) + (t - cvect(i)) * f;

    end

    polyfun = matlabFunction(f);

    end

    另一种方法是采用差商。以下是实现的代码。和之前的说法不同的是,本代码使用的并非递归,而是正向的类似函数值缓存的算法。

    function [ polyfun, vecx ] = recNewtonPI( vect, vecy )

    % 使用差商产生Newton插值多项式系数;

    % 输入两个参数:插值点行向量vect,函数取值cvecy;

    % 输出两个参数:多项式函数polyfun,系数行向量vecx;

    vecsize = size(vect, 2);

    Div = diag(vecy);

    % 差商生成系数行向量vecx

    for k = 1:vecsize-1

    for i = 1:vecsize-k

    Div(i, i+k) = (Div(i+1, i+k) - Div(i, i+k-1))/(vect(i+k) - vect(i));

    end

    end

    vecx = Div(1, :);

    % 生成多项式函数polyfun

    syms f t; f = vecx(vecsize);

    for i = vecsize-1:-1:1

    f = vecx(i) + (t - vect(i)) * f;

    end

    polyfun = matlabFunction(f);

    end

    但不论如何,产生的结果完全一致。用同样的例子:

    vect=[-2, 0, 1]; vecy=[-27, -1, 0];

    % 命令行输入1,调用新增节点方法

    [polyfun, vecx] = newNewtonPI(vect, vecy)

    % 命令行输入2,调用差商方法

    [polyfun, vecx] = recNewtonPI(vect, vecy)

    % 命令行输出1/2,完全相同

    polyfun =

    包含以下值的 function_handle:

    @(t)-(t.*4.0-1.3e1).*(t+2.0)-2.7e1

    vecx =

    -27 13 -4

    容易检验,该多项式函数正是原数据点的多项式插值函数。

    展开全文
  • 双三次插值也是比较古典的个插值方法了,我最近一直在搞图像的插值,搞完双线性插值,终于迈进双三次了。。。这里有个问题 ,望求指点!!! 他的基函数为: 当 |w| < 1; s = 1-2*w*w + W*W*w 当 1=&...
  • 基于径向基函数(RBF)的函数插值

    千次阅读 多人点赞 2020-07-07 01:16:19
    1. 函数插值 2. RBF函数插值 代码实现

    基于径向基函数的函数插值

    1. 函数插值

    函数插值问题: 用形式简单的插值函数 f ^ ( x ) \hat f(x) f^(x) 近似原函数

    ( 1 ) \qquad(1) (1) 设函数 y = f ( x ) y=f(x) y=f(x) 在某个区间上有定义,并且已知该区间上的一些数据点 { x i , y i } \{x_i,y_i\} {xi,yi} 严格满足 y i = f ( x i ) , i = 1 , ⋯   , N y_i=f(x_i),i=1,\cdots,N yi=f(xi),i=1,,N,这些数据点称为“控制节点”或“插值节点

    ( 2 ) \qquad(2) (2) 如果存在一个形式上比较简单(比如 n n n 次多项式)的函数 f ^ ( x ) \hat f(x) f^(x),使得 f ^ ( x i ) = y i , i = 1 , ⋯   , N \hat f(x_i)=y_i,i=1,\cdots,N f^(xi)=yi,i=1,,N 都成立,就称 f ^ ( x ) \hat f(x) f^(x) f ( x ) f(x) f(x) 的插值函数。

    \qquad 典型的函数插值方法:拉格朗日插值和牛顿插值 H e r m i t e Hermite Hermite插值、样条插值等。

    与“函数逼近”的主要区别:
    插值函数 f ^ ( x ) \hat f(x) f^(x) 必须经过“插值节点”,也就是要满足 f ^ ( x i ) = y i , i = 1 , ⋯   , N \hat f(x_i)=y_i,i=1,\cdots,N f^(xi)=yi,i=1,,N

    2. RBF函数插值

    \qquad 与拉格朗日插值之类的常规函数插值不同,基于核函数的函数插值“通过引入核函数”来刻画数据的局部化特征

    \qquad 径向基函数 ( R a d i a l   B a s i s   F u n c t i o n , R B F ) (Radial\ Basis\ Function,RBF) (Radial Basis Function,RBF) 就是一类特殊的基函数,最常用的就是“高斯基函数”,定义为:

    φ ( x ) = e − x 2 2 σ 2 \qquad\qquad\qquad\varphi(x)=e^{-\frac{x^2}{2\sigma^2}} φ(x)=e2σ2x2  (以一维情况为例)

    \qquad 在这里插入图片描述
    RBF函数插值:  f ^ ( x ) = ∑ i = 1 N w i φ ( ∥ x − x i ∥ ) \hat f(x)=\displaystyle\sum_{i=1}^Nw_i\varphi(\parallel x-x_i\parallel) f^(x)=i=1Nwiφ(xxi)

    \qquad 假设有 N N N 个插值节点,也就是已知 { x j , y j } ∣ j = 1 N \{x_j,y_j\}\big |_{j=1}^N {xj,yj}j=1N,其中 f ^ ( x j ) = y j = f ( x j ) \hat f(x_j)=y_j=f(x_j) f^(xj)=yj=f(xj),如下图所示。
    \qquad 在这里插入图片描述

    图中,红色实线为真实函数曲线,绿色空心圆圈代表插值节点 ( x j , y j ) (x_j,y_j) (xj,yj)蓝色实心点为RBF插值所求得的权值 w j w_j wj

    \qquad { x j , y j } ∣ j = 1 N \{x_j,y_j\}\big |_{j=1}^N {xj,yj}j=1N 带入方程 f ^ ( x ) = ∑ i = 1 N w i φ ( ∥ x − x i ∥ ) \hat f(x)=\displaystyle\sum_{i=1}^Nw_i\varphi(\parallel x-x_i\parallel) f^(x)=i=1Nwiφ(xxi),可得到:

    [ φ 11 φ 12 ⋯ φ 1 N φ 21 φ 22 ⋯ φ 2 N ⋮ ⋮ ⋮ φ 11 φ 12 ⋯ φ 1 N ] ⏟ Φ [ w 1 w 2 ⋮ w N ] ⏟ W = [ y 1 y 2 ⋮ y N ] ⏟ y \qquad\qquad\underbrace{\left[ \begin{matrix} \varphi_{11} & \varphi_{12} & \cdots & \varphi_{1N} \\ \varphi_{21} & \varphi_{22} & \cdots & \varphi_{2N} \\ \vdots & \vdots & &\vdots \\ \varphi_{11} & \varphi_{12} & \cdots & \varphi_{1N} \end{matrix} \right] }_{\Phi}\underbrace{ \left[ \begin{matrix} w_1 \\ w_2 \\ \vdots \\ w_N \end{matrix} \right]}_{\bold W}=\underbrace{\left[ \begin{matrix} y_1 \\ y_2 \\ \vdots \\ y_N \end{matrix} \right]}_{\bold y} Φ φ11φ21φ11φ12φ22φ12φ1Nφ2Nφ1NW w1w2wN=y y1y2yN,其中 φ j i = φ ( ∥ x j − x i ∥ ) \varphi_{ji}=\varphi(\parallel x_j-x_i\parallel) φji=φ(xjxi)

    \qquad
    \qquad 其中, Φ = [ φ j i ] \Phi=[\varphi_{ji}] Φ=[φji]插值矩阵。因为 φ j i = φ ( ∥ x j − x i ∥ ) = φ i j \varphi_{ji}=\varphi(\parallel x_j-x_i\parallel)=\varphi_{ij} φji=φ(xjxi)=φij,因此插值矩阵是对称的。对于高斯核函数而言,插值矩阵的对角线元素的值为 1 1 1

    \qquad 将线性方程组记为 Φ W = y \Phi\bold W=\bold y ΦW=y,该方程组的第 j j j 行为:

    f ^ ( x j ) = y j = w 1 φ ( ∥ x j − x 1 ∥ ) + w 2 φ ( ∥ x j − x 2 ∥ ) + ⋯ + w N φ ( ∥ x j − x N ∥ ) \qquad\qquad\hat f(x_j)=y_j=w_1\varphi(\parallel x_j-x_1\parallel)+w_2\varphi(\parallel x_j-x_2\parallel)+\cdots+w_N\varphi(\parallel x_j-x_N\parallel) f^(xj)=yj=w1φ(xjx1)+w2φ(xjx2)++wNφ(xjxN)

    \qquad 因此,可求出 R B F RBF RBF 插值的系数为: W = Φ − 1 y \bold W=\Phi^{-1}\bold y W=Φ1y,其示意图如下图所示。

    Micchelli定理可以保证采用高斯函数时,插值矩阵 Φ \Phi Φ 是可逆的(只要插值节点互不相同)。

    \qquad 在这里插入图片描述

    代码实现

    import numpy as np
    import matplotlib.pyplot as plt
    
    def gen_data(x1,x2):
        y_sample = np.sin(np.pi*x1/2)+np.cos(np.pi*x1/3)
        y_all = np.sin(np.pi*x2/2)+np.cos(np.pi*x2/3)
        return y_sample, y_all
    
    def kernel_interpolation(y_sample,x1,sig):
        gaussian_kernel = lambda x,c,h: np.exp(-(x-x[c])**2/(2*(h**2)))
        num = len(y_sample)
        w = np.zeros(num)
        int_matrix = np.asmatrix(np.zeros((num,num)))
        for i in range(num):
            int_matrix[i,:] = gaussian_kernel(x1,i,sig)  
        w = int_matrix.I * np.asmatrix(y_sample).T      
        return w
    
    def kernel_interpolation_rec(w,x1,x2,sig):
        gkernel = lambda x,xc,h: np.exp(-(x-xc)**2/(2*(h**2)))
        num = len(x2)
        y_rec = np.zeros(num)
        for i in range(num):
            for k in range(len(w)):
                y_rec[i] = y_rec[i] + w[k]*gkernel(x2[i],x1[k],sig)
        return y_rec
    
    if __name__ == '__main__':
        snum = 20   # control point数量
        ratio = 20  # 总数据点数量:snum*ratio
        sig = 1		# 核函数宽度
        xs = -8
        xe = 8
        x1 = np.linspace(xs,xe,snum)
        x2 = np.linspace(xs,xe,(snum-1)*ratio+1)
        y_sample, y_all = gen_data(x1,x2)
        plt.figure(1)
        w = kernel_interpolation(y_sample,x1,sig)   
        y_rec = kernel_interpolation_rec(w,x1,x2,sig)
        plt.plot(x2,y_rec,'k')
        plt.plot(x2,y_all,'r:')        
        plt.ylabel('y')
        plt.xlabel('x')
        for i in range(len(x1)):        
            plt.plot(x1[i],y_sample[i],'go',markerfacecolor='none')        
            
        plt.legend(labels=['reconstruction','original','control point'],loc='lower left')
        plt.title('kernel interpolation:$y=sin(\pi x/2)+cos(\pi x/3)$')       
        plt.show()
    

    运行结果:
    在这里插入图片描述

    在相同区间、分别采用 8 , 12 , 16 , 20 8,12,16,20 8,12,16,20 个控制节点 ( c o n t r o l   p o i n t ) (control\ point) (control point) 进行函数插值的结果
    显然,插值节点过少,无法体现整个函数的特征;插值节点越多,函数插值的结果越精确

    在这里插入图片描述

    扩大插值区间范围,控制节点 ( c o n t r o l   p o i n t ) (control\ point) (control point) 也需要增加数量,才能保持函数插值的准确性

    \quad
    另外,Scipy的插值模块也提供了RBF插值,其实现代码如下:

    import numpy as np
    import matplotlib.pyplot as plt
    from scipy.interpolate import Rbf
    
    f = lambda x: np.sin(np.pi*x/2)+np.cos(np.pi*x/3)
    snum = 20   # control point数量
    ratio = 20  # 总数据点数量:snum*ratio
    xs = -8
    xe = 8
    x1 = np.linspace(xs,xe,snum)				# control points
    x2 = np.linspace(xs,xe,(snum-1)*ratio+1)	# 作图总数据点
    y1 = f(x1)			# control points
    rf = Rbf(x1, y1)	# reconstructed Rbf function
    y2 = rf(x2)  		# Rbf reconstruction
    plt.plot(x2, y2, 'k-', x2, f(x2),'r-', x1, y1, 'go', markerfacecolor='none')
    plt.legend(["Radial basis functions", "Orignal", "control point"],loc='best')
    plt.show()
    

    其运行结果如下:
    在这里插入图片描述

    \quad
    此外,RBF函数插值还可以通过径向基函数网络来实现。

    \quad
    参考:函数插值的python实现——拉格朗日、牛顿插值

    展开全文
  • 通过计算行列式的值,对几种Hermite插值多项式的存在唯一性给出另种证明方法,对带不完全导数的m(m≥4)Hermite插值多项式,给出推广的基函数构造方法,并对带不完全导数的三及四Hermite插值多项式的具体实例,给...
  • 通过引进三角域上的插值基函数,给出了种新的三角域上的二元三次插值样条函数,这种插值样条函数整体达到C’连续,且在各网格点处的参数可由递推公式得到。文中给出的插值样条函数较之Farin提出的分裂三角形方法,...
  • 全域多项式插值指的是在整个插值区域内形成个多项式函数作为插值函数。关于多项式插值的基本知识,见“计算基本理论”。  在单项式基插值和牛顿插值形成的表达式中,求该表达式在某一点处的值使用的Horner嵌套...

      全域多项式插值指的是在整个插值区域内形成一个多项式函数作为插值函数。关于多项式插值的基本知识,见“计算基本理论”

      在单项式基插值和牛顿插值形成的表达式中,求该表达式在某一点处的值使用的Horner嵌套算法啊,见"Horner嵌套算法"。

    1. 单项式(Monomial)基插值

    1)插值函数基  单项式基插值采用的函数基是最简单的单项式:$$\phi_j(t)=t^{j-1}, j=1,2,...n;\quad f(t)=p_{n-1}(t)=x_1+x_2t+x_3t^2+...x_nt^{n-1}=\sum\limits_{j=1}^nx_jt^{j-1}$$  所要求解的系数即为单项式系数 $x_1,x_2,...x_n$ ,在这里仍然采用1,2,...n的下标记号而不采用和单项式指数对应的0,1,2,...,n-1的下标仅仅是出于和前后讨论一致的需要。

    2)叠加系数

      单项式基插值采用单项式函数基,若有m个离散数据点需要插值,设使用n项单项式基底:

    $$x_1+t_1x_2+t_1^2x_3+...+t_1^{n-1}x_n=y_1\\ x_1+t_2x_2+t_2^2x_3+...+t_2^{n-1}x_n=y_2\\ ......   ......   ......   ......   ......   ......\\ x_1+t_mx_2+t_m^2x_3+...+t_m^{n-1}x_n=y_m$$  系数矩阵为一 $m\times n$ 的矩阵($m\leq n$),范德蒙(Vandermonde)矩阵

    $$\begin{bmatrix}1&t_1&t_1^2&...&t_1^{n-1}\\1&t_2&t_2^2&...&t_2^{n-1}\\...&...&...&...&...\\1&t_n&t_n^2&...&t_n^{n-1}\end{bmatrix} \begin{bmatrix}x_1\\x_2\\...\\x_n\end{bmatrix}=\begin{bmatrix}y_1\\y_2\\...\\y_n\end{bmatrix}$$  根据计算基本理论中的讨论,多项式插值的函数基一定线性无关,且只要离散数据点两两不同,所构成的矩阵行也一定线性无关,这保证了矩阵一定行满秩。此时当且仅当m=n时,叠加系数有且仅有一组解。因此,n=插值基底的个数=离散数据点的个数=多项式次数+1。

    3)问题条件和算法复杂度

      生成范德蒙矩阵的复杂度大致在 $O(n^2)$ 量级;由于范德蒙矩阵并没有什么好的性质,既没有稀疏性,也没有对称性,因此只能使用标准的稠密矩阵分解(如LU)来解决,复杂度在 $O(n^3)$ 量级。因此,问题的复杂度在 $O(n^3)$ 量级。 

      范德蒙矩阵存在的问题是,当矩阵的维数越来越高的时候,解线性方程组时问题将越来越病态,条件数越来越大;从另一个角度来说,单项式基底本身趋势相近,次数增大时将越来越趋于平行(见下图)。这将造成随着数据点的增加,确定的叠加系数的不确定度越来越大,因此虽然单项式基很简单,进行插值时却往往不用这一方法。如果仍然采用单项式基底,有时也会进行两种可以改善以上问题的操作:平移(shifting)和缩放(scaling),即将 $((t-c)/d)^{j-1}$ 作为基底。常见的平移和缩放将所有数据点通过线性变换转移到区间[-1,1]中,即:$c=(t_1+t_n)/2,d=(t_n-t_1)/2$ 。

    4)算法实现 

      使用MATLAB实现单项式插值代码如下:

    function [ polyfunc, vecx, condition ] = MonoPI( vect, vecy, shift, scale )
    %    计算单项式型插值多项式系数
    %    输入四个参数:插值点行向量vect,插值点函数值行向量vecy,平移shift,压限scale;
    %    输出两个参数:插值多项式各项系数行向量vecx,矩阵条件数condition;
    
    % 设置缺省值:若只输入两个参数,则不平移不缩放
    if nargin==2
        shift = 0; scale = 1;
    end
    
    % 求解系数
    vecsize = size(vect, 2);
    basis = (vect - shift * ones(1, vecsize))/scale;    % 确定基底在各个数据点的取值向量basis
    Mat = vander(basis); condition = cond(Mat);    % 用vander命令生成basis的范德蒙矩阵并求条件数
    [L, U] = lu(Mat); vecx = (U\(L\vecy.')).'; vecx = fliplr(vecx);    % 标准lu分解解矩阵方程
    
    % 生成句柄函数polyfunc
    syms t;
    monomial = (t - shift)/scale; vecsize = size(vecx, 2); funeval = vecx(vecsize);
    for i = vecsize:-1:2    % 生成函数的算法采用Horner算法提高效率
        funeval = vecx(i - 1) + monomial*funeval;
    end
    polyfunc = matlabFunction(funeval, 'Vars', t);
    
    end
    

      比如对于点:$(-2,-27),(0,-1),(1,0)$ 它具有唯一的二次插值多项式:$p_2(t)=-1+5t-4t^2$ 。调用以上代码:

    % 命令行输入
    [polyfunc, vecx, condition] = MonoPI(vect, vecy)
    % 命令行输出
    polyfunc =
      包含以下值的 function_handle:
        @(t)-t.*(t.*4.0-5.0)-1.0
    vecx =
        -1     5    -4
    condition =
        6.0809
    

      和预期完全一致。

     

    2. 拉格朗日(Lagrange)插值

    1)插值函数基

      拉格朗日插值采用的是一种设计巧妙的多项式基,每个基底都是n-1次多项式,而每个基底函数当且仅当在第i个数据点处取1,在其余数据点均为零。这个多项式基是这样设计的:

    $$l_j(t)=\frac{(t-t_1)(t-t_2)...(t-t_{j-1})(t-t_{j+1})...(t-t_n)}{(t_j-t_1)(t_j-t_2)...(t_j-t_{j-1})(t_j-t_{j+1})...(t_j-t_n)}=\frac{\prod\limits_{k=1,k\neq j}^n(t-t_k)}{\prod\limits_{k=1, k\neq j}^n(t_j-t_k)}$$  因此就有:

    $$l_j(t_i)=\delta_{ij}, i,j=1,2,...n $$  其中,$\delta$ 为克罗内克(Kronecker)记号,当两个下标相等时为1,否则为零;也可以将 $\delta_{ij}$ 理解为一个二阶张量,即单位矩阵。只要将各个$t_i$ 带入定义式,上式是很容易验证的。这意味着拉格朗日插值的叠加系数的求解将会产生很好的性质,即:

    2)叠加系数

      需要求解的插值函数即:$f(t)=\sum\limits_{k=1}^nx_kl_k(t)$ ,而又已知:

    $$l_1(t_1)x_1+l_2(t_1)x_2+...+l_n(t_1)x_n=y_1$$

    $$l_1(t_2)x_1+l_2(t_2)x_2+...+l_n(t_2)x_n=y_2$$

    $$... ... ... ... ...$$

    $$l_1(t_n)x_1+l_2(t_n)x_2+...+l_n(t_n)x_n=y_n$$  写成矩阵形式就是:

    $$\begin{bmatrix}l_1(t_1)&l_2(t_1)&...&l_n(t_1)\\l_1(t_2)&l_2(t_2)&...&l_n(t_2)\\...&...&...&...\\l_1(t_n)&l_2(t_n)&...&l_n(t_n)\end{bmatrix} \begin{bmatrix}x_1\\x_2\\...\\x_n\end{bmatrix}=\begin{bmatrix}1&0&..&0\\0&1&..&0\\..&..&..&..\\0&0&..&1\end{bmatrix} \begin{bmatrix}x_1\\x_2\\...\\x_n\end{bmatrix}=I\boldsymbol{x}=\begin{bmatrix}x_1\\x_2\\...\\x_n\end{bmatrix}=\begin{bmatrix}y_1\\y_2\\...\\y_n\end{bmatrix}$$

      其矩阵即单位矩阵,因此直接得出 $x_i=y_i$ ,$f(t)=p_{n-1}(t)=y_1l_1(t)+y_2l_2(t)+...+y_nl_n(t)=\sum\limits_{k=1}^ny_kl_k(t)$

    3)问题条件和算法复杂度

      拉格朗日插值生成的系数矩阵为单位矩阵,完全不存在条件病态的问题,只需要将各个数据点的取值作为系数即可。同样地,求解系数也将不存在任何复杂度

      但是,作为代价的是,函数求值开销较大。Horner嵌套算法可以用于单项式和牛顿插值表达式的求值,将总运算量大致控制在n次浮点数加法和n次浮点数乘法,但该算法不适用于拉格朗日插值的表达式。拉格朗日插值的求值的复杂度至少也要3n次浮点乘(除)法和2n次浮点加法以上,这还是在所有的系数(将插值系数和基底的分母合并以后的系数)都已经处理完成之后,而处理系数本身可能就需要 $n^2$ 级别的复杂度。此外,拉格朗日插值表达式也不利于求微分和积分;和牛顿插值相比,当新增数据点时不得不将所有的基底都改写,很不方便。总而言之,拉格朗日插值属于非常容易构造的一种插值,很适用于在理论上讨论某些问题,但在数值计算上仍具有很多劣势。

    4)算法实现

       实现拉格朗日多项式插值一种途径的MATLAB代码如下。此处的输出为多项式函数句柄。这些函数(句柄)只需要在函数名后面加括号变量即可调用,即polyfun(3)这样的形式。

    function [ polyfun ] = LagrangePI( vect, vecy )
    %   生成Lagrange插值多项式
    %   输入两个参数:插值点行向量vect,函数行向量vecy;输出一个参数:插值多项式函数句柄polyfun
    vecsize = size(vect, 2);
    syms t f term;
    f = 0;
    for i = 1:vecsize
        term = vecy(i);
        for j = 1:vecsize
            if (j ~= i)
                term = term*(t - vect(j))/(vect(i) - vect(j));
            end
        end
        f = f + term;
    end
    polyfun = matlabFunction(f);
    end
    

      但是,由于多项式形式的函数表达式带入后为符号型变量,这意味着每一项的系数都经历了单独计算,每一项的分子也需要单独计算,这将使得拉格朗日插值表达式的函数求值(function evaluation)的复杂度达到 $O(n^2)$ 量级;如果想要使得每次求值能够控制在 $O(n)$ 量级,就必须实现计算出除了含有未知量的函数基分子以外的全部系数,同时在求值时也需要一些技巧。按照如下的书写方法可以实现这一目的:

    function [ coefficients ] = newLagrangePI( vect, vecy )
    %   生成Lagrange插值多项式的系数(计算分母)
    %   输入两个参数:插值点行向量vect,函数行向量vecy;
    %   输出一个参数:插值多项式的系数行向量coefficients;
    vecsize = size(vect, 2);
    coefficients = zeros(1, vecsize);
    for i = 1:vecsize
        tmp = vecy(i);    % 取得基底函数对应的系数y_i
        for j = 1:vecsize    % 将其除以函数基底的分母
            if (j~=i)
                tmp = tmp/(vect(i) - vect(j));
            end
        end
        coefficients(i) = tmp;
    end
    end
    

      除了求系数的函数还需要一个特别的求值函数:

    function [ funeval ] = evaLagrangePI( coefficients, vect, vecy, t )
    %   Lagrange插值多项式估值
    %   输入四个参数:Lagrange插值的完整系数行向量coefficients,插值点行向量vect,函数行向量vecy,求值点t;
    %   输出一个参数:函数在t处取值funeval
    vecsize = size(vect, 2);
    [found, index] = ismember(t, vect);
    if found    % 如果求值点是原数据点,则直接返回原始信息中数据点的函数值
        funeval = vecy(index);
    else    % 否则,先计算全部(t-t_i)的乘积
        funeval = 0; product = 1;
        for i = 1:vecsize
            product = product*(t - vect(i));
        end
        for i = 1:vecsize    % 然后,计算每一项的值,乘以该项的系数并且除以该项分子不存在的那项(t-t_i)
            funeval = funeval + coefficients(i)*product/(t - vect(i));
        end
    end
    end
    

      同样是对于三点 $(-2,-27),(0,-1),(1,0)$ ,调用Lagrange插值方法:

    vect = [-2, 0, 1]; vecy = [-27, -1, 0];
    % 命令行输入
    coefficients = newLagrangePI(vect, vecy)
    % 命令行输出
    coefficients =
       -4.5000    0.5000         0
    
    % 命令行输入
    val = evaLagrangePI(coefficients, vect, vecy, -2)
    % 命令行输出
    val =
       -27
    
    % 命令行输入
    val = evaLagrangePI(coefficients, vect, vecy, 0.5)
    % 命令行输出
    val =
        0.5000
    

      所有的输出均和实际的多项式插值 $f(t)=p_2(t)=-1+5t-4t^2$ 吻合。

     

    3. 牛顿(Newton)插值

    1)插值函数基底

      单项式基底非常简洁,缺点是求解方程组所用的是稠密的范德蒙矩阵,可能非常病态,复杂度也很高;拉格朗日基底比较精巧复杂,因为求解的系数矩阵是单位矩阵,求解很简单很准确,缺点是生成表达式和函数求值复杂度很高。牛顿插值方法在二者之间提供了一个折衷选项:基底不如拉格朗日的函数基那么复杂,而求解又比单项式基底大大简化,这来源于牛顿插值选取的基底:$$\pi_j(t)=\prod\limits_{k=1}^{j-1}(t-t_k)=(t-t_1)(t-t_2)...(t-t_{j-1}), j=1,...,n$$  相对于拉格朗日基底的特殊性($l_j(t_i)=\delta_{ij}$),牛顿插值基底具有一个弱一点的性质:$$\pi_j(t_i)=0,\forall i<j$$  求出的多项式形如:$f(t)=p_{n-1}(t)=\sum\limits_{j=1}^nx_j\pi_j(t)=x_1+x_2(t-t_1)+...+x_n(t-t_1)(t-t_2)...(t-t_{n-1})$

    2)叠加系数

    $$\pi_1(t_1)x_1+\pi_2(t_1)x_2+...+\pi_n(t_1)x_n=y_1$$

    $$\pi_1(t_2)x_1+\pi_2(t_2)x_2+...+\pi_n(t_2)x_n=y_2$$

    $$............$$

    $$\pi_1(t_n)x_1+\pi_2(t_n)x_2+...+\pi_n(t_n)x_n=y_n$$  写成矩阵形式:

    $$\begin{bmatrix}\pi_1(t_1)&\pi_2(t_1)&...&\pi_n(t_1)\\ \pi_1(t_2)&\pi_2(t_2)&...&\pi_n(t_2)\\...&...&...&...\\ \pi_1(t_n)&\pi_2(t_n)&...&\pi_n(t_n)\end{bmatrix} \begin{bmatrix}x_1\\x_2\\...\\x_n\end{bmatrix}=\begin{bmatrix}1&0&...&0\\1&t_2-t_1&...&0\\...&...&...&...\\1&t_n-t_1&...&(t_n-t_1)..(t_n-t_{n-1})\end{bmatrix}=\begin{bmatrix}y_1\\y_2\\...\\y_n\end{bmatrix}$$  也就是说,牛顿插值的系数求解矩阵为一个下三角矩阵

    3)算法性质和算法复杂度

      我们知道对于下三角矩阵,利用向前代入法可以比较便利地解出,其时间复杂度为 $O(n^2)$ 。再来看生成这个下三角矩阵,复杂度也是 $O(n^2)$ 的运算量。因此求解插值系数的总复杂度即 $O(n^2)$ 量级,比稠密矩阵的求解少一个量级。当然,求解牛顿插值系数不一定需要显示地生成矩阵,然后采用矩阵分解的标准套路;牛顿插值有好几种生成系数的方法可供选择,包括差商方法等,采用递归或者迭代都可以获得良好的效果,还能避免上溢出。

      此外,牛顿插值的表达式在求值时适用Horner嵌套算法(太棒了!),这将把求值的复杂度控制在 $O(n)$ 的量级内,如果带上系数比拉格朗日插值表达式的求值要更高效。

    牛顿插值方法有如下优越的性质:

    3.1 当增加数据点时,可以仅仅通过添加一项函数基和增加一个系数,生成新的牛顿插值多项式。这其实是可以理解的,因为当新增点 $(t_{n+1},y_{n+1})$ 时,下三角系数矩阵所有的元素都没有发生任何变化,仅仅需要新增一列和一行即可,而在该矩阵右侧新增的一列全为零。这意味着所有的系数 $x_1,x_2,...x_n$ 不仅满足原线性方程组,也因此必定是新线性方程组解的部分。而基底部分也只是新增了一个基,所以新的插值多项式仅仅是老的插值多项式加一项而已,即 $f(t)^*=p_n(t)=p_{n-1}(t)+x_{n+1}\pi_{n+1}(t)$ 。对于新的这一项 $x_{n+1}\pi_{n+1}(t)$ 它的系数究竟如何取,只需要将新增的数据点带入即可求得:$$f(t_{n+1})^*=p_{n-1}(t_{n+1})+x_{n+1}\pi_{n+1}(t_{n+1})=y_{n+1}\quad \Rightarrow \quad x_{n+1}=\frac{y_{n+1}-p_{n-1}(t_{n+1})}{\pi_{n+1}(t_{n+1})}$$  生成新系数的复杂度大致需要一次函数求值和一次基底求值,大致为 $O(n)$ 复杂度。如果迭代地使用这一公式,就可以生成全部牛顿插值多项式系数,复杂度 $O(n^2)$ ,和矩阵解法也大致是持平的。

    3.2 差商法也可以实现生成牛顿插值多项式的系数。其中,差商 $f[]$ 的定义为:

    $$f[t_1, t_2,...t_k]=\frac{f[t_2, t_3, ... t_{k}-f[t_1, t_2,...t_{k-1}]}{t_k-t_1}$$  而牛顿多项式的系数则取自:$$x_j=f[t_1, t_2... t_j]$$  对于这个可以有证明:


     

    $$f[t_1]=y_1, x_1=y_1=f[t_1];\quad f[t_1, t_2]=\frac{f[t_2]-f[t_1]}{t_2-t_1}=\frac{y_2-y_1}{t_2-t_1},x_2=\frac{y_2-y_1}{t_2-t_1}=f[t_1, t_2]
    $$
    若$x_j=f[t_1, t_2, ...t_j]=\frac{f[t_2, t_3,...t_j]-f[t_1, t_2,...t_{j-1}]}{t_j-t_1}$ 对于任意 $j\leq k-1$ 成立,当 $j=k$ 时:

    ​ 考虑对于点 $(t_1, y_1), (t_2,y_2),...(t_{k-1},y_{k-1})$ 的 Newton 插值多项式 $p_1(t)$ ;对于点 $(t_2, y_2),(t_3, y_3),...$$(t_k, y_k)$ 的 Newton 插值多项式 $p_2(t)$ ,并且已知分差插值系数对任意 $j\leq k-1$ 均成立,因而有:
    $$
    p_1(t)=\sum\limits_{j=1}^{k-1}a_j\pi_j(t), \quad p_2(t)=\sum\limits_{j=2}^{k}b_j\pi_j(t),\qquad a_j=f[t_1,...t_{j}],b_j=f[t_2,...t_j]
    $$
    由 $p_1(t)$ 过点 $(t_1, y_1)$ 到 $(t_{k-1},y_{k-1})$ ,$p_2(t)$ 过点 $(t_2,y_2)$ 到 $(t_k,y_k)$ ,构造插值多项式:
    $$
    p(t)=\frac{t_k-t}{t_k-t_1}p_1(t)+\frac{t-t_1}{t_k-t_1}p_2(t)
    $$
    就有该多项式通过点 $(t_1, y_1)$ 到 $(t_k,y_k)$ ,因此即为所求的 Newton 插值多项式。带入 $p_1(t),p_2(t)$ 表达式,并比较等式两端最高次项系数即得:
    $$
    p(t)=\sum\limits_{j=1}^kx_j\pi_j(t)=\frac{t_k-t}{t_k-t_1}\sum\limits_{j=1}^{k-1}a_j\pi_j(t)+\frac{t-t_1}{t_k-t_1}\sum\limits_{j=2}^{k}b_j\pi_j'(t)\\
    x_k=\frac{-1}{t_k-t_1}a_{k-1}+\frac{1}{t_k-t_1}b_k=\frac{f[t_2,...t_k]-f[t_1,...t_{k-1}]}{t_k-t_1}=f[t_1, ...t_k]\qquad \square
    $$

     


     这个证明我摘录自奥斯陆大学数学系的 Michael S. Floater 在 Newton Interpolation 讲义里面写的证明。

    4)算法实现

       根据3.1,可以通过新增节点的方法迭代地生成插值系数。利用这种思路的实现代码如下:

    function [ vecx_new, vect_new ] = newNPI( vect, vecx, newPoint )
    %	Newton插值算法新增节点函数;
    %   输入三个参数:原插值点行向量vect,原插值系数行向量vecx,新增节点newPoint;
    %   输入两个参数:新系数行向量vecx_new,新插值点行向量vect_new;
    vecsize = size(vecx, 2);
    vecx_new = zeros(1, vecsize + 1); vecx_new(1:vecsize) = vecx;
    vect_new = zeros(1, vecsize + 1); vect_new(1:vecsize) = vect; vect_new(vecsize + 1) = newPoint(1);
    p_new = HornerNPI(vect, vecx, newPoint(1)); w_new = 1;
    for i = 1:vecsize
        w_new = w_new * (newPoint(1) - vect(i));
    end
    vecx_new(vecsize + 1) = (newPoint(2) - p_new) / w_new;
    end
    

      新增节点函数newNPI可以单独使用;同时也可以反复调用生成牛顿插值系数,如下:

    function [ polyfun, vecx ] = newNewtonPI( cvect, cvecy )
    %    使用新增节点函数逐渐更新产生Newton插值多项式系数;
    %    输入两个参数:插值点行向量cvect,插值系数行向量cvecx;
    %    输出两个参数:多项式函数句柄polyfun,系数行向量vecx;
    
    %    迭代生成系数行向量
    vect = cvect(1); vecx = cvecy(1);
    vecsize = size(cvect, 2);
    for i=2:vecsize
        [vecx, vect] = newNPI(vect, vecx, [cvect(i), cvecy(i)]);
    end
    
    %    采用Horner嵌套算法生成多项式函数句柄
    syms f t; f = vecx(vecsize);
    for i = vecsize-1:-1:1
        f = vecx(i) + (t - cvect(i)) * f;
    end
    polyfun = matlabFunction(f);
    end
    

      另一种方法是采用差商。以下是实现的代码。和之前的说法不同的是,本代码使用的并非递归,而是正向的类似函数值缓存的算法。

    function [ polyfun, vecx ] = recNewtonPI( vect, vecy )
    %    使用差商产生Newton插值多项式系数;
    %    输入两个参数:插值点行向量vect,函数取值cvecy;
    %    输出两个参数:多项式函数polyfun,系数行向量vecx;
    vecsize = size(vect, 2);
    Div = diag(vecy);
    
    % 差商生成系数行向量vecx
    for k = 1:vecsize-1
        for i = 1:vecsize-k
            Div(i, i+k) = (Div(i+1, i+k) - Div(i, i+k-1))/(vect(i+k) - vect(i));
        end
    end
    vecx = Div(1, :);
    
    % 生成多项式函数polyfun
    syms f t; f = vecx(vecsize);
    for i = vecsize-1:-1:1
        f = vecx(i) + (t - vect(i)) * f;
    end
    polyfun = matlabFunction(f);
    end
    

      但不论如何,产生的结果完全一致。用同样的例子:

    vect=[-2, 0, 1]; vecy=[-27, -1, 0];
    
    % 命令行输入1,调用新增节点方法
    [polyfun, vecx] = newNewtonPI(vect, vecy)
    
    % 命令行输入2,调用差商方法
    [polyfun, vecx] = recNewtonPI(vect, vecy)
    
    % 命令行输出1/2,完全相同
    polyfun =
      包含以下值的 function_handle:
        @(t)-(t.*4.0-1.3e1).*(t+2.0)-2.7e1
    vecx =
       -27    13    -4
    

      容易检验,该多项式函数正是原数据点的多项式插值函数。

      

     

    转载于:https://www.cnblogs.com/gentle-min-601/p/9744395.html

    展开全文
  • Aitken(埃特金)逐次插值法 判断离散数据(xi,yi)(i=0,1,2,⋯ ,n)(x_i,y...这种在插值计算精度不够时增加节点(插值多项式的次数一般不宜超过6~8)以提高插值精度方法就是所谓的逐次插值法。 在上述情况中,运用Lagr

    Aitken(埃特金)逐次插值法

    判断离散数据 ( x i , y i ) ( i = 0 , 1 , 2 , ⋯   , n ) (x_i,y_i)(i=0,1,2,\cdots,n) (xi,yi)(i=0,1,2,,n)的插值精度,既可以采用事后误差估计的方法,也可以在插值点x的附近选取部分数据进行插值,然后再增加一些插值节点进行插值。若两次的插值结果之差小于规定的误差,则可认为插值精度复合要求而停止。这种在插值计算精度不够时增加节点(插值多项式的次数一般不宜超过6~8次)以提高插值精度方法就是所谓的逐次插值法。

    在上述情况中,运用Lagrange插值方法存在一个明显缺点,就是当插值节点发生变化和增加时,Lagrange插值公式中的所有基函数都得重新计算,即计算量大。由于插值节点发生变化和增加只是个别节点,因此只对发生变化和增加的结点进行计算以减小计算量十分重要。

    Aitken逐次插值法就是一种可以灵活地增加插值节点数,在前面计算结果的基础上继续进行计算而不必重新开始计算的方法。可见,Aitken逐次插值法具有承袭性的特点。

    先约定表示插值结果的符号,设在插值区间 [ a , b ] [a,b] [a,b]上,有n+1个顺序排列的插值节点 x 0 , ⋯   , x k , ⋯   , x n x_0,\cdots,x_k,\cdots,x_n x0,,xk,,xn,插值点为x。由前k个顺序排列的插值节点 x 0 , x 1 , ⋯   , x k − 1 x_0,x_1,\cdots,x_{k-1} x0,x1,,xk1构成的插值函数是x的k-1次多项式,可以用 f ( x 0 , ⋯   , x k − 1 ; x k − 1 ) f(x_0,\cdots,x_{k-1};x_{k-1}) f(x0,,xk1;xk1)表示,简记为 f k − 1 ( x k − 1 ) f_{k-1}(x_{k-1}) fk1(xk1)。在上述k个插值节点 x 0 , x 1 , ⋯   , x k − 1 x_0,x_1,\cdots,x_{k-1} x0,x1,,xk1的后面,再顺序增加一个新插值节点 x i ( i ≥ k ) x_i(i\geq k) xi(ik),进行k次插值。其插值函数是x的k次多项式,用 f ( x 0 , ⋯   , x k − 1 ; x k ) f(x_0,\cdots,x_{k-1};x_k) f(x0,,xk1;xk)表示,简记为 f k ( x i ) f_k(x_i) fk(xi),其中k表示插值次数, x k x_k xk为新增加的插值节点。在简记符号 f k ( x i ) f_k(x_i) fk(xi)中,k个顺序排列插值节点 x 0 , x 1 , ⋯   , x k − 1 x_0,x_1,\cdots,x_{k-1} x0,x1,,xk1中的最后一个节点 x k − 1 x_{k-1} xk1,由 f k ( x i ) f_k(x_i) fk(xi)下标k隐含地给出。

    1. 一次插值(线性插值)

    首先给出一个固定插值节点 x 0 x_0 x0及其函数值 f ( x 0 ) f(x_0) f(x0),再新增加一个节点 x i ( i ≥ 1 ) x_i(i\geq 1) xi(i1)(自然同时也给出其函数值 f ( x i ) f(x_i) f(xi)),用这两个插值节点进行线性插值,其结果为:
    f 1 ( x i ) = x − x i x 0 − x i f ( x 0 ) + x − x 0 x i − x 0 f ( x i ) ( i ≥ 1 ) f_1(x_i)=\frac{x-x_i}{x_0-x_i}f(x_0)+\frac{x-x_0}{x_i-x_0}f(x_i) \quad (i\geq 1) f1(xi)=x0xixxif(x0)+xix0xx0f(xi)(i1)
    若取 i = 1 i=1 i=1,则表示以 x 0 , x 1 x_0,x_1 x0,x1为节点进行一次插值,结果为:
    f 1 ( x 1 ) = x − x 1 x 0 − x 1 f ( x 0 ) + x − x 0 x 1 − x 0 f ( x 1 ) f_1(x_1)=\frac{x-x_1}{x_0-x_1}f(x_0)+\frac{x-x_0}{x_1-x_0}f(x_1) f1(x1)=x0x1xx1f(x0)+x1x0xx0f(x1)
    若取 i = 2 i=2 i=2,则表示以 x 0 x_0 x0 x 2 x_2 x2为节点进行一次插值,结果为:
    f 1 ( x 2 ) = x − x 2 x 0 − x 2 f ( x 0 ) + x − x 0 x 2 − x 0 f ( x 2 ) f_1(x_2)=\frac{x-x_2}{x_0-x_2}f(x_0)+\frac{x-x_0}{x_2-x_0}f(x_2) f1(x2)=x0x2xx2f(x0)+x2x0xx0f(x2)
    由此可得出如下结论: f 1 ( x i ) f_1(x_i) f1(xi)表示取固定节点 x 0 x_0 x0和变化节点 x i ( i ≥ 1 ) x_i(i\geq 1) xi(i1)及其相应的 f ( x 0 ) , f ( x 1 ) f(x_0),f(x_1) f(x0),f(x1)进行线性插值,得到关于x的一次多项式。

    1. 二次插值

    顺序给出两个插值节点 x 0 , x 1 x_0,x_1 x0,x1,再新增加一个节点 x i ( i ≥ 2 ) x_i(i\geq 2) xi(i2),用这3个插值节点进行插值,其结果为:
    f k ( x i ) = ( x − x 1 ) ( x − x 2 ) ( x 0 − x 1 ) ( x 0 − x 2 ) f ( x 0 ) + ( x − x 0 ) ( x − x 2 ) ( x 1 − x 0 ) ( x 1 − x 2 ) f ( x 1 ) + ( x − x 0 ) ( x − x 1 ) ( x 2 − x 0 ) ( x 2 − x 1 ) f ( x 2 ) f_k(x_i)=\frac{(x-x_1)(x-x_2)}{(x_0-x_1)(x_0-x_2)}f(x_0)+\frac{(x-x_0)(x-x_2)}{(x_1-x_0)(x_1-x_2)}f(x_1)+\frac{(x-x_0)(x-x_1)}{(x_2-x_0)(x_2-x_1)}f(x_2) fk(xi)=(x0x1)(x0x2)(xx1)(xx2)f(x0)+(x1x0)(x1x2)(xx0)(xx2)f(x1)+(x2x0)(x2x1)(xx0)(xx1)f(x2)
    若取i=2,则表示以 x 0 , x 1 x_0,x_1 x0,x1 x 2 x_2 x2为节点进行二次插值,其结果为:
    f 2 ( x 2 ) = ( x − x 1 ) ( x − x 2 ) ( x 0 − x 1 ) ( x 0 − x 2 ) f ( x 0 ) + ( x − x 0 ) ( x − x 2 ) ( x 1 − x 0 ) ( x 1 − x 2 ) f ( x 1 ) + ( x − x 0 ) ( x − x 1 ) ( x 2 − x 0 ) ( x 2 − x 1 ) f ( x 2 ) f_2(x_2)=\frac{(x-x_1)(x-x_2)}{(x_0-x_1)(x_0-x_2)}f(x_0)+\frac{(x-x_0)(x-x_2)}{(x_1-x_0)(x_1-x_2)}f(x_1)+\frac{(x-x_0)(x-x_1)}{(x_2-x_0)(x_2-x_1)}f(x_2) f2(x2)=(x0x1)(x0x2)(xx1)(xx2)f(x0)+(x1x0)(x1x2)(xx0)(xx2)f(x1)+(x2x0)(x2x1)(xx0)(xx1)f(x2)
    若取i=3,则表示以 x 0 , x 1 x_0,x_1 x0x1 x 3 x_3 x3为节点进行二次插值,其结果为:
    f 2 ( x 3 ) = ( x − x 1 ) ( x − x 3 ) ( x 0 − x 1 ) ( x 0 − x 3 ) f ( x 0 ) + ( x − x 0 ) ( x − x 3 ) ( x 1 − x 0 ) ( x 1 − x 3 ) f ( x 1 ) + ( x − x 0 ) ( x − x 3 ) ( x 3 − x 0 ) ( x 2 − x 1 ) f ( x 3 ) f_2(x_3)=\frac{(x-x_1)(x-x_3)}{(x_0-x_1)(x_0-x_3)}f(x_0)+\frac{(x-x_0)(x-x_3)}{(x_1-x_0)(x_1-x_3)}f(x_1)+\frac{(x-x_0)(x-x_3)}{(x_3-x_0)(x_2-x_1)}f(x_3) f2(x3)=(x0x1)(x0x3)(xx1)(xx3)f(x0)+(x1x0)(x1x3)(xx0)(xx3)f(x1)+(x3x0)(x2x1)(xx0)(xx3)f(x3)
    整理 f 2 ( x 2 ) f_2(x_2) f2(x2)得:
    f 2 ( x 2 ) = x − x 2 x 1 − x 2 f 1 ( x 1 ) + x − x 1 x 2 − x 1 f 1 ( x 2 ) f_2(x_2)=\frac{x-x_2}{x_1-x_2}f_1(x_1)+\frac{x-x_1}{x_2-x_1}f_1(x_2) f2(x2)=x1x2xx2f1(x1)+x2x1xx1f1(x2)
    同理
    f 2 ( x i ) = x − x i x 1 − x i f 1 ( x 1 ) + x − x 1 x i − x 1 f 1 ( x i ) f_2(x_i)=\frac{x-x_i}{x_1-x_i}f_1(x_1)+\frac{x-x_1}{x_i-x_1}f_1(x_i) f2(xi)=x1xixxif1(x1)+xix1xx1f1(xi)
    由此可得出如下结论: f 2 ( x i ) f_2(x_i) f2(xi)表示取固定节点 x 1 x_1 x1和变化节点 x i ( i ≥ k ) x_i(i\geq k) xi(ik)及其相应的 f 1 ( x 1 ) , f 1 ( x i ) f_1(x_1),f_1(x_i) f1(x1)f1(xi)进行线性插值,得到关于x的2次多项式。

    1. k次插值

    根据以上分析,可以推出如下结论:用两个 k − 1 k-1 k1次插值的结果 f k − 1 ( x k − 1 ) f_{k-1}(x_{k-1}) fk1(xk1) f k − 1 ( x i ) f_{k-1}(x_i) fk1(xi),进行线性插值,即可得到k次插值的结果 f k ( x i ) f_k(x_i) fk(xi)。即:

    f k ( x i ) f_k(x_i) fk(xi)表示固定节点 x k − 1 x_{k-1} xk1和变化节点 x i ( i ≥ k ) x_i(i\geq k) xi(ik)及其相应的 f k − 1 ( x k − 1 ) , f k − 1 ( x i ) f_{k-1}(x_{k-1}),f_{k-1}(x_i) fk1(xk1),fk1(xi)进行线性插值,从而得到关于x的k次多项式:
    f k ( x i ) = x − x i x k − 1 − x i f k − 1 ( x k − 1 ) + x − x k − 1 x i − x k − 1 f k − 1 ( x i ) , ( i ≥ k ) (10) f_k(x_i)=\frac{x-x_i}{x_{k-1}-x_i}f_{k-1}(x_{k-1})+\frac{x-x_{k-1}}{x_i-x_{k-1}}f_{k-1}(x_i), \quad (i\geq k) \tag{10} fk(xi)=xk1xixxifk1(xk1)+xixk1xxk1fk1(xi),(ik)(10)
    根据(10)式,可以计算出当 k = 1 , 2 , ⋯   , n ( i = k , k + 1 , ⋯   , n ) k=1,2,\cdots,n(i=k,k+1,\cdots,n) k=1,2,,n(i=k,k+1,,n)时的 f k ( x i ) f_k(x_i) fk(xi),由于计算 f k ( x i ) f_k(x_i) fk(xi)有很强的规律性,故将其排列成下表所示,该表称为Aitken插值表。从表中可以看出, f k ( x i ) f_k(x_i) fk(xi)是逐列计算出来的。这种足部提高插值次数以获得更高精度插值结果的插值方法称为Aitken逐次插值方法。

    在这里插入图片描述
    Aitken逐次插值法的特点是:

    (1)将一个高次插值过程归结为线性插值的多次重复。

    (2)插值表中的每个数据均为插值结果,从这些数据的一致程度可判断插值结果的精度,如果未达到精度要求,则在增加一个节点进行插值,直至满意为止。

    展开全文
  • 编程时要用到分段函数曲线的绘制方法:..+.*(分段条件)。 需要注意的是:函数表达式中的乘除和乘方都要加“.”。因为一般的函数都是数在乘变量运算。 x=-2:0.001:2; a=-0.5; w=abs(x); y=(1.5.*w.^3-2.5.*w.^...
  • 昨天在北理,顺道和老铁研究了一下B样条函数插值方法,发现网络上的帖子讲得不太好,看不懂。 今天在B站查找翻看相关内容,发现有个up主对其的讲解深入浅出,栩栩动人,终解疑惑。因此将其转载至B博士的个人CSDN上,...
  • b样条插值函数

    千次阅读 2019-05-07 10:55:34
    b样条插值函数是 spline = spapi(knots,x,y) spapi(k,x,y) spapi({knork1,…,knorkm},{x1,…,xm},y) spapi(…,‘noderiv’) % B样条曲线生成程序 % 说明:给定8个控制顶点{(3 5),(2 4),(3 2),(6 1),(5 8),(10 6)...
  • 常用的插值函数

    万次阅读 2013-05-02 10:04:44
    /*学校的数值分析课程正在讲插值函数,就趁着五一总结一下我所... 据维基百科,科学和工程问题可以通过诸如采样、实验等方法获得若干离散的数据,根据这些数据,我们往往希望得到个连续的函数(也就是曲线)或者更
  • 基函数与函数空间

    万次阅读 2018-09-05 16:10:57
    在学习线性回归模型的时候就会遇到基函数,可能我们会遇到多项式基函数、高斯基函数、sigmoid基函数,当然在高等数学和信号系统中还经常会碰到傅里叶基。有时候,不禁要问,这些基函数为什么这么设计?这些基函数的...
  • Matlab三均匀B样条曲线插值函数

    热门讨论 2014-03-25 11:25:23
    对给定的点进行三B样条插值,得到插值曲线,这里给定的点可以是二维平面上的点或三维点,注意输入的点矩阵要每行为个点坐标,里面都有注释,可以自己简单修改封装成自己想要的带参函数,里面有测试的点数据,...
  • Captain Dialog 2009-09-24 经过一次论坛的报告准备的确是让自己的软件设计思路清晰了不少,这几天一直忙着完成一个插值算法,突然发现,原来自己的数据管理方式已经很久没有审查了,今天大概画画数据流程图,发现...
  • MATLAB 牛顿插值函数

    千次阅读 2017-05-23 13:02:25
    X 为初始值 列向量 Y为初值函数值 列向量 x为插值点 M为插值次数function [y,R,A,C,L]=newdscg(X,Y,x,M) n=length(X);m=length(x); for t=1:m ...第列存Y的转置 s=0.0;p=1.0;q1=1.0;c1=1.0; for j=2:
  • 双三次插值(BiCubic插值)

    千次阅读 2018-11-26 11:13:02
    双三次插值(BiCubic插值 ) 双三次插值又称立方卷积插值。三次卷积插值是种更加复杂的插值方式。该算法利用待采样点周围16个点的灰度值作三次插值,...这种算法需要选取插值基函数来拟合数据,其最常用的插值基...
  • 双三次插值 又称立方卷积插值。 三次卷积插值是种更加复杂的插值方式。...这种算法需要选取插值基函数来拟合数据,其最常用的插值基函数如图1所示,本次实验采用如图所示函数作为基函数。 假设源图像A大...
  • 双三次插值 - 插值图像任意位置亚...这种算法需要选取插值基函数来拟合数据,其最常用的插值基函数如图1所示,本次实验采用如图所示函数作为基函数。 加权算法(a可以不取-0.5):点开链接 根据比例关系x/X=m/M=1/K
  • 高斯径向基函数(RBF)神经网络

    万次阅读 多人点赞 2019-04-03 00:53:52
    高斯径向基函数(RBF)神经网络 牛顿插值法-知乎 ...说径向基网络之前,先聊下径向基函数径向基函数(英语:radial basis function,缩写为RBF)是个取值仅依赖于到原点距离的实值函数,即 ϕ(x...
  • Lagrange插值函数及其Matlab代码

    千次阅读 2019-11-29 19:31:39
    为什么要引进插值函数 在实际问题中,两个变量的关系y=f(x)经常要 靠实验和观测来获得,而在通常的情况下只能得到f(x)在有限个点上的值 =(), i=0,1,2,...,n 人们希望找到f(x)的个近似函数y=(x),使得 ...
  • 利用Bézier曲线的端点插值性质,得到了构造三次插值样条曲线曲面的种改进的基函数——B基函数。由B基函数构造了C1形三次插值样条曲线;构造了C1双三次插值样条曲面
  • 神经网络学习笔记() RBF径向基函数神经网络 2018年08月06日 13:34:26吃机智豆长大的少女乙阅读数:2735 RBF径向基函数神经网络 初学神经网络,以下为综合其他博主学习材料及本人理解所得。 、RBF神经网络的...
  • MATLAB 插值

    千次阅读 2014-09-29 20:14:18
    MATLAB 插值   常用的多项式插值法: 1.拉格朗日插值 拉格朗日插值是簇插值的基本公式,由n+1个n次插值基函数构成,即n次插值多项式。 一般:
  • 在学习线性回归模型的时候就会遇到基函数,可能我们会遇到多项式基函数、高斯基函数、sigmoid基函数,当然在高等数学和信号系统中还经常会碰到傅里叶基。有时候,不禁要问,这些基函数为什么这么设计?这些基函数的...
  • 对正定核函数Φ:X×X→R\Phi:X\times X\to\mathbb{R}Φ:X×X→R和采样点X:={xi}i=1N⊂X\mathcal{X}\coloneqq\left\{x_i\right\}_{i=1}^N\subset XX:={xi​}i=1N​⊂X,其对应的插值函数应为 ui∗=∑j=1N[Φ−1]ijΦ...
  • 我们知道,有限元求解是有流程套路的 1. 伽辽金求解+Green-Guass公式。纽曼边界条件会自然满足,不需要后续边界条件修正 ...3.选择插值函数:一阶插值的公式为 4.将插值函数代入计算即可。 5....
  • 求解n次函数多项式之拉格朗日插值

    千次阅读 2018-08-21 10:44:15
    设F(x)=a0+a1x+a2x2+...+anxnF(x)=a0+a1x+a2x2+...+anxnF(x)=a_0+a_1x+a_2x^2+...+a_nx^n 在一些题目中例如找规律,比如是我上次在区域赛的时候如果可以真名或者猜...这里就用到了个东西叫拉格朗日插值法。 假如...
  • 径向基函数工作原理(样条函数)

    千次阅读 2013-01-17 12:07:00
    =================================...有以下五种基函数: 薄板样条函数 张力样条函数 规则样条函数 高曲面函数 反高曲面函数 在不同的插值表面中,每种基函数都有不同的形状和结果。RBF 方法是样条函数的个特...

空空如也

空空如也

1 2 3 4 5 ... 20
收藏数 3,938
精华内容 1,575
关键字:

一次插值基函数