Matlab数值计算差商与插值

均差定义 若已知函数$f(x)$在点$x_0,x_1,...x_n$处的值$f(x_0),f(x_1),...f(x_n).$如果$i≠j,$则
一阶均差$f[x_j,x_{j+1}]=\frac{f(x_{j+1})-f(x_j)}{x_{j+1}-x_j}(j=0,1,...n-1)$
二阶均差$f[x_j,x_{j+1},x_{j+2}]=\frac{f[x_{j+1},x_{j+2}]-f[x_j,x_{j+1}]}{x_{j+2}-x_j}(j=0,1,...n-2)$
n阶均差$f[x_0,x_1,...,x_n]=\frac{f[x_1,...,x_n]-f[x_0,...,x_{n-1}]}{x_n-x_0}$
例 由函数表求各阶均差

x	-2	-1	0	1	3
y	-56	-16	-2	-2	4

解：按公式计算一阶差商、二阶差商、三阶差商如下

x	f(x)	一阶差商	二阶差商	三阶差商
-2	-56
-1	-16	40
0	-2	14	-13
1	-2	0	-7	2
3	4	3	1	2

Matlab代码

clear
clc
x=[-2 -1 0 1 3]
y=[-56 -16 -2 -2 4]
deltx=diff(x);
delty=diff(y);
firstorder=delty./deltx %一阶
for i=1:length(x)-2
    delt2x(i)=x(i+2)-x(i);
end
delt2y=diff(firstorder);
secondorder=delt2y./delt2x %二阶
for i=1:length(x)-3
    delt3x(i)=x(i+3)-x(i);
end
delt3y=diff(secondorder);
thirdorder=delt3y./delt3x %三阶
for i=1:length(x)-4
    delt4x(i)=x(i+4)-x(i);
end
delt4y=diff(thirdorder);
fourorder=delt4y./delt4x %四阶

结果

x =

    -2    -1     0     1     3


y =

   -56   -16    -2    -2     4


firstorder =

    40    14     0     3


secondorder =

   -13    -7     1


thirdorder =

     2     2


fourorder =

     0

这里用到了diff，就再次介绍一下差分函数
补充：差分函数diff
diff(X) X为向量时(行列均可)，计算相邻两数的差[X(2)-X(1) X(3)-X(2) … X(n)-X(n-1)]
diff(X) X为矩阵时,计算矩阵的2~n行与1~n-1行的差，[X(2:n,:) - X(1:n-1,:)]
diff(X,N) 对上面函数diff(X)的扩充，这里的N指定N阶差分，二阶差分是对一阶差分的结果再做差分运算
DIFF(X,N,DIM) 对上面函数diff(X,N)的扩充，DIM取1或2，取1时按行差分，与上面结果一样，取2时按列差分

把上面的命令用字符串改造了一下，不过太难看懂了,no zuo no die
eval()函数的功能就是将括号内的字符串视为语句并运行，简单记为字符串转语句
num2str()函数的功能就是将括号内的数字转换为字符串，简单记为数字转字符串

clear
clc
x=[-2 -1 0 1 3]
y=[-56 -16 -2 -2 4]
deltx=diff(x);
delty=diff(y);
order1=delty./deltx %一阶
for j=2:4
  str=['for i=1:length(x)-',num2str(j),char(10)];
  str=[str,'delt',num2str(j),'x(i)=x(i+',num2str(j),')-x(i);',char(10)];
  str=[str,'end',char(10)];
  str=[str,'delt',num2str(j),'y=diff(order',num2str(j-1),');',char(10)];
  str=[str,'order',num2str(j),'=delt',num2str(j),'y./delt',num2str(j),'x',char(10)];
  eval(str)
end

结果

x =

    -2    -1     0     1     3


y =

   -56   -16    -2    -2     4


order1 =

    40    14     0     3


order2 =

   -13    -7     1


order3 =

     2     2


order4 =

     0

牛顿插值
牛顿插值公式及其余项
当$n=1$时：
差商$N_1(x)=f(x_0)+f[x_0,x_1](x-x_0)$
余项$R_1(x)=f[x,x_0,x_1](x-x_0)(x-x_1)$
当$n=2$时：
差商$N_2(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)$
余项$R_2(x)=f[x,x_0,x_1,x_2](x-x_0)(x-x_1)(x-x_2)$
$n$阶：
差商$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)$
$+f[x_0,x_1,...,x_n](x-x_0)(x-x_1)...(x-x_{n-1})$
余项$R_n(x)=f[x,x_0,x_1,...,x_n](x-x_0)(x-x_1)...(x-x_n)$
例 已知$x=1,4,9$的平方根为$1,2,3,$利用牛顿基本差商公式$\sqrt 7$的值

解:

$x_i$	$\sqrt {x_i}$	$f[x_i,x_{i+1}]$	$f[x_i,x_{i+1},x_{i+2}]$
1	1
4	2	$\frac{2-1}{4-1}=0.33333$
9	3	$\frac{3-2}{9-4}=0.2$	$\frac{0.2-0.33333}{9-1}=-0.016$

从而得二阶牛顿基本差商公式为
$P_2(x)=1+0.33333(x-1)-0.01667(x-1)(x-4)$
因此计算得$\sqrt 7$的近似值为$P_2(7)=2.69992$

clear
clc
x=[1 4 9]
y=[1 2 3]
deltx=diff(x);
delty=diff(y);
order1=delty./deltx %一阶
for i=1:length(x)-2
    delt2x(i)=x(i+2)-x(i);
end
delt2y=diff(order1);
order2=delt2y./delt2x %二阶
%%牛顿插值需要的值是y(1)、order1(1)、order2(1)、x(1)、x(2)
y(1),order1(1),order2(1),x(1),x(2)
%%构造多项式
P=[0 0 y(1)]+[0 order1(1)*poly(x(1))]+order2(1)*poly(x([1:2]))
%%求值
polyval(P,7)

结果

x =

     1     4     9


y =

     1     2     3


order1 =

    0.3333    0.2000


order2 =

   -0.0167


ans =

     1


ans =

    0.3333


ans =

   -0.0167


ans =

     1


ans =

     4


P =

   -0.0167    0.4167    0.6000


ans =

    2.7000

等距节点、差分与差商的关系
向前差分
一阶：$\varDelta f(k)=f(k+1)-f(k)$，如$\varDelta f(0)=f(1)-f(0)$，$\varDelta f(1)=f(2)-f(1)$
二阶：$\varDelta ^2 f(k)=\varDelta f(k+1)-\varDelta f(k)$，如$\varDelta ^2 f(0)=\varDelta f(1)-\varDelta f(0)$，$\varDelta ^2 f(1)=\varDelta f(2)-\varDelta f(1)$
三阶：$\varDelta ^3 f(k)=\varDelta^2 f(k+1)-\varDelta^2 f(k)$，如$\varDelta ^3 f(0)=\varDelta^2 f(1)-\varDelta^2 f(0)$，$\varDelta ^3 f(1)=\varDelta^2 f(2)-\varDelta^2 f(1)$
一般定义：$\varDelta ^m f(k)=\varDelta^{m-1} f(k+1)-\varDelta^{m-1} f(k)$，$m=2,3...$
此外还有向后差分、中心差分，这里暂时不做介绍
对于等距节点，差分与差商的关系
$f[x_k,x_{k+1},...x_{k+m}=]=\frac{\varDelta^m f(k)}{m!h^m}$
所以原来的牛顿插值公式在等距节点下，写成向前差分的形式就是
差商$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)$
$+f[x_0,x_1,...,x_n](x-x_0)(x-x_1)...(x-x_{n-1})$
$=f(x_0)+\frac{\varDelta f(0)}{h}(x-x_0)+\frac{\varDelta^2 f(0)}{2!h^2}(x-x_0)(x-x_0)(x-x_1)$
$+\frac{\varDelta^n f(0)}{n!h^n}(x-x_0)(x-x_1)...(x-x_{n-1})$