Старожил
Сообщения: 299
Благодарности: 48
|
Профиль
| Цитировать
Дан набор точек:
читать дальше »
Код:  x: y:
50 2.4
55 1.3
60 73
65 40
70 21
75 950
80 430
85 190
90 7600
95 3200
100 1400
Нужно интерполировать функцию, используя кубический сплайн.
На Википедии нашёл формулы для вычисления коэффициентов сплайна (аж в двух вариантах — один на английской Википедии, другой на русской).
Вот, что получилось (оба исходника написаны на C++ с использованием библиотеки OpenBGI).
Программа с коэффициентами взятыми из русской Википедии (почти полностью списана с образца, приведённого там же):
Код: #include <fstream>
#include <cmath>
#include "graphics.h"
using namespace std;
ifstream input("input.txt");
// Display parameters.
const int Width = 800;
const int Height = 600;
const int Border = 10;
// X increment on the screen.
const int dx_display = 3;
class Spline
{
// Spline coefficients at every segment.
struct Spline_chunk {
double a, b, c, d, x;
};
int n; // Number of segments.
Spline_chunk *spline; // Full spline.
double min, max; // Maximum and minimum spline value.
void Free_memory(void)
{
delete[] spline;
spline = NULL;
}
public:
// Constructor.
Spline(const double *x, const double *y, const int n) {
Create_spline(x, y, n);
}
// Destructor.
~Spline(void) {
Free_memory();
}
// Spline creation.
// x — function's arguments, y — function's value,
// n — number of segments.
void Create_spline(const double *x, const double *y, const int n);
// Spline function.
double F(double x);
// Minimum argument's value.
double MinX(void) {
return spline[0].x;
}
// Maximum argument's value.
double MaxX(void) {
return spline[n - 1].x;
}
// Minimum function's value.
double MinY(void) {
return min;
}
// Maximum function's value.
double MaxY(void) {
return max;
}
};
int main()
{
const int n = 11; // Count of points.
// Function's values.
double x[n], y[n];
// Input function's values.
for(int i = 0; i < n ; i++)
input >> x[i] >> y[i];
// Create spline.
Spline spline(x, y, n);
// Scale factors.
double ScaleX = (Width - 2 * Border - 1) /
(spline.MaxX() - spline.MinX());
double ScaleY = (Height - 2 * Border - 1) /
(spline.MaxY() - spline.MinY());
// Set X increment.
double dx = dx_display / ScaleX;
// Graphics initialization.
int gd = CUSTOM, gm = CUSTOM_MODE(Width, Height);
initgraph(&gd, &gm, "");
setcolor(RED);
// Move to first point.
moveto(Border, Height - 1 - Border -
ScaleY * (spline.MaxY() - spline.F(spline.MinX())));
// Draw function.
for(double x = spline.MinX() + dx; x <= spline.MaxX(); x += dx)
lineto(Border + ScaleX * (x - spline.MinX()), Height - Border -
1 - ScaleY * (spline.MaxY() - spline.F(x)));
readkey();
closegraph();
return 0;
}
void Spline::Create_spline(const double *x, const double *y, const int n)
{
this->n = n;
spline = new Spline_chunk[n];
int i; // Counter.
// Initialazing a and x coefficients.
for(i = 0; i < n; i++) {
spline[i].x = x[i];
spline[i].a = y[i];
}
// Calculating c.
spline[0].c = spline[n - 1].c = 0;
// Shuttle coefficients.
double *alpha = new double[n - 1];
double *beta = new double[n - 1];
alpha[0] = beta[0] = 0;
// Calculating shuttle coefficients.
for(i = 1; i < n - 1; i++) {
double h_i = x[i] - x[i - 1], h_i1 = x[i + 1] - x[i];
double A = h_i, C = 2 * (h_i + h_i1), B = h_i1;
double z = A * alpha[i - 1] + C;
alpha[i] = -B / z;
beta[i] = (6 * ((y[i + 1] - y[i]) / h_i1 - (y[i] - y[i - 1]) /
h_i) - A * beta[i - 1]) / z;
}
// Finding solution.
for(i = n - 2; i > 0; i--)
spline[i].c = alpha[i] * spline[i + 1].c + beta[i];
delete[] alpha; delete[] beta;
// Calculating b and d.
for(i = 1; i < n; i++) {
double h_i = x[i] - x[i - 1];
spline[i].b = (y[i] - y[i - 1]) / h_i - h_i * (4 *
spline[i].c - spline[i - 1].c) / 6;
spline[i].d = (spline[i].c - spline[i - 1].c) / h_i;
}
// Coefficients for first spline.
spline[0].b = spline[1].b - spline[1].c * (spline[1].x - spline[0].x) *
(spline[1].x - spline[0].x);
spline[0].d = spline[1].c;
/******************************************************
* Searching minimum and maximum values of function. *
******************************************************/
// Initializing values
// (spline[].a is a function value on boundary point).
min = max = spline[0].a;
// Find minimum and maximum values of function at the boundary points.
for(Spline_chunk *s = spline; s < spline + n; s++) {
if(s->a < min)
min = s->a;
else if(s->a > max)
max = s->a;
}
// Finding extremum on each segment by using derivative.
for(Spline_chunk *s = spline; s < spline + n; s++) {
// Discriminant of the function's derivative.
double D = (s->c - s->d * s->x) * (s->c - s->d * s->x) -
2. * s->d * (s->b + (s->d * s->x / 2. - s->c) *
s->x);
if(D == 0) {
double dx = -s->c / s->d; // Argument's increment.
// Calculating function's value at this point.
double y = s->a + (s->b + (s->c / 2. + s->d *
dx / 6.) * dx) * dx;
// Compare with avaliable values.
if(y < min)
min = y;
else if(y > max)
max = y;
}
else if(D > 0) {
// Argument's increments.
double dx1 = (sqrt(D) - s->c) / s->d;
double dx2 = -(sqrt(D) + s->c) / s->d;
// Calculating function's values at this points.
double y1 = s->a + (s->b + (s->c / 2. + s->d *
dx1 / 6.) * dx1) * dx1;
double y2 = s->a + (s->b + (s->c / 2. + s->d *
dx2 / 6.) * dx2) * dx2;
// Sorting these values:
// y1 — minimum, y2 — maximum from them.
if(y1 > y2) {
y1 = y1 + y2;
y2 = y1 - y2;
y1 = y1 - y2;
}
// Compare with avaliable values.
if(y1 < min)
min = y1;
if(y2 > max)
max = y2;
}
}
/*********************
* End of searching. *
*********************/
}
double Spline::F(double x)
{
// Point to a corresponding spline segment.
Spline_chunk *s;
// If x is less than or equal to first array element,
// use first spline chunk.
if(x <= spline[0].x)
s = spline;
// If x is greater than or equal to last array element,
// use last spline chunk.
else if(x >= spline[n - 1].x)
s = spline + (n - 1);
// If x lies between boundary points,
// finding corresponding segment by using binary search.
else {
int a = 0, b = n - 1;
while(b - a > 1) {
int i = (a + b) / 2;
if(x <= spline[i].x)
b = i;
else
a = i;
}
s = spline + a;
}
double dx = x - s->x;
// Return function value.
return s->a + (s->b + (s->c / 2. + s->d * dx / 6.) * dx) * dx;
}
Программа с коэффициентами взятыми из английской Википедии:
Код: #include <stdio.h>
#include "graphics.h"
FILE *input; // Input stream.
// Display parameters.
const int Width = 800;
const int Height = 600;
const int Border = 10;
// X increment on the screen.
const int dx_display = 3;
const int n = 10; // Count of spline segments.
// Struct containing spline.
struct{
double a, b, c, d, x;
}spline[n];
// Function returning spline value at point x.
double F(double x)
{
int a = 0, b = n - 1, i; // Variables using in binary search.
int s; // Number of corresponding spline segment.
// If x is less than second array element,
// use first spline chunk.
if(x < spline[1].x)
s = 0;
// If x is greater than or equal to last array element,
// use last spline chunk.
else if(x >= spline[n - 1].x)
s = n - 1;
// If x lies between boundary points,
// finding corresponding segment by using binary search.
else {
while(b - a > 1) {
i = (a + b) / 2;
if(x <= spline[i].x)
b = i;
else
a = i;
}
s = a;
}
double dx = x - spline[s].x;
// Return function value.
return spline[s].a + spline[s].b * dx + spline[s].c * dx * dx +
spline[s].d * dx * dx * dx;
}
int main()
{
// Function's values: x — arguments, a — values.
double X[n + 1], a[n + 1];
double b[n], c[n + 1], d[n]; // Spline coefficients.
double h[n]; // x difference.
// Auxiliary coefficients.
double alpha[n - 1], l[n + 1], mu[n + 1], z[n + 1];
// Minimum and maximum spline values.
double Min, Max;
// Scale factors for drawing.
double ScaleX, ScaleY;
// x — argument, dx — increment used by some functions.
double x, dx;
// Open stream.
input = fopen("input.txt", "rt");
// Input function's values.
for(int i = 0; i < n + 1; i++)
fscanf(input, "%lf%lf", &X[i], &a[i]);
// Calculating differences.
for(i = 0; i < n; i++)
h[i] = X[i + 1] - X[i];
// Calculating alpha coefficients.
alpha[0] = 0;
for(i = 1; i < n - 1; i++)
alpha[i] = 3 / h[i] * (a[i + 1] - a[i]) - 3 / h[i - 1] *
(a[i] - a[i - 1]);
// Calculating other auxiliary coefficients.
l[0] = 1; mu[0] = z[0] = 0;
for(i = 1; i < n - 1; i++) {
l[i] = 2 * (X[i + 1] - X[i - 1]) - h[i - 1] * mu[i - 1];
mu[i] = h[i] / l[i];
z[i] = (alpha[i] - h[i - 1] * z[i - 1]) / l[i];
}
l[n] = 1; z[n] = c[n] = 0;
// Calculating spline coefficients.
for(i = n - 1; i >= 0; i--) {
c[i] = z[i] - mu[i] * c[i + 1];
b[i] = (a[i + 1] - a[i]) / h[i] - h[i] *
(c[i + 1] + 2 * c[i]) / 3;
d[i] = (c[i + 1] - c[i]) / 3 * h[i];
}
// Populates them into spline structure.
for(int i = 0; i < n; i++) {
spline[i].a = a[i];
spline[i].b = b[i];
spline[i].c = c[i];
spline[i].d = d[i];
spline[i].x = X[i];
}
// Set scale for x axis.
ScaleX = (Width - 1 - 2 * Border) / (spline[n - 1].x - spline[0].x);
// Set x increment.
dx = dx_display / ScaleX;
// Finding minimum and maxim spline values.
Min = Max = spline[0].a;
for(x = spline[0].x + dx; x < spline[n - 1].x; x+=dx) {
// Calculating function value at this point.
ScaleY = F(x);
if(ScaleY < Min)
Min = ScaleY;
else if(ScaleY > Max)
Max = ScaleY;
}
// Scale factors.
ScaleY = (Height - 1 - 2 * Border) / (Max - Min);
// Graphics initialization.
int gd = CUSTOM, gm = CUSTOM_MODE(Width, Height);
initgraph(&gd, &gm, "");
setcolor(RED);
// Move to the first point.
moveto(Border, Height - 1 - Border - ScaleY * (Max - F(spline[0].x)));
// Draw function.
for(x = spline[0].x + dx; x <= spline[n - 1].x; x += dx)
lineto(Border + ScaleX * (x - spline[0].x), Height - 1 -
Border - ScaleY * (Max - F(x)));
readkey();
closegraph();
fclose(input);
return 0;
}
Правильный вариант сплайна, построенный с помощью пакета Mathematica:
Сплайны получаются разные, при этом в моей реализации сплайн даже не является гладкой кривой. Не могу понять, где ошибся, разве что формулы для коэффициентов не верны. Не мог бы кто-нибудь помочь разобраться с задачей. Заранее спасибо. |
