Demonstration of the interpolated derivative, linear vs. cubic. The linear derivative is not smooth, whereas the cubic derivative is second-order smooth.
Size: 13x10x23 value at {5.1,6.2,6.7} Ipol: Linear: 0.72 Cubic: 0.975172 Size: 13x10x23 value at {5.1,6.2,6.8} Ipol: Linear: 0.72 Cubic: 0.730266 Size: 13x10x23 value at {5.1,6.2,6.9} Ipol: Linear: 0.72 Cubic: 0.405208 Size: 13x10x23 value at {5.1,6.2,7} Ipol: Linear: -0.72 Cubic: 0 Size: 13x10x23 value at {5.1,6.2,7.1} Ipol: Linear: -0.72 Cubic: -0.405208 Size: 13x10x23 value at {5.1,6.2,7.2} Ipol: Linear: -0.72 Cubic: -0.730266
#include <vector/DeriveArray.hpp> #include <iostream> #include <vector/CubicIpol.hpp> #include <vector/LinearIpol.hpp> #include <vector/NearestNeighborIpol.hpp> #include <vector/MultiArray.hpp> #include <vector/Interpolate.hpp> #include <time.h> using namespace std; using namespace Fiber; /* Demonstration on how to compute the partial derivative of a three-dimensional scalar field, derived by the third index. */ static void dMA3() { double Data[13*10*23]; MultiArray<3, double> D(Data, MIndex(13,10,23) ); D.set(0.0); MultiIndex<3> I; I = 5,6,7; D[ I ] = 1.0; double value = D[ I ]; assert(value==1.0); // symbolic names for coordinate functions enum { x=0, y=1, z=2 }; for(double fr=-1.2; fr<=1.2; fr+=0.2) { //FixedArray<double,3> point; point = 5.1, 6.2, 7+fr; FixedArray<double,3> point; point = 5., 6., 7+fr; // computes derivative by coordinate z Interpolate<3, double, LinearIpol<double> , double, NoDelimiter<double>, z> LI(D, point); Interpolate<3, double, CubicIpol<double> , double, NoDelimiter<double>, z> LIc(D, point); cout << "\t" << fr << "\t df/dz " << point << " =\t" << LI.eval() << "\t|\t " << LIc.eval() << endl; } } static void printArray(const MultiArray<2, double> &V) { MultiIndex<2> idx; do { if (idx[1] == 0) cout << endl; cout << V[idx] << " "; } while( idx.inc( V.Size() ) ); cout << endl; } static void dMA2() { double Data[7*11], dData[7*11]; MultiArray<2, double> D(Data, MIndex(7,11) ), dDdx(dData, MIndex(7,11) ); D.set(0.0); D[ MIndex(3,4) ] = 1.0; D[ MIndex(4,4) ] = 0.3; // symbolic names for coordinate functions enum { x=0, y=1}; cout << "2D Array: \n"; printArray(D); cout << "Compute DeriveArray \n"; DeriveArray<LinearIpol<double>, x>::compute(dDdx, D); cout << "Derivative: \n"; printArray(dDdx); } int dInterpolCheck() { dMA3(); dMA2(); return 1; }