Matlab数值计算差商与插值
2015-06-25 11:24
369 查看
均差定义 若已知函数f(x)f(x)在点x0,x1,...xnx_0,x_1,...x_n处的值f(x0),f(x1),...f(xn).f(x_0),f(x_1),...f(x_n).如果i≠j,i≠j,则
一阶均差f[xj,xj+1]=f(xj+1)−f(xj)xj+1−xj(j=0,1,...n−1)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[xj,xj+1,xj+2]=f[xj+1,xj+2]−f[xj,xj+1]xj+2−xj(j=0,1,...n−2)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[x0,x1,...,xn]=f[x1,...,xn]−f[x0,...,xn−1]xn−x0f[x_0,x_1,...,x_n]=\frac{f[x_1,...,x_n]-f[x_0,...,x_{n-1}]}{x_n-x_0}
例 由函数表求各阶均差
解:按公式计算一阶差商、二阶差商、三阶差商如下
Matlab代码
结果
这里用到了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()函数的功能就是将括号内的数字转换为字符串,简单记为数字转字符串
结果
牛顿插值
牛顿插值公式及其余项
当n=1n=1时:
差商N1(x)=f(x0)+f[x0,x1](x−x0)N_1(x)=f(x_0)+f[x_0,x_1](x-x_0)
余项R1(x)=f[x,x0,x1](x−x0)(x−x1)R_1(x)=f[x,x_0,x_1](x-x_0)(x-x_1)
当n=2n=2时:
差商N2(x)=f(x0)+f[x0,x1](x−x0)+f[x0,x1,x2](x−x0)(x−x1)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)
余项R2(x)=f[x,x0,x1,x2](x−x0)(x−x1)(x−x2)R_2(x)=f[x,x_0,x_1,x_2](x-x_0)(x-x_1)(x-x_2)
nn阶:
差商Nn(x)=f(x0)+f[x0,x1](x−x0)+f[x0,x1,x2](x−x0)(x−x1)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[x0,x1,...,xn](x−x0)(x−x1)...(x−xn−1)+f[x_0,x_1,...,x_n](x-x_0)(x-x_1)...(x-x_{n-1})
余项Rn(x)=f[x,x0,x1,...,xn](x−x0)(x−x1)...(x−xn)R_n(x)=f[x,x_0,x_1,...,x_n](x-x_0)(x-x_1)...(x-x_n)
例 已知x=1,4,9x=1,4,9的平方根为1,2,3,1,2,3,利用牛顿基本差商公式7√\sqrt 7的值
解:
从而得二阶牛顿基本差商公式为
P2(x)=1+0.33333(x−1)−0.01667(x−1)(x−4)P_2(x)=1+0.33333(x-1)-0.01667(x-1)(x-4)
因此计算得7√\sqrt 7的近似值为P2(7)=2.69992P_2(7)=2.69992
结果
等距节点、差分与差商的关系
向前差分
一阶:Δf(k)=f(k+1)−f(k)\varDelta f(k)=f(k+1)-f(k),如Δf(0)=f(1)−f(0)\varDelta f(0)=f(1)-f(0),Δf(1)=f(2)−f(1)\varDelta f(1)=f(2)-f(1)
二阶:Δ2f(k)=Δf(k+1)−Δf(k)\varDelta ^2 f(k)=\varDelta f(k+1)-\varDelta f(k),如Δ2f(0)=Δf(1)−Δf(0)\varDelta ^2 f(0)=\varDelta f(1)-\varDelta f(0),Δ2f(1)=Δf(2)−Δf(1)\varDelta ^2 f(1)=\varDelta f(2)-\varDelta f(1)
三阶:Δ3f(k)=Δ2f(k+1)−Δ2f(k)\varDelta ^3 f(k)=\varDelta^2 f(k+1)-\varDelta^2 f(k),如Δ3f(0)=Δ2f(1)−Δ2f(0)\varDelta ^3 f(0)=\varDelta^2 f(1)-\varDelta^2 f(0),Δ3f(1)=Δ2f(2)−Δ2f(1)\varDelta ^3 f(1)=\varDelta^2 f(2)-\varDelta^2 f(1)
一般定义:Δmf(k)=Δm−1f(k+1)−Δm−1f(k)\varDelta ^m f(k)=\varDelta^{m-1} f(k+1)-\varDelta^{m-1} f(k),m=2,3...m=2,3...
此外还有向后差分、中心差分,这里暂时不做介绍
对于等距节点,差分与差商的关系
f[xk,xk+1,...xk+m=]=Δmf(k)m!hmf[x_k,x_{k+1},...x_{k+m}=]=\frac{\varDelta^m f(k)}{m!h^m}
所以原来的牛顿插值公式在等距节点下,写成向前差分的形式就是
差商Nn(x)=f(x0)+f[x0,x1](x−x0)+f[x0,x1,x2](x−x0)(x−x1)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[x0,x1,...,xn](x−x0)(x−x1)...(x−xn−1)+f[x_0,x_1,...,x_n](x-x_0)(x-x_1)...(x-x_{n-1})
=f(x0)+Δf(0)h(x−x0)+Δ2f(0)2!h2(x−x0)(x−x0)(x−x1)=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)
+Δnf(0)n!hn(x−x0)(x−x1)...(x−xn−1)+\frac{\varDelta^n f(0)}{n!h^n}(x-x_0)(x-x_1)...(x-x_{n-1})
一阶均差f[xj,xj+1]=f(xj+1)−f(xj)xj+1−xj(j=0,1,...n−1)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[xj,xj+1,xj+2]=f[xj+1,xj+2]−f[xj,xj+1]xj+2−xj(j=0,1,...n−2)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[x0,x1,...,xn]=f[x1,...,xn]−f[x0,...,xn−1]xn−x0f[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 |
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=1n=1时:
差商N1(x)=f(x0)+f[x0,x1](x−x0)N_1(x)=f(x_0)+f[x_0,x_1](x-x_0)
余项R1(x)=f[x,x0,x1](x−x0)(x−x1)R_1(x)=f[x,x_0,x_1](x-x_0)(x-x_1)
当n=2n=2时:
差商N2(x)=f(x0)+f[x0,x1](x−x0)+f[x0,x1,x2](x−x0)(x−x1)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)
余项R2(x)=f[x,x0,x1,x2](x−x0)(x−x1)(x−x2)R_2(x)=f[x,x_0,x_1,x_2](x-x_0)(x-x_1)(x-x_2)
nn阶:
差商Nn(x)=f(x0)+f[x0,x1](x−x0)+f[x0,x1,x2](x−x0)(x−x1)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[x0,x1,...,xn](x−x0)(x−x1)...(x−xn−1)+f[x_0,x_1,...,x_n](x-x_0)(x-x_1)...(x-x_{n-1})
余项Rn(x)=f[x,x0,x1,...,xn](x−x0)(x−x1)...(x−xn)R_n(x)=f[x,x_0,x_1,...,x_n](x-x_0)(x-x_1)...(x-x_n)
例 已知x=1,4,9x=1,4,9的平方根为1,2,3,1,2,3,利用牛顿基本差商公式7√\sqrt 7的值
解:
xix_i | xi−−√\sqrt {x_i} | f[xi,xi+1]f[x_i,x_{i+1}] | f[xi,xi+1,xi+2]f[x_i,x_{i+1},x_{i+2}] |
1 | 1 | ||
4 | 2 | 2−14−1=0.33333\frac{2-1}{4-1}=0.33333 | |
9 | 3 | 3−29−4=0.2\frac{3-2}{9-4}=0.2 | 0.2−0.333339−1=−0.016\frac{0.2-0.33333}{9-1}=-0.016 |
P2(x)=1+0.33333(x−1)−0.01667(x−1)(x−4)P_2(x)=1+0.33333(x-1)-0.01667(x-1)(x-4)
因此计算得7√\sqrt 7的近似值为P2(7)=2.69992P_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
等距节点、差分与差商的关系
向前差分
一阶:Δf(k)=f(k+1)−f(k)\varDelta f(k)=f(k+1)-f(k),如Δf(0)=f(1)−f(0)\varDelta f(0)=f(1)-f(0),Δf(1)=f(2)−f(1)\varDelta f(1)=f(2)-f(1)
二阶:Δ2f(k)=Δf(k+1)−Δf(k)\varDelta ^2 f(k)=\varDelta f(k+1)-\varDelta f(k),如Δ2f(0)=Δf(1)−Δf(0)\varDelta ^2 f(0)=\varDelta f(1)-\varDelta f(0),Δ2f(1)=Δf(2)−Δf(1)\varDelta ^2 f(1)=\varDelta f(2)-\varDelta f(1)
三阶:Δ3f(k)=Δ2f(k+1)−Δ2f(k)\varDelta ^3 f(k)=\varDelta^2 f(k+1)-\varDelta^2 f(k),如Δ3f(0)=Δ2f(1)−Δ2f(0)\varDelta ^3 f(0)=\varDelta^2 f(1)-\varDelta^2 f(0),Δ3f(1)=Δ2f(2)−Δ2f(1)\varDelta ^3 f(1)=\varDelta^2 f(2)-\varDelta^2 f(1)
一般定义:Δmf(k)=Δm−1f(k+1)−Δm−1f(k)\varDelta ^m f(k)=\varDelta^{m-1} f(k+1)-\varDelta^{m-1} f(k),m=2,3...m=2,3...
此外还有向后差分、中心差分,这里暂时不做介绍
对于等距节点,差分与差商的关系
f[xk,xk+1,...xk+m=]=Δmf(k)m!hmf[x_k,x_{k+1},...x_{k+m}=]=\frac{\varDelta^m f(k)}{m!h^m}
所以原来的牛顿插值公式在等距节点下,写成向前差分的形式就是
差商Nn(x)=f(x0)+f[x0,x1](x−x0)+f[x0,x1,x2](x−x0)(x−x1)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[x0,x1,...,xn](x−x0)(x−x1)...(x−xn−1)+f[x_0,x_1,...,x_n](x-x_0)(x-x_1)...(x-x_{n-1})
=f(x0)+Δf(0)h(x−x0)+Δ2f(0)2!h2(x−x0)(x−x0)(x−x1)=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)
+Δnf(0)n!hn(x−x0)(x−x1)...(x−xn−1)+\frac{\varDelta^n f(0)}{n!h^n}(x-x_0)(x-x_1)...(x-x_{n-1})
相关文章推荐
- matlab乘与点乘的区别 (*与.* ^与.^)
- matlab中用imshow()显示double类型图像中出现的问题
- matlab中norm函数的用法
- matlab echo 的用法
- 实验二:FFT算法的MATLAB实现
- 实验一 离散时间序列卷积和MATLAB实现
- Huffman编码用MTLAB的实现及编码注释----------Matlab
- MIMO信道容量及注水算法---------Matlab
- MATLAB实用源代码
- 协方差,方差,期望的意义
- matlab中sum()求和函数
- PCA ( 主成分分析) 详解 ( 写给初学者) 结合matlab
- Matlab常用函数流水账
- 几种简单常用的镜头边缘检测算法(matlab实现)
- matlab中 tic,toc函数的用法
- 利用MATLAB绘制信号时域波形和信号的频谱
- [图像]Canny检测的Matlab实现(含代码)
- matlab中clc,close,close all,clear,clear all作用区别
- Matlab——zeros函数和ones函数
- win7运行matlab7.0闪一下就消失