//整理 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