#include <blitz/array.h> #include <stdio.h> BZ_USING_NAMESPACE(blitz) BZ_USING_NAMESPACE(blitz::tensor) BZ_DECLARE_STENCIL3(deriv14,A,dA,h) dA = central14n(A) / h; BZ_END_STENCIL BZ_DECLARE_STENCIL3(deriv12,A,dA,h) dA = central12n(A) / h; BZ_END_STENCIL template<typename T_stencil> double calculateError(double h, T_stencil stencil) { const int N = 1024; Array<double,1> A(N), dA(N), x(N), hv(N), exactdA(N); x = 1/3.0 + (i - N/2.0) * h; hv = h; // Calculate the sampled function values A = cos(x*50); // Apply the finite difference approximation applyStencil(stencil, A, dA, hv); // Calculate the exact derivative exactdA = -50*sin(x*50); // Find the root mean-squared error (RMS), ignoring values near // the boundaries Range I(5,N-5); double rms = sqrt(sum(pow2(dA(I)-exactdA(I))) / I.length()); return rms; } template<typename T_stencil> double calculateAccuracyOrder(T_stencil stencil) { double h = 1/32.0; double oldrms = 0.0, oldh = 0.0; double order; for (int i=0; i < 9; ++i) { double rms = calculateError(h, stencil); if (i > 0) { order = (log(rms)-log(oldrms))/(log(h)-log(oldh)); printf("%18.10lf %18.10lf %18.10lf\n", h, rms, order); } oldh = h; h /= 2.0; oldrms = rms; } return order; } int main() { cout << "This program illustrates the O(h^2) and O(h^4) finite difference" << endl << "approximations to the first derivative." << endl << endl; cout << "central12: first derivative, O(h^2) error" << endl; cout.flush(); printf(" %-15s %-15s %-15s\n", "h", "Error", "Order"); calculateAccuracyOrder(deriv12()); cout << endl << endl << "central14: first derivative, O(h^4) error" << endl; cout.flush(); printf(" %-15s %-15s %-15s\n", "h", "Error", "Order"); calculateAccuracyOrder(deriv14()); return 0; }