#include "stdafx.h"
#include <cmath>
#include <iostream>
#include <fstream>
#include <iomanip>
using namespace std;
double a = 20;
double b = 30;
const int N = 5;
double A0 = 2;
double B0 = 8;
int flag = 1; //1 - равномерная сетка
//2 - неравномерная сетка
double eps = 1e-3;
double func(double x)
{
if (x > 24) return 100;
else return 101;
}
double X[N + 1];
double Y[N + 1];
double H[N];
double C[N + 1];
void make_spline()
{
C[0] = A0;
C[N] = B0;
if (flag == 1)
{
for (int i = 0; i < N; i++)
{
X[i] = a + (b - a) / N * i;
}
X[N] = b;
}
else
{
X[0] = a;
X[N] = b;
for (int i = 1; i <= N - 1; i++)
{
X[i] = a + (b - a) / N * ((i + 1) + 0.4*cos(5 * (i + 1)) - 1);
}
}
for (int i = 0; i <= N; i++)
{
Y[i] = func(X[i]);
}
for (int i = 0; i < N; i++)
{
H[i] = X[i + 1] - X[i];
}
double alphas[N];
double betas[N];
alphas[N - 1] = 0;
betas[N - 1] = B0;
for (int i = N - 2; i >= 0; i--)
{
double A = H[i] / 6;
double B = (H[i] + H[i + 1]) / 3;
double D = H[i + 1] / 6;
alphas[i] = -A / (B + alphas[i + 1] * D);
}
for (int i = N - 2; i >= 0; i--)
{
double D = H[i + 1] / 6;
double B = (H[i] + H[i + 1]) / 3;
double F = (Y[i + 2] - Y[i + 1]) / H[i + 1] - (Y[i + 1] - Y[i]) / H[i];
betas[i] = (F - betas[i + 1] * D) / (B + D * alphas[i + 1]);
}
for (int i = 1; i <= N - 1; i++)
{
C[i] = alphas[i - 1] * C[i - 1] + betas[i - 1];
}
}
double spline(double x)
{
for (int i = 0; i <= N - 1; i++)
{
if (x <= X[i + 1])
{
return Y[i] * (X[i + 1] - x) / H[i] + Y[i + 1] * (x - X[i]) / H[i] + C[i] * (pow((X[i + 1] - x), 3) - H[i] * H[i] * (X[i + 1] - x)) / (6 * H[i]) + C[i + 1] * (pow((x - X[i]), 3) - H[i] * H[i] * (x - X[i])) / (6 * H[i]);
}
}
}
int main()
{
cout << setprecision(12);
make_spline();
double d2f_left = (spline(a) - 2 * spline(a + eps) + spline(a + 2 * eps)) / (eps*eps);
double d2f_right = (spline(b - 2 * eps) - 2 * spline(b - eps) + spline(b)) / (eps*eps);
cout << "expected A: " << A0 << endl;
cout << "evaluated A: " << d2f_left << endl;
cout << "expected B: " << B0 << endl;
cout << "evaluated B: " << d2f_right << endl;
double r, y_spl;
double eps = pow(10, -3);
FILE *f1 = fopen("C:\\Octave\\SPL\\Spline_S.txt", "w");
for (r = X[0]; r <= X[N]; r = r + eps)
{
y_spl = spline(r);
fprintf(f1, "%.3Lf %.7Lf \n", r, y_spl);
}
fclose(f1);
FILE *f2 = fopen("C:\\Octave\\SPL\\Spline_Dots.txt", "w");
for (int i = 0; i <= N; i++)
{
fprintf(f2, "%.3Lf %.7Lf \n", X[i], Y[i]);
}
fclose(f2);
FILE *f3 = fopen("C:\\Octave\\SPL\\Spline_F.txt", "w");
for (r = X[0]; r <= X[N]; r = r + eps)
{
y_spl = spline(r);
fprintf(f1, "%.3Lf %.7Lf \n", r, func(r));
}
fclose(f3);
//запись в универе
/*
FILE *f1 = fopen("H:\\spl\\Spline_S.txt", "w");
for (r = X[0]; r <= X[N - 1]; r = r + eps)
{
y_spl = spline(X, Y, h, G, r);
fprintf(f1, "%.3Lf %.7Lf \n", r, y_spl);
}
fclose(f1);
FILE *f2 = fopen("H:\\spl\\Spline_Dots.txt", "w");
for (int i = 0; i < N; i++)
{
fprintf(f2, "%.3Lf %.7Lf \n", X[i], Y[i]);
}
fclose(f2);
FILE *f3 = fopen("H:\\spl\\Spline_F.txt", "w");
for (r = X[0]; r <= X[N - 1]; r = r + eps)
{
y_spl = spline(X, Y, h, G, r);
fprintf(f1, "%.3Lf %.7Lf \n", r, func(r));
}
fclose(f1);
*/
system("pause");
return 0;
}
Be the first to comment
You can use [html][/html], [css][/css], [php][/php] and more to embed the code. Urls are automatically hyperlinked. Line breaks and paragraphs are automatically generated.