在工程测量和科学实验中，所得到的数据通常都是离散的。如果要得到这些离散点以外的其他点的数值，就需要根据这些已知数据进行插值。例如，测量得 n n n 个点的数据为 ( x 1 , y 1 ) ， ( x 2 , y 2 ) ， … ， ( x n , y n ) (x_{1},y_{1})，(x_{2},y_{2})，\dots ，(x_{n},y_{n}) (x1,y1)，(x2,y2)，…，(xn,yn)，这些数据点反映了一个函数关系 y = f ( x ) y=f(x) y=f(x)，然而并不知道 f ( x ) f(x) f(x) 的解析式。
数据插值的任务就是根据上述条件构造一个函数 y = g ( x ) y=g(x) y=g(x)， 使得在 x ( i = 1 , 2 , . . . n ) x(i=1, 2, ... n) x(i=1,2,...n)，有 g ( x ) = f ( x ) g(x)=f(x) g(x)=f(x)，且在两个相邻的采样点 ( x i , x i + 1 ) ( i = 1 , 2 , … , n − 1 ) (x_{i},x_{i+1})(i=1,2,\dots ,n-1) (xi,xi+1)(i=1,2,…,n−1) 之间， g ( x ) g(x) g(x) 光滑过渡。如果被插值函数 f ( x ) f(x) f(x) 是光滑的，并且采样点足够密，一般在采样区间内， f ( x ) f(x) f(x) 与 g ( x ) g(x) g(x) 比较接近。插值函数 g ( x ) g(x) g(x) 一般由线性函数、多项式、样条函数或这些函数的分段函数充当。
根据被插值函数的自变量个数，插值问题分为一维插值、 二维插值和多维插值等；根据选用的插值方法，插值问题又分为线性插值、多项式插值和样条插值等。MATLAB 提供了一维、二维、 N N N 维数据插值函数 interp1、interp2 和 interpn，以及 3 次埃尔米特（Hermite）插值函数 pchip 和 3 次样条插值函数 spline 等。
函数根据 X 、 Y X、Y X、Y 的值，计算函数在 X 1 X1 X1 处的值。其中， X 、 Y X、Y X、Y 是两个等长的已知向量，分别描述采样点和采样值。若同一个采样点有多种采样值，则 Y Y Y 可以为矩阵， Y Y Y 的每一列对应一组采样。 X 1 X1 X1 是一个向量或标量，描述欲插值的点， Y 1 Y1 Y1 是一个与 X 1 X1 X1 等长的插值结果。method 用于指定插值方法，允许的取值有以下几种。
注意：对于超出 X 范围的插值点，使用 “linear” 和 “nearest” 插值方法，会给出 “NaN” 错误。对于其他的方法，将对超出范围的插值点执行外插值算法。
例如，以下概率积分的数据如下表所示，我们用不同的插值方法计算 f ( 0.472 ) f(0.472) f(0.472)。 f ( x ) = 2 π ∫ 0 x e − x 2 d x f(x)=\frac{2}{\sqrt{\pi}}\int_{0}^{x}e^{-x^{2}}\mathrm{d}x f(x)=π2∫0xe−x2dx
x x x
f ( x ) f(x) f(x)
0.46
0.4846555
0.47
0.4937542
0.48
0.5027498
0.49
0.5116683
这是一个一维插值问题，程序如下：
>> x=0.46:0.01:0.49;%给出x
>> f=[0.4846555,0.4937542,0.5027498,0.5116683];%给出f(x)>> format long%显示小数点后16位
>>interp1(x,f,0.472)%用默认方法，即线性插值计算f(0.472)ans =0.495553320000000>>interp1(x,f,0.472,'nearest')%用最近点插值计算f(0.472)ans =0.493754200000000>>interp1(x,f,0.472,'pchip')%用3次埃尔米特计算f(0.472)ans =0.495561119712056>>interp1(x,f,0.472,'spline')%用3次样条插值计算f(0.472)ans =0.495560736000000>> format short%小数点后显示4位
其中， X 、 Y X、Y X、Y 是两个向量，分别描述两个参数的采样点， Z Z Z 是与参数采样点对应的函数值， X 1 、 Y 1 X1、Y1 X1、Y1 是两个向量或标量，描述欲插值的点。 Z 1 Z1 Z1 是根据相应的插值方法得到的插值结果。
method 的取值与一维插值函数相同，但二维插值不支持 'pchip' 方法。 X 、 Y 、 Z X、Y、Z X、Y、Z 也可以是矩阵形式。同样，对于超出 X 、 Y X、Y X、Y 范围的插值点，使用 'linear' 和 'nearest' 插值方法，会给出 NaN 错误。对于 'spline' 方法，将对超出范围的插值点执行外插值算法。
例如，设 z = x 2 + y 2 z=x^{2}+y^{2} z=x2+y2，我们对 z z z 函数在 [ 0 ， 1 ] × [ 0 ， 2 ] [0，1]×[0，2] [0，1]×[0，2] 区域内进行插值。
ci =95.000030.20005.600000090.333347.400027.466716.00007.20004.000081.000058.866744.933333.200022.733317.666767.000064.600058.000051.600046.600041.0000
下图是根据插值结果 [ d i , t i , c i ] [d_{i},t_{i},c_{i}] [di,ti,ci]，用绘图函数 surf(di,ti,ci) 绘制的钢轨温度分布图。如果加密插值点，则绘制的分布图更理想。
为此，构造函数 y = g ( x ) y=g(x) y=g(x) 去逼近 f ( x ) f(x) f(x)，这里不要求曲线 g ( x ) g(x) g(x) 严格通过采样点，但希望 g ( x ) g(x) g(x) 能尽量地靠近这些点，就是使误差 δ i = g ( x ) − f ( x ) ( i = 1 , 2 ， … ， n ) δ_{i}=g(x)-f(x) (i=1, 2，…，n) δi=g(x)−f(x)(i=1,2，…，n) 在某种意义下达到最小。
MATLAB 曲线拟合的最优标准是采用常见的最小二乘原理，所构造的 g ( x ) g(x) g(x) 是一个次数小于插值节点个数的多项式。
我们假设测得 n n n 个离散数据点 ( x i , y i ) ( i = 1 , 2 ， . . . n ) (x_{i}, y_{i})(i=1, 2，... n) (xi,yi)(i=1,2，...n)，欲构造一个 m ( m ≤ n ) m(m≤n) m(m≤n) 次多项式 p ( x ) p(x) p(x)： p ( x ) = a m x m + a m − 1 x m − 1 + … + a 1 x + a 0 p(x)=a_{m}x^{m}+a_{m-1}x^{m-1}+…+a_{1}x+a_{0} p(x)=amxm+am−1xm−1+…+a1x+a0
所谓曲线拟合的最小二乘原理，就是使上述拟合多项式在各节点处的偏差 p ( x i ) − y i p(x_{i})-y_{i} p(xi)−yi 的平方和 ∑ i = 1 n ( p ( x i ) − y i ) 2 \sum_{i=1}^{n}(p(x_{i})-y_{i})^{2} ∑i=1n(p(xi)−yi)2 达到最小。数学上已经证明，上述最小二乘逼近问题的解总是确定的。
函数根据采样点 X X X 和采样点函数值 Y Y Y，产生一个 m m m 次多项式 P P P 及其在采样点的误差向量 S S S。
其中 X 、 Y X、Y X、Y 是两个等长的向量， P P P 是一个长度为 m + 1 m+1 m+1 的向量， P P P 的元素为多项式系数。mu 是一个二元向量，mu(1) 是函数 mean(X) 求向量 X X X 的平均值，而 mu(2) 是函数 std(X) 求向量 X X X 的方差。
可以用 polyval 函数按所得的多项式计算 X X X 各点上多项式的值。
例如，我们用一个 3 次多项式在区间 [ 0 , 2 π ] [0,2\pi] [0,2π] 内逼近函数 sin x \sin x sinx。
以上求得了 3 次拟合多项式 p ( x ) p(x) p(x) 的系数，得 p ( x ) = 0.0912 x 3 − 0.8596 x 2 + 1.8527 x − 0.1649 p(x)=0.0912x^{3}-0.8596x^{2}+1.8527x-0.1649 p(x)=0.0912x3−0.8596x2+1.8527x−0.1649。
下面，我们利用绘图函数得方法将多项式 p ( x ) p(x) p(x) 和 sin x \sin x sinx进行比较，程序如下：
绘制出 sin x \sin x sinx 和多项式 p ( x ) p(x) p(x) 在给定区间得函数曲线，如下图所示。其中虚线是 sin x \sin x sinx，实线是 p ( x ) p(x) p(x)。函数 polyval(P,X) 返回多项式 p ( x ) p(x) p(x) 得值。
三、数值微分
一般来说，函数的导数依然是一个函数。设函数 f ( x ) f(x) f(x) 的导函数 f ′ ( x ) = g ( x ) f'(x)=g(x) f′(x)=g(x)，高等数学关心的是 g ( x ) g(x) g(x) 的形式和性质，而数值分析关心的问题是怎样计算 g ( x ) g(x) g(x) 在一串离散点 X = ( x 1 , x 2 ， … ， x n ) X=(x_{1}, x_{2}，…，x_{n}) X=(x1,x2，…，xn) 的近似值 G = ( g 1 ， g 2 ， … ， g n ) G=(g_{1}，g_{2}，…，g_{n}) G=(g1，g2，…，gn) 以及所计算的近似值有多大误差。
1. 数值差分与差商
任意函数 f ( x ) f(x) f(x) 在 x x x 点的导数是通过极限定义的： f ′ ( x ) = lim h → 0 f ( x + h ) − f ( x ) h f'(x)=\lim_{h \to 0}\frac{f(x+h)-f(x)}{h} f′(x)=h→0limhf(x+h)−f(x) f ′ ( x ) = lim h → 0 f ( x ) − f ( x − h ) h f'(x)=\lim_{h \to 0}\frac{f(x)-f(x-h)}{h} f′(x)=h→0limhf(x)−f(x−h) f ′ ( x ) = lim h → 0 f ( x + h / 2 ) − f ( x − h / 2 ) h f'(x)=\lim_{h \to 0}\frac{f(x+h/2)-f(x-h/2)}{h} f′(x)=h→0limhf(x+h/2)−f(x−h/2)
上述式子中，均假设 h > 0 h>0 h>0，如果去掉上述等式右端的 h ⟶ 0 h\longrightarrow 0 h⟶0 的极限过程，并引进记号： △ f ( x ) = f ( x + h ) − f ( x ) \bigtriangleup f(x)=f(x+h)-f(x) △f(x)=f(x+h)−f(x) ▽ f ( x ) = f ( x ) − f ( x − h ) \bigtriangledown f(x)=f(x)-f(x-h) ▽f(x)=f(x)−f(x−h) δ f ( x ) = f ( x + h / 2 ) − f ( x − h / 2 ) \delta f(x)=f(x+h/2)-f(x-h/2) δf(x)=f(x+h/2)−f(x−h/2)
称 △ f ( x ) 、 ▽ f ( x ) \bigtriangleup f(x)、\bigtriangledown f(x) △f(x)、▽f(x) 及 δ f ( x ) \delta f(x) δf(x) 分别在函数 x x x 点处以 h ( h > 0 ) h(h>0) h(h>0) 为步长的向前差分、向后差分和中心差分。当步长为 h h h 充分小时，有 f ′ ( x ) ≈ △ f ( x ) h f'(x)\approx \frac{\bigtriangleup f(x)}{h} f′(x)≈h△f(x) f ′ ( x ) ≈ ▽ f ( x ) h f'(x)\approx \frac{\bigtriangledown f(x)}{h} f′(x)≈h▽f(x) f ′ ( x ) ≈ δ f ( x ) h f'(x)\approx \frac{\delta f(x)}{h} f′(x)≈hδf(x)
和差分一样，称 △ f ( x ) / h 、 ▽ f ( x ) / h \bigtriangleup f(x)/h、\bigtriangledown f(x)/h △f(x)/h、▽f(x)/h 及 δ f ( x ) / h \delta f(x)/h δf(x)/h 分别为函数在 x x x 点处以 h ( h > 0 ) h(h>0) h(h>0) 为步长的向前差商、向后差商和中心差商。
当步长 h ( h > 0 ) h(h>0) h(h>0) 充分小时，函数 f f f 在点 x x x 的微分接近于函数在该点的任意种差分，而 f f f 在点 x x x 的导数接近于函数在该点的任意种差商。
2. 数值微分的实现
有两种方式计算任意函数 f ( x ) f(x) f(x) 在给定点 x x x 的数值导数，第一种方式是用多项式或样条函数 g ( x ) g(x) g(x) 对 f ( x ) f(x) f(x) 进行逼近（插值或拟合），然后用逼近函数 g ( x ) g(x) g(x) 在点 x x x 处的导数作为 f ( x ) f(x) f(x) 在点 x x x 处的导数，第二种方式是用 f ( x ) f(x) f(x) 在点 x x x 处的某种差商作为其导数。
（1） DX=dif(X)：计算向量 X X X 的向前差分， D X ( i ) = X ( i + 1 ) − X ( i ) ， i = 1 ， 2 ， … ， n − 1 DX(i)=X(i+1)-X(i)， i=1，2，…，n-1 DX(i)=X(i+1)−X(i)，i=1，2，…，n−1。
（2） DX= diff(X,n)：计算 X X X 的 n n n 阶向前差分。例如，diff(X,2)=diff(diff(X))。
（3） DX-dif(A,n,dim)：计算矩阵 A A A 的 n n n 阶差分，dim=1 时（默认状态），按列计算差分；dim=2 时，按行计算差分。
例如，我们设 x x x 由 [ 0 , 2 π ] [0, 2\pi] [0,2π] 间均匀分布的 6 个点组成，求 sin x \sin x sinx 的 1~3 阶差分。
例如，设 f ( x ) = x 3 + 2 x 2 − x + 12 + x + 5 6 + 5 x + 2 f(x)=\sqrt{x^{3}+2x^{2}-x+12}+\sqrt[6]{x+5}+5x+2 f(x)=x3+2x2−x+12+6x+5+5x+2 我们用不同的方法求函数 f ( x ) f(x) f(x) 的数值导数，并在用一个坐标系中画出 f ′ ( x ) f'(x) f′(x) 的图像。
首先用一个 5 次多项式 p ( x ) p(x) p(x) 拟合函数 f ( x ) f(x) f(x)，并对 p ( x ) p(x) p(x) 求一般意义下的导数 d p ( x ) dp(x) dp(x)，求出 d p ( x ) dp(x) dp(x) 在假设点的值。
第二种方法直接求 f ( x ) f(x) f(x) 在假设点的数值导数。
第三种方法求出 f ′ ( x ) f'(x) f′(x)： f ′ ( x ) = 3 x 2 + 4 x − 1 2 x 3 + 2 x 2 − x + 12 + 1 6 x + 5 6 + 5 f'(x)=\frac{3x^{2}+4x-1}{2\sqrt{x^{3}+2x^{2}-x+12}}+\frac{1}{6\sqrt[6]{x+5}}+5 f′(x)=2x3+2x2−x+123x2+4x−1+66x+51+5
然后直接求 f ′ ( x ) f'(x) f′(x) 在假设点的导数，最后在同一坐标轴显示这 3 条曲线。
