Use C++ Classes Generated from MATLAB Classes to Model Simple and Damped Oscillators

Open Live Script

This example uses a MATLAB® class that represents a simple oscillator and a subclass that adds mechanical damping to the simple oscillator.

When you use a MATLAB class to represent a physical system, you can:

Specify the system parameters by using private class properties.
Create an instance of the system by using the class constructor.
Capture the time evolution of the system by using a public method that return the trajectory of the system for a given initial state.
Modularize the mathematical analysis by creating private or protected helper methods.

When you model a physical system, you often begin with a simple system and then introduce additional effects, such as mechanical damping, to increase the accuracy of your analysis. In MATLAB, you can model the simple system as a base class and model the enhanced system as a subclass that inherits from the base class. The subclass can define private properties for additional system parameters. Additionally, the subclass can inherit certain methods from the base class and can overload the other methods.

This example shows how to generate C++ classes for MATLAB classes and how to integrate the generated C++ classes into a C++ executable.

Examine the Oscillator Base Class

A simple harmonic oscillator has two parameters, the mass $m$ and the spring constant $k$ . The angular frequency of the oscillator is $ω = \sqrt{\frac{k}{m}}$ . This equation defines the position of the oscillator $x$ as a function of time $t$ :

$x (t) = A \sin (ω t + ϕ) .$

The initial position $x_{0}$ and initial velocity $v_{0}$ determine the amplitude $A$ and the phase constant $ϕ$ .

Examine the MATLAB class simpleOscillator. This class models a one-dimensional simple harmonic oscillator in the MATLAB namespace mySystem. This class has two properties, mass and springConstant. The dynamics method returns the final position of the oscillator after a specified time interval. To evolve the system over a specified number of time steps, the evolution method calls the dynamics method iteratively.

type +mySystem/simpleOscillator.m

classdef simpleOscillator
    properties (SetAccess = private, GetAccess = protected)
        mass
        springConstant
    end

    methods
        function obj = simpleOscillator(m,k)
            obj.mass = m;
            obj.springConstant = k;
        end

        function [time,position] = evolution(obj,initialPosition,initialVelocity,timeInterval,timeStep)
            numSteps = floor(timeInterval/timeStep);
            position = zeros(numSteps + 1,1);
            time = zeros(numSteps + 1,1);
            position(1) = initialPosition;
            for i = 1:numSteps
                currentTime = i*timeStep;
                position(i+1) = obj.dynamics(initialPosition,initialVelocity,currentTime);
                time(i+1) = currentTime;
            end
        end
    end

    methods (Access = protected)
        function omega = angularFrequency(obj)
            omega = sqrt(obj.springConstant/obj.mass);
        end

        function amplitudeValue = amplitude(obj,initialPosition,initialVelocity)
            omega = obj.angularFrequency;
            positionSquared = initialPosition^2;
            velocityTerm = (initialVelocity / omega)^2;
            amplitudeValue = sqrt(positionSquared + velocityTerm);
        end

        function phi = phase(obj,initialPosition,initialVelocity)
            omega = obj.angularFrequency;
            phi = atan2(omega*initialPosition,initialVelocity);
        end

        function finalPosition = dynamics(obj,initialPosition,initialVelocity,timeInterval)
            omega = obj.angularFrequency;
            amplitudeValue = obj.amplitude(initialPosition, initialVelocity);
            phi = obj.phase(initialPosition,initialVelocity);
            finalPosition = amplitudeValue*sin(omega*timeInterval+phi);
        end
    end
end

Examine the Damped Oscillator Subclass

To model the effects of mechanical damping on a harmonic oscillator, you need one additional parameter, the damping constant $b$ . The equation $γ = \frac{b}{2 m}$ represents the damping parameter. Because the damping parameter is small compared to the angular frequency $ω$ , only the first-order damping effects are significant. The position of the damped oscillator $x_{d}$ as a function of time $t$ is:

$x_{d} (t) = A e^{- γ t} \sin (ω t + ϕ_{d})$ .

Like the simple oscillator, the initial position $x_{0}$ and initial velocity $v_{0}$ determine the amplitude $A$ and the phase constant $ϕ_{d}$ . The damping constant causes the amplitude to decay exponentially.

Examine the MATLAB class dampedOscillator. This class is a subclass of the simpleOscillator base class and is in the MATLAB namespace mySystem. The dampedOscillator subclass has one additional property compared to the simpleOscillator base class, dampingConstant. The subclass overloads the phase and dynamics methods to include the effects of damping and defines an additional method, dampingParameter, that calculates the normalized damping parameter in the dynamical equation.

type +mySystem/dampedOscillator.m

classdef dampedOscillator < mySystem.simpleOscillator
    properties (SetAccess = private, GetAccess = protected)
        dampingConstant
    end
    
    methods
        function obj = dampedOscillator(m,b,k)
            obj@mySystem.simpleOscillator(m,k);
            obj.dampingConstant = b;
        end
    end

    methods (Access = protected)        
        function gamma = dampingParameter(obj)
            gamma = obj.dampingConstant/(2*obj.mass);
        end
               
        function phi = phase(obj,initialPosition,initialVelocity)
            omega = obj.angularFrequency();
            gamma = obj.dampingParameter();
            phi = atan2(omega*initialPosition,initialVelocity+gamma*initialPosition);
        end
        
        function finalPosition = dynamics(obj,initialPosition,initialVelocity,timeInterval)
            gamma = obj.dampingParameter();
            omega = obj.angularFrequency();
            amplitudeValue = obj.amplitude(initialPosition, initialVelocity);
            phi = obj.phase(initialPosition, initialVelocity);
            expDecay = exp(-gamma * timeInterval);
            finalPosition = amplitudeValue*expDecay*sin(omega*timeInterval+phi);
        end
    end
end

Examine and Run MATLAB Function

Examine the MATLAB function effectOfDamping. This function creates a simple oscillator and a damped oscillator with identical starting conditions and calculates the trajectories of the two oscillators over time.

type effectOfDamping.m

function [time_simple,position_simple,time_damped,position_damped] = effectOfDamping(params,initialPosition,initialVelocity,timeInterval,timeStep)
% Create and evolve simple oscillator
myOscillator = mySystem.simpleOscillator(params.springConstant,params.mass);
[time_simple,position_simple] = myOscillator.evolution(initialPosition,initialVelocity,timeInterval,timeStep);

% Create and evolve damped oscillator
myDampedOscillator = mySystem.dampedOscillator(params.springConstant,params.dampingConstant,params.mass);
[time_damped,position_damped] = myDampedOscillator.evolution(initialPosition,initialVelocity,timeInterval,timeStep);
end

The effectOfDamping function takes these input arguments:

params, which is a structure with three fields that correspond to the parameters springConstant, dampingConstant, and mass in normalized units.
initialPosition, which corresponds to the initial position of the oscillator, $x_{0}$ .
initialVelocity, which corresponds to the initial velocity of the oscillator, $v_{0}$ .
timeInterval, which corresponds to the time interval over which the system evolves, $t$ .
timeStep, which corresponds to the time step that the evolution method uses to evolve the system.

The function returns four row vectors, time_simple, position_simple, time_damped, and position_damped, which describe the positions of the simple and damped oscillator at each time step.

To visualize the system in MATLAB, first create a structure to hold the system parameters.

params.springConstant = 1;
params.dampingConstant = 0.1;
params.mass = 1;

Use the effectOfDamping function to calculate the trajectories of the simple and damped oscillators. Specify an initialPosition of 1, an initialVelocity of 0, a timeInterval of 100, and a timeStep of 0.01.

[time_s,position_s,time_d,position_d] = effectOfDamping(params,1,0,100,0.01);

Plot the position of each oscillator over time. Observe that the amplitude of the damped oscillator decays exponentially with time.

plot(time_s,position_s)
hold on
plot(time_d,position_d)

Figure contains an axes object. The axes object contains 2 objects of type line.

Display the final position of the simple oscillator.

disp(position_s(end))

    0.8623

Display the final position of the damped oscillator. Damping causes the final position of this oscillator to be closer to the mean position, 0.

disp(position_d(end))

    0.0056

Generate and Run C++ MEX Function

Generate a C++ MEX function for the effectOfDamping function by using the codegen command. Then, run the generated MEX function to check that the generated code has the same behavior as the original MATLAB code.

It is a best practice to perform this step because you can run the generated MEX function to detect run-time errors that are harder to diagnose in standalone code. For example, the MEX function includes memory integrity checks by default.

By default, the codegen command generates a C MEX function in the working folder. Use the -lang:c++ option to produce a C++ MEX function. Use the -args option to specify input arguments. Use params as an example value to instruct the code generator to accept a structure with the fields springConstant, dampingConstant, and mass. Specify the remaining input arguments as scalar doubles by using 0 as an example value.

codegen -lang:c++ effectOfDamping -args {params,0,0,0,0}

Code generation successful.

Test the MEX function with the same input that you passed to the original MATLAB function and compare the results. The MEX function produces the same output.

[time_s_mex,position_s_mex,time_d_mex,position_d_mex] = effectOfDamping_mex(params,1,0,100,0.01);
disp(position_s_mex(end))

    0.8623

disp(position_d_mex(end))

    0.0056

Generate Standalone C++ Code

To generate a C++ class interface, first create a code configuration object for a static library.

cfg = coder.config("lib");

Then, modify these settings:

To generate C++ code, set cfg.TargetLang to "C++". See Language.
To instruct the code generator to generate methods in a C++ class from the effectOfDamping function, set cfg.CppInterfaceStyle to "Methods". See Interface style.
To specify a name for the interface class, set cfg.CppInterfaceClassName to "myOscillators". See C++ interface class name.
The code generator produces C++ classes from MATLAB classes by default. However, optimizations later in the code generation process can result in the inlining of these C++ classes. For this example, you can prevent the code generator from inlining the class definitions by setting cfg.InlineBetweenUserFunctions to "Never". See Inline strategy for user written MATLAB functions.

cfg.TargetLang = "C++";
cfg.CppInterfaceStyle = "Methods";
cfg.CppInterfaceClassName = "myOscillators";
cfg.InlineBetweenUserFunctions = "Never";

Generate a static C++ library by using the codegen command and specify the code configuration object by using the -config option. Use the same -args syntax that you used to generate the MEX function.

codegen -config cfg effectOfDamping -args {params,0,0,0,0} -report

Code generation successful: View report

Inspect Generated C++ Code

In the generated code, the simpleOscillator and dampedOscillator classes are defined in the mySystem C++ namespace. For more information, see Organize Generated C++ Code into Namespaces.

To inspect the generated code, open the code generation report.

C++ Class Implementations

The files simpleOscillator.cpp and dampedOscillator.cpp contain the implementations of the C++ classes for the simple and damped oscillators, respectively. Because the generated code does not preserve the inheritance structure of the MATLAB class and subclass, dampedOscillator is not a subclass of simpleOscillator. Instead, the dampedOscillator class reimplements the methods that the corresponding MATLAB class inherits. To learn more about how the code generator maps MATLAB classes to C++ classes, see Generate C++ Code for MATLAB Classes.

As an example, examine the declaration and definition of the generated mySystem::dampedOscillator class. The declaration of the mySystem::simpleOscillator class is very similar, except that it does not include the dampingParameter and b_simpleOscillator methods and the dampingConstant property.

file = fullfile("codegen","lib","effectOfDamping","dampedOscillator.h");
coder.example.extractLines(file,"namespace","};",1,1)

namespace mySystem {
class dampedOscillator {
public:
  void init(double m, double b, double k);
  void evolution(double initialPosition, double initialVelocity,
                 double timeInterval, double timeStep,
                 coder::array<double, 1U> &b_time,
                 coder::array<double, 1U> &position) const;

protected:
  double dynamics(double initialPosition, double initialVelocity,
                  double timeInterval) const;
  double dampingParameter() const;
  double angularFrequency() const;
  double amplitude(double initialPosition, double initialVelocity) const;
  double phase(double initialPosition, double initialVelocity) const;

private:
  void b_simpleOscillator(double m, double k);

protected:
  double mass;
  double springConstant;
  double dampingConstant;
};

file = fullfile("codegen","lib","effectOfDamping","dampedOscillator.cpp");
coder.example.extractLines(file,"// Arguments    : double initialPosition","} // namespace mySystem",1,1)

// Arguments    : double initialPosition
//                double initialVelocity
// Return Type  : double
//
namespace mySystem {
double dampedOscillator::amplitude(double initialPosition,
                                   double initialVelocity) const
{
  double amplitudeValue;
  amplitudeValue = initialVelocity / angularFrequency();
  return std::sqrt(initialPosition * initialPosition +
                   amplitudeValue * amplitudeValue);
}

//
// Arguments    : void
// Return Type  : double
//
double dampedOscillator::angularFrequency() const
{
  return std::sqrt(springConstant / mass);
}

//
// Arguments    : double m
//                double k
// Return Type  : void
//
void dampedOscillator::b_simpleOscillator(double m, double k)
{
  mass = m;
  springConstant = k;
}

//
// Arguments    : void
// Return Type  : double
//
double dampedOscillator::dampingParameter() const
{
  return dampingConstant / (2.0 * mass);
}

//
// Arguments    : double initialPosition
//                double initialVelocity
//                double timeInterval
// Return Type  : double
//
double dampedOscillator::dynamics(double initialPosition,
                                  double initialVelocity,
                                  double timeInterval) const
{
  return amplitude(initialPosition, initialVelocity) *
         std::exp(-dampingParameter() * timeInterval) *
         std::sin(angularFrequency() * timeInterval +
                  phase(initialPosition, initialVelocity));
}

//
// Arguments    : double initialPosition
//                double initialVelocity
// Return Type  : double
//
double dampedOscillator::phase(double initialPosition,
                               double initialVelocity) const
{
  return coder::b_atan2(angularFrequency() * initialPosition,
                        initialVelocity + dampingParameter() * initialPosition);
}

//
// Arguments    : double initialPosition
//                double initialVelocity
//                double timeInterval
//                double timeStep
//                coder::array<double, 1U> &b_time
//                coder::array<double, 1U> &position
// Return Type  : void
//
void dampedOscillator::evolution(double initialPosition, double initialVelocity,
                                 double timeInterval, double timeStep,
                                 coder::array<double, 1U> &b_time,
                                 coder::array<double, 1U> &position) const
{
  double numSteps;
  int loop_ub;
  numSteps = std::floor(timeInterval / timeStep);
  loop_ub = static_cast<int>(numSteps + 1.0);
  position.set_size(static_cast<int>(numSteps + 1.0));
  for (int i{0}; i < loop_ub; i++) {
    position[i] = 0.0;
  }
  b_time.set_size(static_cast<int>(numSteps + 1.0));
  for (int i{0}; i < loop_ub; i++) {
    b_time[i] = 0.0;
  }
  position[0] = initialPosition;
  loop_ub = static_cast<int>(numSteps);
  for (int i{0}; i < loop_ub; i++) {
    numSteps = (static_cast<double>(i) + 1.0) * timeStep;
    position[i + 1] = dynamics(initialPosition, initialVelocity, numSteps);
    b_time[i + 1] = numSteps;
  }
}

//
// Arguments    : double m
//                double b
//                double k
// Return Type  : void
//
void dampedOscillator::init(double m, double b, double k)
{
  b_simpleOscillator(m, k);
  dampingConstant = b;
}

} // namespace mySystem

Interface Class Implementation

The files myOscillators.h and myOscillators.cpp contain the implementation of the interface class myOscillators. The method myOscillators::effectOfDamping implements the entry-point function effectOfDamping. Examine the declaration of this method in the header file myOscillators.h.

file = fullfile("codegen","lib","effectOfDamping","myOscillators.h");
coder.example.extractLines(file,"class","#endif",1,0)

class myOscillators {
public:
  myOscillators();
  ~myOscillators();
  void effectOfDamping(const struct0_T *params, double initialPosition,
                       double initialVelocity, double timeInterval,
                       double timeStep, coder::array<double, 1U> &time_simple,
                       coder::array<double, 1U> &position_simple,
                       coder::array<double, 1U> &time_damped,
                       coder::array<double, 1U> &position_damped);
};

The class constructor and the class destructor implement the initialize and terminate functions, respectively. The run-time inputs timeInterval and timeStep determine the size of the output arguments of the effectOfDamping function. The generated code represents these arguments as dynamic C++ arrays by using the coder::array class template. For more information, see Use Dynamically Allocated C++ Arrays in Generated Function Interfaces.

Example `main` Function

Examine the example C++ main function. When you generate standalone code, the code generator produces an example C or C++ main function. Because the example main function demonstrates how to call the generated C or C++ function, you can use it as a template for your application.

The code generator creates the example CPP and H files in the folder codegen/lib/effectOfDamping/examples. The code generator overwrites the files in the examples folder every time you generate code. Copy the example main source and header files to a location outside of the codegen folder before you modify these files.

Examine the example C++ main function.

file = fullfile("codegen","lib","effectOfDamping","examples","main.cpp");
coder.example.extractLines(file,"int main","return 0;",1,1)

int main(int, char **)
{
  myOscillators *classInstance;
  classInstance = new myOscillators;
  // Invoke the entry-point functions.
  // You can call entry-point functions multiple times.
  main_effectOfDamping(classInstance);
  delete classInstance;
  return 0;

The example main function uses the C++ new operator to allocate memory for an instance of myOscillators, invokes the main_effectOfDamping function, and frees the memory by using the C++ delete operator.

Examine the example main_effectOfDamping function.

file = fullfile("codegen","lib","effectOfDamping","examples","main.cpp");
coder.example.extractLines(file,"void main_effectOfDamping","// File trailer",1,0)

void main_effectOfDamping(myOscillators *instancePtr)
{
  coder::array<double, 1U> position_damped;
  coder::array<double, 1U> position_simple;
  coder::array<double, 1U> time_damped;
  coder::array<double, 1U> time_simple;
  struct0_T r;
  double initialPosition_tmp;
  // Initialize function 'effectOfDamping' input arguments.
  // Initialize function input argument 'params'.
  initialPosition_tmp = argInit_real_T();
  // Call the entry-point 'effectOfDamping'.
  r = argInit_struct0_T();
  instancePtr->effectOfDamping(&r, initialPosition_tmp, initialPosition_tmp,
                               initialPosition_tmp, initialPosition_tmp,
                               time_simple, position_simple, time_damped,
                               position_damped);
}

//

The size of output arguments of the main_effectOfDamping function are determined by the run-time inputs timeInterval and timeStep. So, the generated code represents these arguments as dynamic C++ arrays that are implemented by using the coder::array class template. For more information, see Use Dynamically Allocated C++ Arrays in Generated Function Interfaces.

Examine Modified `main` Function

For this example, the files main_damped_oscillator.h and main_damped_oscillator.cpp in the working directory are edited versions of the example files that have been modified to interact with the generated code. Examine the main function in the file main_damped_oscillator.cpp.

coder.example.extractLines("main_damped_oscillator.cpp","int main","",1,0)

int main(int, const char * const [])
{
  myOscillators classInstance;

  // Invoke the entry-point functions.
  // You can call entry-point functions multiple times.
  main_effectOfDamping(&classInstance);
  return 0;
}

The main function uses the interface class myOscillators to interact with the generated code.

Examine the C++ function main_effectOfDamping.

coder.example.extractLines("main_damped_oscillator.cpp","// Return Type  : void","}",0,1)

//
static void main_effectOfDamping(myOscillators *instancePtr)
{
  coder::array<double, 1U> position1;
  coder::array<double, 1U> position2;
  coder::array<double, 1U> time1;
  coder::array<double, 1U> time2;
  struct0_T r;

  // Initialize function 'effectOfDamping' input arguments.
  // Initialize function input argument 'params'.
  r.springConstant = 1;
  r.dampingConstant = 0.1;
  r.mass = 1;

  // Call the entry-point 'effectOfDamping'.
  instancePtr->effectOfDamping((&r), 1, 0, 100, 0.01, time1, position1, time2, position2);

  std::cout << position1[position1.size(0) - 1] << std::endl;
  std::cout << position2[position2.size(0) - 1] << std::endl;
}

The main_effectOfDamping function performs the same computation as the MATLAB function effectOfDamping. The C++ function uses the coder::array API to interact with the dynamic arrays that the generated effectOfDamping function returns. At the end of its execution, the main_effectOfDamping function prints the final positions of the two oscillators.

Generate and Test C++ Executable

To generate a C++ executable for the effectOfDamping function, modify the code configuration object cfg:

To instruct the code generator to generate an executable file, set the cfg object property OutputType to "EXE". See Build type.
To instruct to code generator to compile the source file main_damped_oscillator.cpp and link it with the generated code, set the cfg object property CustomSource to "main_damped_oscillator.cpp". See Additional source files.

cfg.OutputType = "EXE";
cfg.CustomSource = "main_damped_oscillator.cpp";

Generate an executable by using the codegen command and specify the code configuration object by using the -config option. Use the same -args syntax that you used to generate the MEX function. The code generator creates the effectOfDamping C++ executable in the working folder.

codegen -config cfg effectOfDamping -args {params,0,0,0,0}

Code generation successful.

To run the generated executable, use the system command. Use the ispc and isunix functions to select the appropriate system command. The output of the C++ executable is the same as the output of the original MATLAB function.

if isunix
    system('./effectOfDamping');
elseif ispc
    system('effectOfDamping.exe');
else
    disp('Platform is not supported');
end

0.862319
0.00563263

Use C++ Classes Generated from MATLAB Classes to Model Simple and Damped Oscillators

Examine the Oscillator Base Class

Examine the Damped Oscillator Subclass

Examine and Run MATLAB Function

Generate and Run C++ MEX Function

Generate Standalone C++ Code

Inspect Generated C++ Code

C++ Class Implementations

Interface Class Implementation

Example `main` Function

Examine Modified `main` Function

Generate and Test C++ Executable

See Also

Topics

Use C++ Classes Generated from MATLAB Classes to Model Simple and Damped Oscillators

Examine the Oscillator Base Class

Examine the Damped Oscillator Subclass

Examine and Run MATLAB Function

Generate and Run C++ MEX Function

Generate Standalone C++ Code

Inspect Generated C++ Code

C++ Class Implementations

Interface Class Implementation

Example main Function

Examine Modified main Function

Generate and Test C++ Executable

See Also

Topics

Example `main` Function

Examine Modified `main` Function