This code shows how numerical integrals and derivatives are taken of a given dataset. It is used as an example and a testing place for a relatively more complicated application, the PID Controller.
To get the code running, create a build folder, and go into it.
mkdir build
cd build/
Then use CMake to generate make files.
cmake ../CMakeLists.txt
And finally use make to build
make .
Then launch the executable
./numerical_calc.exe
The code starts with initializing the number of samples in the dataset, and the resolution factor. The resolution factor dictates the time delta, and therefore the accuracy of these calculations.
It then calculates the size for a ring buffer of an appropriate size.
const uint16_t sample_count = 10;
const uint16_t resolution_factor = 3000;
const double time_const = 1.0 / (double)resolution_factor;
const uint16_t buffer_size = (uint16_t)((double)sample_count / time_const);
It should be noted that specifically using a ring buffer instead of a usual buffer has no purpose in these calculations, but due to the application, it has to be used int the PID Controller project so it is being used here as well.
A function double f (double x)
is also
declared, upon which derivation and integration would be performed. In this
code example it is set to
The ring buffer is then initialized, and filled with output values of the
function double f (double x)
, with inputs starting from 0 to sample_count
.
ring_buffer rb;
ring_buffer_init(&rb, buffer_size);
for (uint16_t i = 1; i <= rb.size; i++)
{
double number = (double)i * time_const;
ring_buffer_add(&rb, f(number));
}
Derivatives are taken using the simplest method possible
This is implemented as
double calculate_derivative(const ring_buffer *rb, uint16_t index, double time_interval)
{
double derivative = 0.0;
double x_this = ring_buffer_get_item(rb, index);
double x_last = ring_buffer_get_item(rb, index - 1);
double delta = x_this - x_last;
derivative = delta / time_interval;
return derivative;
}
Integrals can be taken using the typical trapezoidal method.
This is implemented as
double trapezoidal(const ring_buffer *rb, double time_interval)
{
double integral = 0.0;
double x_first = ring_buffer_get_item(rb, 0);
double x_sum_between = 0.0;
for (uint16_t i = 1; i < rb->size; i++)
{
double x_this = ring_buffer_get_item(rb, i);
x_sum_between += x_this;
}
double x_last = ring_buffer_get_item(rb, rb->size);
integral = (time_interval / 2.0) * (x_first + (2.0 * x_sum_between) + x_last);
return integral;
}
However a more accurate integral is taken using this modified method
This is implemented as
double calculate_integral(const ring_buffer *rb, double time_interval)
{
double integral = 0.0;
for (uint16_t i = 0; i < rb->size; i++)
{
double x_this = ring_buffer_get_item(rb, i);
double x_last = ring_buffer_get_item(rb, i - 1);
integral += (x_this + x_last);
}
integral *= time_interval;
integral /= 2.0;
return integral;
}
In this code example, the input data set is 1 to 10
The function used is
Which makes the function output dataset used for calculations
The exact integral should be
When the time delta is 1 unit, the numerical integral is
When the time delta is reduced from 1 unit to 0.001 unit, the numerical integral becomes 2500.000025 and the numerical derivative at each point becomes: