17 #include "algebra/dense_matrix_operations.h"
25 CubicSpline::CubicSpline(std::vector<double> points, std::vector<double> values,
26 double bc_left,
double bc_right,
27 bool first_der_left,
bool first_der_right,
28 bool extrapolate_der_left,
bool extrapolate_der_right) :
Hermes::Hermes1DFunction<double>(), points(points), values(values),
29 bc_left(bc_left), bc_right(bc_right), first_der_left(first_der_left),
30 first_der_right(first_der_right), extrapolate_der_left(extrapolate_der_left),
31 extrapolate_der_right(extrapolate_der_right)
33 for (
unsigned short i = 1; i < points.size(); i++)
34 if (points[i] <= points[i - 1])
35 throw Exceptions::Exception(
"Points must be in ascending order when constructing a spline.");
36 this->is_const =
false;
48 void CubicSpline::free()
80 if (extrapolate_der_right ==
false)
109 else return derivative_left;
116 if (extrapolate_der_right ==
false)
return 0;
119 else return derivative_right;
127 double derivative_end,
double x_in)
const
129 return value_end + derivative_end * (x_in - point_end);
134 double x2 = x_in * x_in;
135 double x3 = x2 * x_in;
142 double x2 = x_in * x_in;
144 + 3 * this->
coeffs[m].d * x2;
150 int i_right =
points.size() - 1;
154 if (x_in <
points[i_left])
return false;
155 if (x_in >
points[i_right])
return false;
157 while (i_left + 1 < i_right)
159 int i_mid = (i_left + i_right) / 2;
160 if (
points[i_mid] < x_in) i_left = i_mid;
161 else i_right = i_mid;
168 void CubicSpline::plot(
const char* filename,
double extension,
bool plot_derivative,
int subdiv)
const
170 FILE *f = fopen(filename,
"wb");
172 throw Hermes::Exceptions::Exception(
"Could not open a spline file for writing.");
177 throw Hermes::Exceptions::Exception(
"The cubic spline has no coefficients. Calculate using calculate_coeffs.");
182 double h = extension / subdiv;
183 for (
int j = 0; j < subdiv; j++)
185 double x = x_left + j * h;
187 if (!plot_derivative) val =
value(x);
189 fprintf(f,
"%g %g\n", x, val);
193 if (!plot_derivative) val_last =
value(x_last);
195 fprintf(f,
"%g %g\n", x_last, val_last);
198 for (
unsigned int i = 0; i <
points.size() - 1; i++)
201 for (
int j = 0; j < subdiv; j++)
203 double x =
points[i] + j * h;
206 if (!plot_derivative) val = this->
value(x);
208 fprintf(f,
"%g %g\n", x, val);
213 if (!plot_derivative) val_last =
value(x_last);
215 fprintf(f,
"%g %g\n", x_last, val_last);
218 double x_right = point_right + extension;
219 h = extension / subdiv;
220 for (
int j = 0; j < subdiv; j++)
222 double x = point_right + j * h;
224 if (!plot_derivative) val =
value(x);
226 fprintf(f,
"%g %g\n", x, val);
229 if (!plot_derivative) val_last =
value(x_last);
231 fprintf(f,
"%g %g\n", x_last, val_last);
238 int nelem =
points.size() - 1;
243 this->warn(
"Empty points or values vector in CubicSpline, cancelling coefficients calculation.");
248 this->warn(
"At least two points and values required in CubicSpline, cancelling coefficients calculation.");
253 this->warn(
"Mismatched number of points and values in CubicSpline, cancelling coefficients calculation.");
258 double eps = Hermes::HermesEpsilon;
259 for (
int i = 0; i < nelem; i++)
263 this->warn(
"Duplicated or improperly ordered points in CubicSpline detected, cancelling coefficients calculation.");
271 const int n = 4 * nelem;
272 double** matrix = new_matrix<double>(n, n);
273 for (
int i = 0; i < n; i++)
275 for (
int j = 0; j < n; j++)
280 double* rhs = malloc_with_check<CubicSpline, double>(n,
this);
281 for (
int j = 0; j < n; j++)
287 for (
int i = 0; i < nelem; i++)
290 rhs[2 * i + 1] =
values[i + 1];
295 for (
int i = 0; i < nelem; i++) {
298 double xx3 = xx2 * xx;
299 matrix[2 * i][4 * i + 0] = 1;
300 matrix[2 * i][4 * i + 1] = xx;
301 matrix[2 * i][4 * i + 2] = xx2;
302 matrix[2 * i][4 * i + 3] = xx3;
306 matrix[2 * i + 1][4 * i + 0] = 1.0;
307 matrix[2 * i + 1][4 * i + 1] = xx;
308 matrix[2 * i + 1][4 * i + 2] = xx2;
309 matrix[2 * i + 1][4 * i + 3] = xx3;
314 int offset = 2 * nelem;
315 for (
int i = 1; i < nelem; i++) {
318 matrix[offset + i - 1][4 * (i - 1) + 1] = 1;
319 matrix[offset + i - 1][4 * (i - 1) + 2] = 2 * xx;
320 matrix[offset + i - 1][4 * (i - 1) + 3] = 3 * xx2;
321 matrix[offset + i - 1][4 * (i - 1) + 5] = -1;
322 matrix[offset + i - 1][4 * (i - 1) + 6] = -2 * xx;
323 matrix[offset + i - 1][4 * (i - 1) + 7] = -3 * xx2;
328 offset = 2 * nelem + nelem - 1;
329 for (
int i = 1; i < nelem; i++) {
331 matrix[offset + i - 1][4 * (i - 1) + 2] = 2;
332 matrix[offset + i - 1][4 * (i - 1) + 3] = 6 * xx;
333 matrix[offset + i - 1][4 * (i - 1) + 6] = -2;
334 matrix[offset + i - 1][4 * (i - 1) + 7] = -6 * xx;
345 offset = 2 * nelem + 2 * (nelem - 1);
350 matrix[offset + 0][2] = 2;
351 matrix[offset + 0][3] = 6 * xx;
357 matrix[offset + 0][1] = 1;
358 matrix[offset + 0][2] = 2 * xx;
359 matrix[offset + 0][3] = 3 * xx*xx;
365 if (first_der_right ==
false)
367 matrix[offset + 1][n - 2] = 2;
368 matrix[offset + 1][n - 1] = 6 * xx;
370 rhs[n - 1] = bc_right;
374 matrix[offset + 1][n - 3] = 1;
375 matrix[offset + 1][n - 2] = 2 * xx;
376 matrix[offset + 1][n - 1] = 3 * xx*xx;
378 rhs[n - 1] = bc_right;
383 int* perm = malloc_with_check<CubicSpline, int>(n,
this);
384 ludcmp(matrix, n, perm, &d);
385 lubksb<double>(matrix, n, perm, rhs);
386 free_with_check(perm);
390 for (
int i = 0; i < nelem; i++)
393 coeff.a = rhs[4 * i + 0];
394 coeff.b = rhs[4 * i + 1];
395 coeff.c = rhs[4 * i + 2];
396 coeff.d = rhs[4 * i + 3];
411 free_with_check(matrix);
412 free_with_check(rhs);
std::vector< double > points
Grid points, ordered.
std::vector< SplineCoeff > coeffs
A set of four coefficients a, b, c, d for an elementary cubic spline.
bool find_interval(double x_in, int &m) const
double get_value_from_interval(double x_in, int m) const
Gets value at a point that lies in interval 'm'.
void plot(const char *filename, double extension, bool plot_derivative=false, int subdiv=50) const
~CubicSpline()
Destructor.
double point_left
Values and derivatives at end points for extrapolation purposes.
double extrapolate_value(double point_end, double value_end, double derivative_end, double x_in) const
Extrapolate the value of the spline outside of its interval of definition.
void calculate_coeffs()
Calculates coefficients.
CubicSpline(std::vector< double > points, std::vector< double > values, double bc_left, double bc_right, bool first_der_left=true, bool first_der_right=true, bool extend_der_left=true, bool extend_der_right=true)
Constructor (general case).
bool extrapolate_der_left
double bc_left
Boundary conditions.
double derivative(double x) const
One-dimensional function derivative value.
double get_derivative_from_interval(double x_in, int m) const
Gets derivative at a point that lies in interval 'm'.
std::vector< double > values
Values at the grid points.
double value(double x) const
One-dimensional function value.