91 lines
2.1 KiB
C
91 lines
2.1 KiB
C
#include "integrate.h"
|
|
|
|
static inline int bad_args(integrand_f f, float *a, float *b, unsigned *n,
|
|
float *out)
|
|
{
|
|
return (f == NULL || out == NULL || *n == 0u || !(*b > *a));
|
|
}
|
|
|
|
integ_status_t midpoint(integrand_f f, void *ctx, float a, float b, unsigned n,
|
|
float *out)
|
|
{
|
|
/*check if all necessary parameters have been given */
|
|
if (bad_args(f, &a, &b, &n, out))
|
|
return INTEG_ERR_BAD_ARGS;
|
|
|
|
float sum = 0.0f;
|
|
float fx = 0.0f;
|
|
|
|
/*intervall width h: */
|
|
const float h = (b - a) / (float)n;
|
|
float x = a + 0.5f * h;
|
|
|
|
for (unsigned i = 0; i < n; ++i)
|
|
{
|
|
fx = f(x, ctx);
|
|
sum = sum + fx;
|
|
x = x + h;
|
|
}
|
|
|
|
*out = sum * h;
|
|
return INTEG_OK;
|
|
}
|
|
|
|
integ_status_t trapezoid(integrand_f f, void *ctx, float a, float b, unsigned n,
|
|
float *out)
|
|
{
|
|
/*check if all necessary parameters have been given */
|
|
if (bad_args(f, &a, &b, &n, out))
|
|
{
|
|
return INTEG_ERR_BAD_ARGS;
|
|
}
|
|
|
|
/*intervall width h: */
|
|
const float h = (b - a) / (float)n;
|
|
float sum = 0.5f * (f(a, ctx) + f(b, ctx));
|
|
|
|
float x = a + h;
|
|
for (unsigned i = 1; i < n; ++i)
|
|
{
|
|
sum = sum + f(x, ctx);
|
|
x = x + h;
|
|
}
|
|
|
|
*out = sum * h;
|
|
return INTEG_OK;
|
|
}
|
|
|
|
integ_status_t simpson(integrand_f f, void *ctx, float a, float b, unsigned n,
|
|
float *out)
|
|
{
|
|
/*check if all necessary parameters have been given */
|
|
if (bad_args(f, &a, &b, &n, out))
|
|
{
|
|
return INTEG_ERR_BAD_ARGS;
|
|
}
|
|
if (n & 1u) /*n must be even*/
|
|
{
|
|
return INTEG_ERR_N_EVEN_REQUIRED;
|
|
}
|
|
|
|
float sum_odd = 0.0f;
|
|
float sum_even = 0.0f;
|
|
/*intervall width h: */
|
|
const float h = (b - a) / (float)n;
|
|
|
|
float x = a + h;
|
|
for (unsigned i = 1; i < n; ++i)
|
|
{
|
|
const float fx = f(x, ctx);
|
|
if (i & 1u)
|
|
sum_odd = sum_odd + fx;
|
|
else
|
|
sum_even = sum_even + fx;
|
|
x = x + h;
|
|
}
|
|
|
|
*out =
|
|
(h / 3.0f) * (f(a, ctx) + 4.0f * sum_odd + 2.0f * sum_even + f(b, ctx));
|
|
return INTEG_OK;
|
|
}
|