11 int main(int argc, char **argv)
13 double *fp, *u, *y, *fp2;
18 double bi, Ai, Bi, Ci, Di;
20 fd = open("plot.dat", O_WRONLY|O_CREAT|O_TRUNC);
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));
31 for(i = 0; i < N+1; i++)
34 fp[i] = ( sin(xi) - xi * cos(xi) )/(xi * xi * xi * xi);
35 dprintf(fd, "%f %f \n", xi, fp[i]);
42 for(i = 1; i < N+1; i++)
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. );
50 for(i = N; i > 0 ; i--)
52 fp2[i] = u[i] + y[i]*fp2[i+1];
55 for(i = 0; i < N; i++)
57 for(k = 1; k < N2; k++)
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];
66 dprintf(fd, "%f %f \n", x, fx);