数值微积分
引言
  在微分中，函数的导数是用极限来定义的，如果一个函数是以数值给出的离散形式，那么它的导数就无法用极限运算方法求得，当然也就更无法用求道方法去计算函数在某点处的导数。
  一般来说，函数的导数依然是一个函数。设函数f(x)的导数f′(x)=g(x),高等数学关心的是g（x）的形式和性质，而数值分析关心的问题是怎样的计算g（x）在一串离散点X=（x1,x2,…xn）的近似值G=（g1,g2,..gn）以及所计算的近似值有多大的误差。
1.    数值差分与差商
任意函数f(x)在x点的导数是通过极限定义的：
            f′(x)=lim f(x+h)-f(x)/h
            f′(x)=lim f(x)-f(x-h)/h
            f′(x)=lim f(x+h/2)-f(x-h/2)/h
 
上述式子中，均假设h﹥0，如果去掉上述等式右端的h→0的极限过程，并引进记号：
△f(x)=f(x+h)-f(x)
▽f(x)=f(x)-f(x-h)
δf(x)=f(x+h/2)-f(x-h/2)
 
称△f(x) ▽f(x)、δf(x)分别为函数在x点的以h（h>0）为步长的向前差分、向后差分、向中差分。当步长h充分的小时，有
 
            f′(x）≈△f(x)/h
diff函数
            f′(x）≈▽f(x)/h
            f′(x）≈δf(x)/h
和差分一样，称△f(x)/h, ▽f(x)/h,、δf(x)/h分别为函数在x点一h（h>0）为步长的向前差商、向后差商、向中差商。当步长h充分的小时，函数f在x点微分接近于函数在该点的任意种差分，而f在点x的导数接近于函数在该点的任意种差商。
 
2.    数值微分的实现
有两种方式计算任意函数f（x）在给定点x点数值导数。第一种方式是用多项式或样条函数g（x）对f(x)进行逼近，然后用逼近函数g(x)在点x的导数作为f(x)在点x的导数。
第二种方式是用f（x）在点x处的某种差商作为其导数。在MATLAB中，没有直接提供数值导数的函数，只有计算向前差分的函数diff，其调用格式为：
  Diff(X),输入参数x可为矢量或矩阵。若X为矢量，则返回[X(2)-X(1),X(3)-X(2),…,X(n)-X(n-
1)];若X为矩阵，则返回矩阵每行的差分，即[X(2:m，：)-X(1:m-1,:)].通常，diff(X)返回沿X的第一个成对维的差值。
  Y=diff(X,n):递归调用diff函数n次，生成n次差分。
  Y=diff(X,n,dim):  返回X中第dim维的n次差分或导数值
 
设计1如下：
  X=[3 7 5;0 4 2];
  Diff(x)
  ans=
      -3  -3  -3
 
设计2如下：
  x=[3 7 5;0 4 2];
  Diff=(x,2)
  ans=
0                0
 
 
 
设计3如下：
x=[3 7 5；0 4 2]；
diff=(x,2,1)
ans=
empty matrix: 0-by-3
 
 
设：
     
用不同的方法求函数f(x)的数值导数，并在同一个坐标系中画出f′(x)的图形
 
程序如下：
f=inline(‘sqrt(x.^3+2*x.^2-x+12)+(x+5).^(1/6)+5*x+2’);
g=inline(‘(3*x.^2+4*x-1)./sqrt(x.^3+2*x.^2-x+12)/2+1/6./(x+5).^(5/6)+5’);
x=-3:0.01:3;
p=polyfit(x,f(x),5);                            ﹪用5次多项式p拟合f(x)
dp=polyder(p);                                  ﹪对拟合多项式p求导数dp
dpx=polyval(dp,x);                              ﹪求dp在假设点的函数值
dx=diff(f([x,3.01]))/0.01;                      ﹪直接对f(x)求数值导数
gx=g(x);                                        ﹪求函数f的导函数g在假设点的导数
plot(x,dpx,x,dx,’.’,x,gx,’-‘);              ﹪作图
 
  程序运行后得到图6.3所示的图形。结果表明，用3种方法求得的数值导数比较近似。
 
 
 
 
 
 
 
 
 
 
 
MATLAB语言提供了计算给定向量差分的函数diff(),其调用的方法是dy=diff(y).假设向量y是｛yi｝,i=1,2,3,…,n构成，则经过diff()函数处理后将得出一个新的向量：
｛yi+1-yi｝,i=1,2,3,…,n,显然新得出的向量长度比原来向量的长度y要小1.如：
 
>>v=vander(1:6)
  v=
      1    1    1    1    1    1 
      32    16    8    4    2    1
    243    81    27    9    3    1
    1024    256  64    16    4    1
    3125    625  125  25    5    1
    7776    1296  216  36    6    1
 
  >>diff(v)
  ans=
    31    15    7    3      1    0
    211    65    19    5      1    0
    781    175    37    7      1    0
    2101    369    61    9      1    0
    4651    671    91    11    1    0

版权声明：本站内容均来自互联网，仅供演示用，请勿用于商业和其他非法用途。如果侵犯了您的权益请与我们联系QQ:729038198，我们将在24小时内删除。