polynom_interpolation.c

   1 #define _GNU_SOURCE
   2 #include <stdio.h>
   3 #include <stdlib.h>
   4 #include <math.h>
   5 #include "g_plot.h"
   6 #include <string.h>
   7
   8 #define I_START 5
   9 #define I_END 20
  10 #define I_LENGTH (I_END - I_START)
  11 #define N_I 3
  12
  13 int main(int argc,char **argv) {
  14  double *buf,*buf_i;
  15  int N;
  16  double step,step_i;
  17  double x,x_i,x_j,x_k,w;
  18  int i,j,k,l;
  19  char file[64];
  20  int fd;
  21
  22  if(argc!=2) {
  23   printf("usage: %s <N>\n",argv[0]);
  24   return 1;
  25  }
  26
  27  N=atoi(argv[1]);
  28  step=(double)I_LENGTH/(N+1);
  29
  30  /* savefile init */
  31  strcpy(file,"polynom_interpolation.plot");
  32  fd=gp_init(file);
  33
  34  /* calculate fixpoints */
  35  buf=(double *)malloc((N+2)*sizeof(double));
  36  for(i=0;i<N+2;i++) {
  37   x=I_START+i*step;
  38   buf[i]=(sin(x)-x*cos(x))/(x*x*x*x);
  39   /* manually, gp_api sux! */
  40   dprintf(fd,"%f %f\n",x,buf[i]);
  41  }
  42
  43  /* do interpolation */
  44  step_i=step/(N_I+1);
  45  buf_i=(double *)malloc((N+1)*N_I*sizeof(double));
  46  for(i=0;i<(N+1);i++) {
  47   for(l=1;l<N_I+1;l++) {
  48    buf_i[i*(N_I+1)+l]=0;
  49    x_i=I_START+i*step+l*step_i;
  50    for(j=0;j<(N+2);j++) {
  51     x_j=I_START+j*step;
  52     w=1;
  53     for(k=0;k<(N+2);k++) {
  54      x_k=I_START+k*step;
  55      if(j!=k) w=w*((x_i-x_k)/(x_j-x_k));
  56     }
  57     buf_i[i*(N_I+1)+l]+=w*buf[j];
  58    }
  59    dprintf(fd,"%f %f\n",x_i,buf_i[i*(N_I+1)+l]);
  60   }
  61  }
  62
  63  return 1;
  64 }