*** empty log message ***

[physik/computational_physics.git] / cubic.c

diff --git a/cubic.c b/cubic.c
new file mode 100644 (file)
index 0000000..3e108d5
--- /dev/null
+++ b/cubic.c
@@ -0,0 +1,73 @@
+#define _GNU_SOURCE
+
+#include <stdio.h>
+#include <math.h>
+#include <stdlib.h>
+#include <unistd.h>
+#include <sys/types.h>
+#include <sys/stat.h>
+#include <fcntl.h>
+
+int main(int argc, char **argv)
+{
+   double *fp, *u, *y, *fp2;
+   int fd;
+   int N, N2, i, k;
+   double STEP, STEP2;
+   double xi, x, fx; 
+   double bi, Ai, Bi, Ci, Di;   
+   
+   fd = open("plot.dat", O_WRONLY|O_CREAT|O_TRUNC);
+   N = atoi(argv[1]);
+   STEP = (20. - 5.)/N;
+   N2 = 3;
+   STEP2 = STEP / N2;
+
+   fp = (double *)malloc((N+1)*sizeof(double));
+   fp2 = (double *)malloc((N+1)*sizeof(double));
+   u = (double *)malloc((N+1)*sizeof(double));
+   y = (double *)malloc((N+1)*sizeof(double));
+
+   for(i = 0; i < N+1; i++)
+     {
+       xi = 5. + i*STEP;
+       fp[i] = ( sin(xi) -  xi * cos(xi) )/(xi * xi * xi * xi);
+       dprintf(fd, "%f   %f \n", xi, fp[i]);       
+     }
+
+   dprintf(fd, "\n");
+  
+   y[0] = 0;
+   u[0] = 0;
+   for(i = 1; i < N+1; i++)
+     { 
+       y[i] = (-0.5) / ( 0.5 * y[i-1] + 2. );
+       bi = ( 6 / ( 2*STEP ) ) * ( (fp[i+1] - fp[i]) / STEP  - (fp[i] - fp[i-1])/ STEP );
+       u[i] =  ( bi - 0.5 * u[i-1] ) / ( 0.5 * y[i-1] + 2. );
+     }
+
+   fp2[N+1] = 0;
+   for(i = N; i > 0 ; i--)
+     {
+       fp2[i] = u[i] + y[i]*fp2[i+1];
+     }
+
+   for(i = 0; i < N; i++)
+     {
+       for(k = 1; k < N2; k++)
+        {
+          x = 5. + i*STEP + k*STEP2;
+           Ai = (- x + (5. + (i+1)*STEP))/STEP;
+           Bi = (x - (5. + i*STEP))/STEP;
+           Ci = (STEP*STEP/6) * (Ai * Ai * Ai - Ai);
+           Di = (STEP*STEP/6) * (Bi * Bi * Bi - Bi);
+           fx = Ai * fp[i] + Bi * fp[i+1] + Ci * fp2[i] + Di * fp2[i+1];
+
+           dprintf(fd, "%f  %f \n", x, fx); 
+         }
+ 
+     }
+   close(fd);
+
+   return 0;
+}