InterpolatedDerivative.cpp

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;
}

---