cubic.c

   1 #define _GNU_SOURCE
   2
   3 #include <stdio.h>
   4 #include <math.h>
   5 #include <stdlib.h>
   6 #include <unistd.h>
   7 #include <sys/types.h>
   8 #include <sys/stat.h>
   9 #include <fcntl.h>
  10
  11 int main(int argc, char **argv)
  12 {
  13    double *fp, *u, *y, *fp2;
  14    int fd;
  15    int N, N2, i, k;
  16    double STEP, STEP2;
  17    double xi, x, fx;
  18    double bi, Ai, Bi, Ci, Di;
  19
  20    fd = open("plot.dat", O_WRONLY|O_CREAT|O_TRUNC);
  21    N = atoi(argv[1]);
  22    STEP = (20. - 5.)/N;
  23    N2 = 3;
  24    STEP2 = STEP / N2;
  25
  26    fp = (double *)malloc((N+1)*sizeof(double));
  27    fp2 = (double *)malloc((N+1)*sizeof(double));
  28    u = (double *)malloc((N+1)*sizeof(double));
  29    y = (double *)malloc((N+1)*sizeof(double));
  30
  31    for(i = 0; i < N+1; i++)
  32      {
  33        xi = 5. + i*STEP;
  34        fp[i] = ( sin(xi) -  xi * cos(xi) )/(xi * xi * xi * xi);
  35        dprintf(fd, "%f   %f \n", xi, fp[i]);
  36      }
  37
  38    dprintf(fd, "\n");
  39
  40    y[0] = 0;
  41    u[0] = 0;
  42    for(i = 1; i < N+1; i++)
  43      {
  44        y[i] = (-0.5) / ( 0.5 * y[i-1] + 2. );
  45        bi = ( 6 / ( 2*STEP ) ) * ( (fp[i+1] - fp[i]) / STEP  - (fp[i] - fp[i-1])/ STEP );
  46        u[i] =  ( bi - 0.5 * u[i-1] ) / ( 0.5 * y[i-1] + 2. );
  47      }
  48
  49    fp2[N+1] = 0;
  50    for(i = N; i > 0 ; i--)
  51      {
  52        fp2[i] = u[i] + y[i]*fp2[i+1];
  53      }
  54
  55    for(i = 0; i < N; i++)
  56      {
  57        for(k = 1; k < N2; k++)
  58          {
  59            x = 5. + i*STEP + k*STEP2;
  60            Ai = (- x + (5. + (i+1)*STEP))/STEP;
  61            Bi = (x - (5. + i*STEP))/STEP;
  62            Ci = (STEP*STEP/6) * (Ai * Ai * Ai - Ai);
  63            Di = (STEP*STEP/6) * (Bi * Bi * Bi - Bi);
  64            fx = Ai * fp[i] + Bi * fp[i+1] + Ci * fp2[i] + Di * fp2[i+1];
  65
  66            dprintf(fd, "%f  %f \n", x, fx);
  67          }
  68
  69      }
  70    close(fd);
  71
  72    return 0;
  73 }