123456789101112131415161718192021222324252627282930313233343536373839404142434445464748495051525354555657585960616263646566676869707172737475767778798081828384858687888990919293949596979899100101102103104105106107108109110111112113114115116117118119120121122123124125126127128129130131132133134135136137138139140141142143144145146147148149150151152153154 |
- #include <stdio.h>
- #include <stdlib.h>
- #include <math.h>
- typedef struct{
- double a, b, c, d;
- double x;
- }CubicSpline;
- double *EquallySpaced(double a, double b, int node){
- int i;
- double step;
- step = (b-a)/(node-1);
- double* x = (double*)malloc(node*sizeof(double));
- for (i=0; i<node; i++){
- x[i] = a + i*step;
- }
- return x;
- }
- void FindCubicSplines(CubicSpline* splines, int n, double *xk, double *yk, double dfa, double dfb)
- {
- double* h = (double*)malloc((n-1)*sizeof(double));
-
- double* b = (double*)malloc(n * sizeof(double));
- double* c = (double*)malloc(n * sizeof(double));
- double* d = (double*)malloc(n * sizeof(double));
- double* f = (double*)malloc(n * sizeof(double));
- double* M = (double*)malloc(n * sizeof(double));
- int i;
- b[0] = 0;
- d[n-1] = 0;
- for (i = 0; i < n - 1; i++)
- {
- h[i] = xk[i+1] - xk[i];
- }
- c[0] = -h[0] / 3;
- d[0] = -h[0] / 6;
- b[n-1] = h[n-2] / 6;
- c[n-1] = h[n-2] / 3;
- f[0] = dfa - (yk[1]-yk[0]) / h[0];
- f[n-1] = dfb - (yk[n-1]-yk[n-2]) / h[n-2];
- for(i = 1; i < n - 1; i++)
- {
- b[i] = h[i-1] / (h[i-1] + h[i]);
- c[i] = 2;
- d[i] = h[i] / (h[i-1] + h[i]);
- f[i] = (6 / (h[i-1] + h[i])) * ((yk[i+1] - yk[i]) / h[i] - (yk[i] - yk[i-1]) / h[i-1]);
- }
- //Найти Mi
- for(i = 1; i < n; i++)
- {
- c[i] = c[i] - b[i] * d[i-1] / c[i-1];
- f[i] = f[i] - b[i] * f[i-1]/c[i-1];
- b[i] = 0;
- }
- M[n-1] = f[n-1] / c[n-1];
- for(i = n - 2; i >= 0; i--)
- {
- M[i] = (f[i] - d[i] * M[i+1]) / c[i];
- }
-
- //Найти коэффициенты
- for(i = 0; i < n - 1; i++){
- splines[i].a = M[i] / (6*h[i]);
- splines[i].b = M[i+1] / (6*h[i]);
- splines[i].c = (yk[i+1] - yk[i]) / h[i] - h[i] * (M[i+1]-M[i]) / 6;
- splines[i].d = yk[i] - M[i] * h[i] * h[i] / 6;
- splines[i].x = xk[i];
- }
- splines[n-1].x = xk[n-1];
- free(h);
- free(b);
- free(c);
- free(d);
- free(f);
- free(M);
- }
- double cubicSplineInterpolate(CubicSpline* splines, int n, double x) {
- int i;
- double dx1, dx2;
- for (i = 0; i < n - 1; i++) {
- if (x >= splines[i].x && x <= splines[i+1].x) {
- dx1 = x - splines[i].x;
- dx2 = splines[i+1].x - x;
- return splines[i].a*dx2*dx2*dx2 + splines[i].b*dx1*dx1*dx1 + splines[i].c*dx1 + splines[i].d;
- }
- }
- return 0;
- }
- double Func1(double x) {
- return 1.0*(exp(x) - exp(-x))/(exp(x)+exp(-x));
- }
- double Func2(double x) {
- double sign_x;
- if (x>0) sign_x =1;
- else if (x==0) sign_x =0;
- else sign_x = -1;
- return 3 * sign_x * pow(x*1.0, 4) + 4 * pow(x*1.0, 3) - 12 * pow(x*1.0, 2) + 4;
- }
- int main() {
- FILE *f = fopen("data.txt", "w");
- if (f == NULL) {
- printf("Error opening file!\n");
- return 1;
- }
- int nodes[] = {5,7,10};
- int i,j,n;
- double dfa, dfb, a, b, x, y;
- dfa = 48;
- dfb = 31.5;
- a = -1; b = 1.5;
- for(i=0; i<3;i++){
- n = nodes[i];
- double *xk;
- CubicSpline* splines = (CubicSpline*)malloc(n * sizeof(CubicSpline));
- double* yk = (double*)malloc(n*sizeof(double));
-
- fprintf(f,"%d ", n);
- xk = EquallySpaced(a, b, n);
- for(j =0; j< n; j++) yk[j] = Func2(xk[j]);
- FindCubicSplines(splines, n, xk, yk, dfa, dfb);
-
- for(j=0; j<1000;j++){
- x = a+j*(b-a)/999;
- y = cubicSplineInterpolate(splines, n, x);
- fprintf(f, "%.13f ", y);
- }
- fprintf(f,"\n");
- free(splines);
- free(xk);
- free(yk);
- }
- fclose(f);
- return 0;
- }
|