X-Git-Url: https://hackdaworld.org/gitweb/?p=physik%2Fcomputational_physics.git;a=blobdiff_plain;f=polynom_interpolation.c;fp=polynom_interpolation.c;h=587faa9ccf67cc5d482feac5267c261353193245;hp=0000000000000000000000000000000000000000;hb=8fc948825d4563572a1172aba6c0b998bc1d907c;hpb=e8839193842d29f1b0fa31ee8c8622d257ceb041

diff --git a/polynom_interpolation.c b/polynom_interpolation.c
new file mode 100644
index 0000000..587faa9
--- /dev/null
+++ b/polynom_interpolation.c
@@ -0,0 +1,64 @@
+#define _GNU_SOURCE
+#include <stdio.h>
+#include <stdlib.h>
+#include <math.h>
+#include "g_plot.h"
+#include <string.h>
+
+#define I_START 5
+#define I_END 20
+#define I_LENGTH (I_END - I_START)
+#define N_I 3
+
+int main(int argc,char **argv) {
+ double *buf,*buf_i;
+ int N;
+ double step,step_i;
+ double x,x_i,x_j,x_k,w;
+ int i,j,k,l;
+ char file[64];
+ int fd;
+
+ if(argc!=2) {
+  printf("usage: %s <N>\n",argv[0]);
+  return 1;
+ }
+
+ N=atoi(argv[1]);
+ step=(double)I_LENGTH/(N+1);
+
+ /* savefile init */
+ strcpy(file,"polynom_interpolation.plot");
+ fd=gp_init(file);
+
+ /* calculate fixpoints */
+ buf=(double *)malloc((N+2)*sizeof(double));
+ for(i=0;i<N+2;i++) {
+  x=I_START+i*step;
+  buf[i]=(sin(x)-x*cos(x))/(x*x*x*x);
+  /* manually, gp_api sux! */
+  dprintf(fd,"%f %f\n",x,buf[i]);
+ }
+
+ /* do interpolation */
+ step_i=step/(N_I+1);
+ buf_i=(double *)malloc((N+1)*N_I*sizeof(double));
+ for(i=0;i<(N+1);i++) {
+  for(l=1;l<N_I+1;l++) {
+   buf_i[i*(N_I+1)+l]=0;
+   x_i=I_START+i*step+l*step_i;
+   for(j=0;j<(N+2);j++) {
+    x_j=I_START+j*step;
+    w=1;
+    for(k=0;k<(N+2);k++) {
+     x_k=I_START+k*step;
+     if(j!=k) w=w*((x_i-x_k)/(x_j-x_k));
+    }
+    buf_i[i*(N_I+1)+l]+=w*buf[j];
+   }
+   dprintf(fd,"%f %f\n",x_i,buf_i[i*(N_I+1)+l]);
+  }
+ }
+     
+ return 1;
+}