/*****************************************************************************
 *
 * deriv.cpp     Blitz++ Vector<T> example, illustrating expressions,
 *               ranges, and subvectors.
 *
 * $Id: deriv.cpp,v 1.1.1.1 2005/06/15 14:31:20 tveldhui Exp $
 *
 * $Log: deriv.cpp,v $
 * Revision 1.1.1.1  2005/06/15 14:31:20  tveldhui
 *
 *
 * Revision 1.1  1997/01/17 17:54:06  tveldhui
 * Initial revision
 */

#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.

    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