include iostream include vector using namespace std typedef vector dou

  1
  2
  3
  4
  5
  6
  7
  8
  9
 10
 11
 12
 13
 14
 15
 16
 17
 18
 19
 20
 21
 22
 23
 24
 25
 26
 27
 28
 29
 30
 31
 32
 33
 34
 35
 36
 37
 38
 39
 40
 41
 42
 43
 44
 45
 46
 47
 48
 49
 50
 51
 52
 53
 54
 55
 56
 57
 58
 59
 60
 61
 62
 63
 64
 65
 66
 67
 68
 69
 70
 71
 72
 73
 74
 75
 76
 77
 78
 79
 80
 81
 82
 83
 84
 85
 86
 87
 88
 89
 90
 91
 92
 93
 94
 95
 96
 97
 98
 99
100
101
102
103
104
105
106
107
108
#include <iostream>
#include <vector>
using namespace std;
typedef vector<double> Vec;
typedef vector<Vec> Matrix;
int n;
void BuildSpline(Vec knotX, Vec knotY, Matrix &Coef)
{
Vec a(n);
Vec b(n);
Vec c(n);
Vec d(n);
Vec alpha(n - 1);
Vec betta(n - 1);
for (int i = 0; i < n; ++i)
{
a[i] = knotY[i];
}
double h1, h2, A, B, C, F, temp;
alpha[0] = betta[0] = 0;
for (int i = 1; i < n - 1; ++i)
{
h1 = knotX[i] - knotX[i - 1];
h2 = knotX[i + 1] - knotX[i];
A = h1;
C = 2 * (h1 + h2);
B = h2;
F = 6 * ((knotY[i + 1] - knotY[i]) / h2 - (knotY[i] - knotY[i - 1]) / h1);
temp = (A * alpha[i - 1] + C);
alpha[i] = -B / temp;
betta[i] = (F - A * betta[i - 1]) / temp;
}
c[n - 1] = (F - A * betta[n - 2]) / (C + A * alpha[n - 2]);
for (int i = n - 2; i > 0; --i)
{
c[i] = alpha[i] * c[i + 1] + betta[i];
}
for (int i = n - 1; i > 0; --i)
{
double h = knotX[i] - knotX[i - 1];
d[i] = (c[i] - c[i-1]) /h;
b[i] = h * (2 * c[i] + c[i-1]) / 6.0 + (knotY[i] - knotY[i - 1]) / h;
}
for (int j = 0; j< n; j++)
{
Coef[j][0] = a[j];
Coef[j][1] = b[j];
Coef[j][2] = c[j];
Coef[j][3] = d[j];
}
}
double Interpolate(double x, Vec knotX, Matrix Coef)
{
int ind = 0;
if (x <= knotX[0])
ind = 1;
else if (x >= knotX[n - 1])
ind = n-1;
else
{
int i = 0, j = n - 1;
while (i + 1 < j)
{
std::size_t k = i + (j - i) / 2;
if (x <= knotX[k])
j = k;
else
i = k;
}
ind = j;
}
double dx = x - knotX[ind];
return Coef[ind][0] + (Coef[ind][1] + (Coef[ind][2] / 2 + Coef[ind][3] * dx / 6)*dx)*dx;
}
int main()
{
freopen("output.txt", "w", stdout);
Vec knotX = { -3.14, -2.44, -1.74, -1.04,-0.34, 0.36, 1.06, 1.76, 2.46};
Vec knotY = { -0.0016, -0.6454, -0.9857, -0.8624, -0.3335, 0.3523, 0.8724, 0.9822, 0.6300 };
n = knotX.size();
Matrix Coef(n, Vec(4, 0));
BuildSpline(knotX, knotY, Coef);
cout << "[ ";
for (double x = -3.14; x <= 3.14; x += 0.1)
{
cout << x << " ";
double temp = Interpolate(x, knotX, Coef);
cout << temp << endl;
}
cout << " ]";
return 0;
}