spl4

#include "stdafx.h" #include <iostream> #include "math.h" #include <fstream> using namespace std; const int N = 20; const double eps_0 = 1e-2;//pow(10, -); //function to interpolate double func(double x) { return x * sin(x); } double dfunc(double x) { return 2 * cos(x) - x * sin(x); } //spline itself double spline(double* X, double* Y, double* h, double* G, double var) { int j = -1; double spline_val = 0; for (int i = 0; i < N - 1; i++) { if ((var >= X[i]) && (var <= X[i + 1])) j = i; } spline_val = Y[j] * (X[j + 1] - var) / h[j + 1]; spline_val += Y[j + 1] * (var - X[j]) / h[j + 1]; spline_val += G[j] * (pow(X[j + 1] - var, 3) - h[j + 1] * h[j + 1] * (X[j + 1] - var)) / (6 * h[j + 1]); spline_val += G[j + 1] * (pow(var - X[j], 3) - h[j + 1] * h[j + 1] * (var - X[j])) / (6 * h[j + 1]); return spline_val; } int main() { double A, B; double X[N], Y[N], G[N], h[N - 1]; double X_0 = -14, X_N = 15; A = dfunc(X_0); B = dfunc(X_N); // even grid double step_0 = (X_N - X_0) / (N - 1); for (int i = 0; i < N; i++) { X[i] = X_0 + i * step_0; Y[i] = func(X[i]); cout << X[i] << " " << Y[i] << endl; } // for (int i = 1; i < N; i++) { h[i] = X[i] - X[i - 1]; } // //TDMS for (N-1)x(N-1) matrix for G[i] int n = N - 1; double alpha[N - 1], beta[N - 1]; alpha[1] = 0; beta[1] = 0; for (int i = 1; i < N - 1; i++) //back { double a, b, c, d; a = h[i] / 6; b = (h[i] + h[i + 1]) / 3; c = h[i + 1] / 6; d = (Y[i + 1] - Y[i]) / h[i + 1] - (Y[i] - Y[i - 1]) / h[i]; alpha[i + 1] = -c / (a*alpha[i] + b); beta[i + 1] = (d - a * beta[i]) / (a*alpha[i] + b); // cout << a << " " << b << " " << c << " " << d << " " << endl; } G[0] = A; G[N - 1] = B; G[N] = beta[N - 1]; for (int i = N - 1; i-- > 0;) //forward { G[i] = alpha[i + 1] * G[i + 1] + beta[i + 1]; } //G[0] = A; костыль!! //допустим работает cout << endl; //spline value //double var = 9.4444; //or any other value from A to B //double spline_val = 0; //cout << "for x = " << var << " S = " << spline(X, Y, h, G, var) << endl; //DFS double d2fBB, d2fAA; d2fAA = (spline(X, Y, h, G, X[0]) - 2 * spline(X, Y, h, G, X[0] + eps_0) + spline(X, Y, h, G, X[0] + 2 * eps_0)) / pow(eps_0, 2); d2fBB = (spline(X, Y, h, G, X[N - 1] - 2 * eps_0) - 2 * spline(X, Y, h, G, X[N - 1] - eps_0) + spline(X, Y, h, G, X[N - 1])) / pow(eps_0, 2); cout << spline(X, Y, h, G, X[0]) - 2 * spline(X, Y, h, G, X[0] + eps_0) + spline(X, Y, h, G, X[0] + 2 * eps_0) << endl; cout << spline(X, Y, h, G, X[N - 1] - 2 * eps_0) - 2 * spline(X, Y, h, G, X[N - 1] - eps_0) + spline(X, Y, h, G, X[N - 1]) << endl; cout << "A: " << d2fAA << " " << A << endl; cout << "B: " << d2fBB << " " << B << endl; //in .txt file double r, y_spl; double eps = pow(10, -1); FILE *f1 = fopen("Spline_Graph.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 %.7Lf \n", r, func(r), y_spl); } fclose(f1); FILE *f2 = fopen("Spline_Dots.txt", "w"); for (int i = 0; i < N; i++) { fprintf(f2, "%.3Lf %.7Lf \n", X[i], Y[i]); } fclose(f2); system("pause"); return 0; }