牵着老婆满街逛

严以律己,宽以待人. 三思而后行.
GMail/GTalk: yanglinbo#google.com;
MSN/Email: tx7do#yahoo.com.cn;
QQ: 3 0 3 3 9 6 9 2 0 .

Blitz++ 计算二阶导数

// 整理 by RobinKin

#include  < blitz / vector.h >

using   namespace  blitz;

 int  main()
 {
     //  In this example, the function cos(x)^2 and its second derivative
     //  2 (sin(x)^2 - cos(x)^2) are sampled over the range [0,1).
     //  The second derivative is approximated numerically using a
     //  [ 1 -2  1 ] mask, and the approximation error is computed.

/*  cos(x)^2  的二阶导数 是2 (sin(x)^2 - cos(x)^2)

下面在[0,1) 的范围中，以delta为步长 ，用[ 1 -2  1 ] 的方法计算二阶导数

看看和精确值的误差
 */


     const   int  numSamples  =   100 ;               //  Number of samples
     double  delta  =   1 .  /  numSamples;           //  Spacing of samples
    Range R( 0 , numSamples  -   1 );               //  Index set of the vector

     //  Sample the function y = cos(x)^2 over [0,1)
     //
     //  An object of type Range can be treated as a vector, and used
     //  as a term in vector expressions.
     //
     //  The initialization for y (below) will be translated via expression
     //  templates into something of the flavour
     //
     //  for (unsigned i=0; i < 99; ++i)
     //  {
     //      double _t1 = cos(i * delta);
     //      y[i] = _t1 * _t1;
     //  }
    
    Vector < double >  y  =  sqr(cos(R  *  delta));

     //  Sample the exact second derivative
    Vector < double >  y2exact  =   2.0   *  (sqr(sin(R  *  delta))  -  sqr(cos(R  *  delta)));

     //  Approximate the 2nd derivative using a [ 1 -2  1 ] mask
     //  We can only apply this mask to the elements 1 .. 98, since
     //  we need one element on either side to apply the mask.
    Range I( 1 ,numSamples - 2 );
    Vector < double >  y2(numSamples);

    y2(I)  =  (y(I - 1 )  -   2   *  y(I)  +  y(I + 1 ))  /  (delta * delta);
  
     //  The above difference equation will be transformed into
     //  something along the lines of
     //
     //  double _t2 = delta*delta;
     //  for (int i=1; i < 99; ++i)
     //      y2[i] = (y[i-1] - 2 * y[i] + y[i+1]) / _t2;
 
     //  Now calculate the root mean square approximation error:

     double  error  =  sqrt(mean(sqr(y2(I)  -  y2exact(I))));
 
     //  Display a few elements from the vectors.
     //  This range constructor means elements 1 to 91 in increments
     //  of 15.
    Range displayRange( 1 ,  91 ,  15 );
 
    cout  <<   " Exact derivative: "   <<  y2exact(displayRange)  <<  endl
          <<   " Approximation:    "   <<  y2(Range(displayRange))  <<  endl
          <<   " RMS Error:        "   <<  error  <<  endl;

     return   0 ;
}


Output:

Exact derivative:[     - 1.9996    - 1.89847    - 1.62776    - 1.21164   - 0.687291   - 0.1015495
   ]
Approximation:   [    - 1.99953    - 1.89841     - 1.6277     - 1.2116   - 0.687269   - 0.1015468
   ]
RMS Error:        4.24826e-05

posted on 2006-07-01 13:23 杨粼波 阅读(693) 评论(0)  编辑 收藏 引用


只有注册用户登录后才能发表评论。
网站导航: 博客园   IT新闻   BlogJava   博问   Chat2DB   管理