计算数学

content

插值法

拉格朗日插值

通过平面上不同两点可以确定一条直线经过这两点，这就是拉格朗日线性插值问题，对于不在同一条直线的三个点得到的插值多项式为抛物线。这里给出一般的插值公式，拉格朗日插值的基多项式为：

l_{i} (x) = \prod_{\begin{matrix} j = 0 \\ j \neq i \end{matrix}}^{n} \frac{x - x_{i}}{x_{i} - x_{j}}, i = 0, 1, 2, \dots, n

对于没有接触过插值方法的读者，一个比较好的快速熟悉拉格朗日插值方法是把上面的多项式展开一下，看看如何计算基函数。有了基函数以后就可以直接构造插值多项式，插值多项式为

p_{n} = \sum_{i = 0}^{n} f (x_{i}) l_{i} (x)

虽然拉格朗日插值法较为基础，但因算法实现简单，且在一定条件下精度能得到保证，故而很多实际应用问题就采用该方法实现插值。

下面给出拉格朗日插值的C++代码。


void lagrange(const VEC &x,const VEC &y,const VEC &xx,VEC &yy)
/*------------------------------------------------------
   Versions and Changes :
      v1.0----2015-8-11
        lagrange 插值方法， 对多点一次性完成插值
        C++语言中，允许 xx为 VEC(1)   yy为 VEC(1)
--------------------------------------------------------
paramtegers :
 *  x-----------样本数据
 *  y-----------样本函数值
 *  xx----------要计算的点
 *  返回值为-----插值结果，维数同xx
Ref . :
     scientific computation in C#  Tsinghua Univ. Press
------------------------------------------------------
   Author      : Song Yezhi  
   Copyrigt(C) : Shanghai Astronomical Observatory, CAS
                (All rights reserved)             2015
-------------------------------------------------------*/
{
    //获得样本维数
    int N = x.size();
    //获得被插值点的数据
    int M = xx.size();    
    int p,i;
    for ( p = 0; p < M; p++)
    //对所有的要插值的点做循环
    {
        //至各分式之和为0
        double sum = 0;
        for (i = 0; i < N; i++)
        {
            double tmp = y(i);
            for (int k = 0; k < N; k++)
            {
                if (k == i) continue;
                //如果k与i相等，则退出当次循环
                //但并不退出循环体
                tmp = tmp * (xx(p) - x(k)) / (x(i) - x(k));
            }
            sum = sum + tmp;        }
        //结果赋值给要计算的点
        yy(p) = sum;
     }

牛顿插值

拉格朗日插值多项式每增加一个插值基点，原先计算的插值多项式

p_{n} (x)

对

p_{n + 1} (x)

没有用，这样必然增加计算工作量，尤其在没有计算机的情况下这个问题更加突出。我们期望增加插值基点时原先的计算结果对后面的计算仍然有用. 在介绍牛顿插值之前先要了解一下差商的基本概念。

\begin{array}{l} f [x_{k}] = f (x_{k}) \\ f [x_{k - 1}, x_{k}] = \frac{f [x_{k}] - f [x_{k - 1}]}{x_{k - 1} - x_{k}} \\ f [x_{k - 2}, x_{k - 1}, x_{k}] = \frac{f [x_{k - 1}, x_{k}] - f [x_{k - 2}, x_{k - 1}]}{x_{k - 2} - x_{k}} \\ ⋮ \\ f [x_{k - j}, x_{k - j + 1}, \dots, x_{k}] = \frac{f [x_{k - j + 1}, \dots, x_{k}] - f [x_{k - j}, \dots, x_{k - 1}]}{x_{k} - x_{k - j}} \end{array}

有了上面的准备，便可以给出牛顿插值多项式为：

\begin{array}{l} N_{n} (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}) + \dots + f [x_{0}, x_{1}, \dots, x_{n}] \\ (x - x_{0}) (x - x_{1}) \dots (x - x_{n - 1}) \end{array}

如果记

w_{k} = (x - x_{0}) (x - x_{1}) \dots (x - x_{k - 1}), k = 1, 2, \dots, n

则牛顿插值可以表达为：

\begin{array}{l} N_{n} (x) = f [x_{0}] + f [x_{0}, x_{1}] w_{1} + f [x_{0}, x_{1}, x_{2}] w_{2} \\ + \dots + f [x_{0}, x_{1}, \dots, x_{n}] w_{n} \end{array}

下面给出牛顿插值的C#源代码。


using System;
using System.Collections.Generic;
using System.Linq;
using System.Text;
using System.IO;
class cnewton
/*--------------------------------------class comment
 Version   :  V1.0
 Coded by  :  syz
 Date      :  2011-07-16 21:33:04          *星期六*
----------------------------------------------------
Desciption :  牛顿插值法
 *
parameters :
 *   fout--------输出文件
 * 
Methods    :
 *  newton-------牛顿插值法（对多数据一次性完成插值）
 *  newton0------对单个点插值法
 *  drive--------驱动程序
--------------------------------------------------*/
{
    public StreamWriter fout = new StreamWriter("fout.txt");
    public VEC newton(VEC x, VEC y, VEC xx)
    /*------------------------------------------ comment
    Author      : syz 
    Date        : 2011-07-16 21:34:30         
    ---------------------------------------------------
    Desciption  : 牛顿插值法
     *
    Post Script :
     *   详细算法参看《Fortran科学计算与工程》
    paramtegers :
     *  x--------样本自变量
     *  y--------样本函数值
     *  xx-------要插值的向量
     *  返回值-------插值结果（维数同xx）
    -------------------------------------------------*/
    {
        int N = x.dim;
        int M = xx.dim;
        VEC yy=new VEC(M);
        for (int i = 0; i < M; i++)
            yy[i] = newton0(x, y, xx[i]);
        return yy;

    }
    public double newton0(VEC x, VEC y, double t)
    /*------------------------------------------ comment
    Author      : syz 
    Date        : 2011-07-16 21:36:41         
    ---------------------------------------------------
    Desciption  :  对单点牛顿插值方法
     *
    Post Script :
     *     编程时需要注意向量与矩阵的指标
    paramtegers :
     * x,y-------样本自变量与应变量
     * t---------要插值的点
    -------------------------------------------------*/
    {
        //获取样本维数
        int N = x.dim;        
        //声明矩阵
        MAT Q = new MAT(N, N);
        VEC b = new VEC(N);
        //把样本函数值复制给Q第一列
        for (int i = 0; i < N; i++)
            Q[i, 0] = y[i];
        for (int i = 1; i < N; i++)
            for (int j = 1; j <= i; j++)
                Q[i, j] = (Q[i, j - 1] - Q[i - 1, j - 1]) / (x[i] - x[i - j]);
        b[N - 1] = Q[N - 1, N - 1];
        for (int k = N - 1; k >= 1; k--)
            b[k - 1] = Q[k - 1, k - 1] + b[k] * (t - x[k - 1]);
        //返回运算结果
        return b[0];

    }
}

Hermit插值

有些实际问题，不仅要求插值多项式在插值点与函数值相同，而且还要求导数相同，这类问题就是Hermit插值问题。

如果 $f (x)$ 在区间[a,b]上连续可以导， $x_{0}, x_{1}, \dots x_{n} \in [a, b]$ 是互异的，那么存在唯一的多项式 $H_{2 n + 1} (x)$ 满足多项式在这些点上的值与函数 $f (x)$ 的值相等、多项式在这些点的一阶导数值与函数的一阶导数值相等。

这个多项式可以表示为

\begin{array}{l} H_{2 n + 1} (x) = \sum_{i = 0}^{n} f (x_{i}) [1 - 2 (x - x_{i}) {l_{i}}^{'} (x_{i})] {l_{i}}^{2} (x) \\ + \sum_{i = 0}^{n} f^{'} (x_{i}) (x - x_{i}) {l_{i}}^{2} (x) \end{array}

其中

l_{i} (x) = \prod_{\begin{matrix} j = 0 \\ j \neq i \end{matrix}}^{n} \frac{(x - x_{i})}{(x_{i} - x_{j})}, i = 0, 1, \dots, n

{l_{i}}^{'} (x_{i}) = \sum_{\begin{matrix} j = 0 \\ j \neq i \end{matrix}}^{n} \frac{1}{x_{i} - x_{j}}, i = 0, 1, \dots, n

这就是Hermite插值。


using System;
using System.Collections.Generic;
using System.Linq;
using System.Text;
using System.IO;
class chermite
/*--------------------------------------class comment
 Version   :  V1.0
 Coded by  :  syz
 Date      :  2011-07-16 22:31:09          *星期六*
----------------------------------------------------
Desciption : Hermite插值
 *
parameters :
 *  fout-----------输出文件对象
 *
Methods    :
 *  Hermite--------插值方法
 *  Hermite--------对一点的插值方法
 *  drive----------驱动函数
 *
--------------------------------------------------*/
{
    public StreamWriter fout = new StreamWriter("fout.txt");
    public VEC hermite(VEC x, VEC y, VEC dy, VEC xx)
    /*------------------------------------------ comment
    Author      : syz
    Date        : 2011-07-16 22:32:24
    ---------------------------------------------------
    Desciption  : Hermite插值方法
     *
    Post Script :
     *
    paramtegers :
     *  x--------样本自变量
     *  y--------样本应变量
     *  dy-------样本导数值
     *  xx-------待插值点
     *  返回值--------插值结果（维数同xx）
    -------------------------------------------------*/
    {
        //获取要插值数据向量的维数
        int M = xx.dim;

        //声明yy为新的向量并分配地址
        VEC yy=new VEC(M);

        for (int i = 0; i < M; i++)
            yy[i] = hermite0(x, y, dy,xx[i]);

        return yy;

    }
    public double hermite0(VEC x, VEC y, VEC dy, double xx)
    /*------------------------------------------ comment
    Author      : syz
    Date        : 2011-07-16 22:33:41
    ---------------------------------------------------
    Desciption  : 对单点的插值方法
     *
    Post Script :
     *
    paramtegers :
     *   x-------样本自变量
     *   y-------样本应变量
     *   dy------样本导数值
     *   xx------要插值点
     *   返回值-------插值结果
    -------------------------------------------------*/
    {
        //获取x的维数
        int N = x.dim;

        //声明L和dL并分配地址
        VEC L=new VEC(N);
        VEC dL=new VEC(N);

        //该段循环求得L
        for (int i = 0; i < N ; i++)
        {
            L[i] = 1.0;
            for (int j = 0; j < N; j++)
            {
                //如果j和i相等，退出当前循环，
                //但并不退出循环体
                if (j == i)
                    continue;

                L[i] = L[i] * (xx - x[j]) / (x[i] - x[j]);

            }

        }
        //该段循环求得dL
        for (int i = 0; i < N; i++)
        {
            dL[i] = 0;
            for (int j = 0; j < N; j++)
            {
                if (j == i)
                    continue;
                dL[i] = dL[i]+1.0 / (x[i] - x[j]);
            }
        }
        double  yy=0;
        for (int i = 0; i < N; i++)
            yy = yy + y[i] * (1.0 - 2.0* (xx - x[i]) *
             dL[i]) * L[i] * L[i]
                + dy[i] * (xx - x[i]) * L[i] * L[i];
        //返回计算结果
        return yy;
    }
}