analytical_spring.cc

This example shows how to use the ARGUMENT, TARGET, ERROR_LOG and THIS macro definitions to implement a multi-valued inverse function with private data in an atomic class.

In particular, this atomic class implements the following callbacks:

• CONSTRUCT
• DESTRUCT

Consult the following website for a detailed explanation of how the equations implemented in this class have been derived: http://oregonstate.edu/dept/math/CalculusQuestStudyGuides/ode/second/so_lin_homocc/so_lin_homocc.html

/// analytical_spring.cc
/// Atomic class that computes the analytical solution (a.k.a. closed-form solution)
/// of the second-order, constant coefficient, linear, homogeneous ODE.
///
/// Second order differential equation modeling simple harmonic motion,
/// as might be described by a vibrating spring, whose motion is resisted by 
/// a force proportional to velocity.
///
/// ODE:
///    ddy + c*dy + k*y = 0
///
/// where :
///    y   : displacement of a unit mass attached to the end of the spring
///    dy  : first derivative of y
///    ddy : second derivative of y
///
///    c   : damping constant
///    k   : stiffness constant
///
///////////////////////////////////////////////////////////////////////////////

#ifdef SPARK_PARSER

PORT time "time" [s];

PORT k "stiffness constant" [1/s^2];
PORT c "damping constant" [1/s];

// Note it would be possible to use any other distance unit in place of [m]
// as long as the same unit is consistently used for y, dy, and ddy
PORT y   "displacement" [m];
PORT dy  "velocity" [m/s];
PORT ddy "acceleration" [m/s^2];

PORT y_IC  "displacement (initial condition)" [m];
PORT dy_IC "velocity (initial condition)" [m/s];


EQUATIONS {
    y, dy, ddy = analytical_spring__evaluate( k, c );
}

FUNCTIONS {
    y, dy, ddy = analytical_spring__evaluate( time )
        CONSTRUCT = analytical_spring__construct( time, y_IC, dy_IC, k, c )
        DESTRUCT  = analytical_spring__destruct( )
        ;
}

#endif  /*SPARK_PARSER*/ 

// System includes
#include <iostream>
using std::endl;

#include <strstream>
using std::ostrstream;
using std::ends;

#include <string>
using std::string;


#include "spark.h"


///////////////////////////////////////////////////////////////////////////////
namespace analytical_spring {

///////////////////////////////////////////////////////////////////////////
/// Class that computes the closed-form solution of the linear, 2nd order ODE
/// assuming that the initial conditions passed to the constructor are
/// consistent and that the constants k and c remain unchanged throughout the 
/// simulation. If the constants k and c change, then you should construct a 
/// anew instance of the class TData with the new values passed to the constructor.
///
/// General form of ths solution when characterictic polynomial has 2 distinct 
/// real roots ALPHA and BETA:
///    y(t) = C1*exp(ALPHA*t) + C2*exp(BETA*t)
///
/// Initial conditions:        ddy(t0) + c*dy(t0) + k*y(t0) = 0
/// Characteristic polynomial: c^2-4k > 0
///
/// Distinct real roots (both negative):
///    ALPHA = 0.5*( -c + sqrt(c^2-4k) )
///    BETA  = 0.5*( -c - sqrt(c^2-4k) )
///
class TData {
public:
    ///////////////////////////////////////////////////////////////////////
    // Structors
    ///////////////////////////////////////////////////////////////////////
    TData(double time0, double y0, double dy0, double k, double c);
    ~TData();

    ///////////////////////////////////////////////////////////////////////
    // Access methods
    ///////////////////////////////////////////////////////////////////////
    double GetALPHA() const { return ALPHA; }
    double GetBETA() const { return BETA; }

    double GetC1() const { return C1; }
    double GetC2() const { return C2; }

    double y(double time) const ;
    double dy(double time) const ;
    double ddy(double time) const ;

private:
    double ALPHA;
    double BETA;
    double C1;
    double C2;

};
///////////////////////////////////////////////////////////////////////////


}; // namespace analytical_spring 
///////////////////////////////////////////////////////////////////////////////



///////////////
// CALLBACKS //
///////////////


CONSTRUCT( analytical_spring__construct )
{
    ARGUMENT( 0, time );
    ARGUMENT( 1, y_IC );
    ARGUMENT( 2, dy_IC );
    ARGUMENT( 3, k );
    ARGUMENT( 4, c );

    analytical_spring::TData* MyData = 0;

    try {
        MyData = new analytical_spring::TData(
            time,
            y_IC.GetInit(),
            dy_IC.GetInit(),
            k,
            c
        );
    }

    // Process exception thrown by constructor if any
    catch (const string& Msg) { 
        REQUEST__ABORT( Msg.c_str() );
    }


    // Abort simulation if memory could not be allocated!
    if ( !MyData ) { 
        REQUEST__ABORT( "Could not allocate memory for my private data!" );
    }


    // Store pointer to private data within this object
    THIS->SetData( MyData );


    // Report equation string for analytical solution to run log file
    ostrstream Text;
    Text << "y(t) = (" << MyData->GetC1() << " * exp(" << MyData->GetALPHA() 
         << "*t)) + (" << MyData->GetC2() << " * exp(" << MyData->GetBETA() 
         << "*t))" << ends;

    RUN_LOG( Text.str() );
}


EVALUATE( analytical_spring__evaluate )
{
    ARGUMENT( 0, time );
    TARGET( 0, y );
    TARGET( 1, dy );
    TARGET( 2, ddy );

    // Cast void* to private data type
    analytical_spring::TData* MyData = static_cast<analytical_spring::TData*>( 
        THIS->GetData()
    );

    // Set target values
    y   = MyData->y( time );
    dy  = MyData->dy( time );
    ddy = MyData->ddy( time );
}


DESTRUCT( analytical_spring__destruct )
{
    // Cast void* to private data type
    analytical_spring::TData* MyData = static_cast<analytical_spring::TData*>(
        THIS->GetData()
    );

    // Release allocated memory
    if ( MyData ) {
        delete MyData;
    }
}



//////////////////
// PRIVATE DATA //
//////////////////

using analytical_spring::TData;


// Function name        : TData::TData
// Description      : 
// Return type          : 
// Argument         : double time0
// Argument         : double y0
// Argument         : double dy0
// Argument         : double k
// Argument         : double c
TData::TData(double time0, double y0, double dy0, double k, double c)
{
    double Discriminant = c*c - 4.0*k;

    if ( Discriminant > 0.0 ) { // strongly damped spring
        //
        // y(t) = C1*exp(ALPHA*t) + C2*exp(BETA*t)
        //
        double Temp = sqrt(Discriminant);
        ALPHA = 0.5*(-c + Temp);
        BETA  = 0.5*(-c - Temp);

        // Initial conditions
        //
        // C1*exp(ALPHA*time0) + C2*exp(BETA*time0) = y(time0) = y0
        // C1*ALPHA*exp(ALPHA*time0) + C2*BETA*exp(BETA*time0) = dy(time0) = dy0
        //
        const double Determinant = BETA - ALPHA;
        const double A_11 = exp(ALPHA*time0);
        const double A_12 = exp(BETA*time0);
        const double A_21 = ALPHA*exp(ALPHA*time0);
        const double A_22 = BETA*exp(BETA*time0);

        C1 = 1.0/Determinant * ( A_22*y0 - A_12*dy0 );
        C2 = 1.0/Determinant * (-A_21*y0 + A_11*dy0 );
    }
    else {
        string ErrMsg("Could not construct analytical_spring::TData object: ");
        ErrMsg += "characteristic polynomial does not have 2 distinct real roots!";
        throw ErrMsg;
    }
}


// Function name        : TData::~TData
// Description      : 
// Return type          : 
TData::~TData()
{
    // Nothing to do
}


// Function name        : TData::y
// Description      : 
// Return type          : double 
// Argument         : double time
double TData::y(double time) const 
{
    return (
        C1 * exp(ALPHA * time) + 
        C2 * exp(BETA * time)
    );
}


// Function name        : TData::dy
// Description      : 
// Return type          : double 
// Argument         : double time
double TData::dy(double time) const 
{
    return (
        C1 * ALPHA * exp(ALPHA * time) + 
        C2 * BETA * exp(BETA * time)
    );
}


// Function name        : TData::ddy
// Description      : 
// Return type          : double 
// Argument         : double time
double TData::ddy(double time) const 
{
    return (
        C1 * ALPHA * ALPHA * exp(ALPHA * time) + 
        C2 * BETA * BETA * exp(BETA * time)
    );
}

---