Script to solve second order (polynomial) interpolation

Shell Programming and Scripting
1

Currently I have awk command to do linear interpolation

awk '
{
     P[$1]=$2
     I[i++]=$1
}
END {
     j=0; s=I[j]; t=I[j+1]
     for(i=m;i<=n;i++) {
          if(I[j+2] && i>t) {
               j++; s=I[j]; t=I[j+1]
          }
          print i, P+(i-s)*(P[t]-P)/(t-s)
     }
} ' m=1 n=8 infile

FILE CONTENT INPUT# a.txt

#X Y
1 22.3125
4 22.5
8 22.1875

Any idea how to change the code to polynomial interpolation?
Assume that Y of X0 = 0. Please help me...

2

What do you mean by polynomial interpolation?

3

Polynomial interpolation is an interpolation that based on three value points (two previous points and next point). So, if we see the data, there are an ID 1, 4, and 8. So, we need to find out the value of 2, 3, 5, 6, and 7 based on value 1, 4, and 8. The formula is that
(((x-x2) * (x-x3)) / ((x1-x2) * (x1-x3))) * y1 + (((x-x1) * (x-x3)) / ((x2-x1) * (x2-x3))) * y2 + (((x-x1) * (x-x2)) / ((x3-x1) * (x3-x1))) * y3
x = current ID;
x1 = the first known ID (second previous known ID); ---> 1
x2 = the second known ID (first previous known ID); ---> 4
x3 = the third known ID (next known ID); ---> 8
y1 = the first known value (the value of ID x1)
y2 = the second known value (the value of ID x2)
y3 = the third known value (the value of ID x3)

4

How about this:

awk '
{ P[$1]=$2 ; m=$1>m?$1:m; x[NR]=$1 ; y[NR]=$2 }
END {
    for(i=1;i<=m;i++) {
       if (i in P) printf "%0.4f %0.4f\n", i, P;
       else printf "%0.4f %0.4f\n", i, \
           (((i-x[2]) * (i-x[3])) / ((x[1]-x[2]) * (x[1]-x[3]))) * y[1] + \
           (((i-x[1]) * (i-x[3])) / ((x[2]-x[1]) * (x[2]-x[3]))) * y[2] + \
           (((i-x[1]) * (i-x[2])) / ((x[3]-x[1]) * (x[3]-x[1]))) * y[3] ;
    }
}' infile
1 Like
5

It is working but shows the wrong result, thus it remains dissolved. Any how, thanks.

6

Damn,

Here is a debug version that show the formula with expanded values, it might help us find the issue.

Note data file should not contain the heading ( #X Y ) line:

awk '
{ P[$1]=$2 ; m=$1>m?$1:m; x[NR]=$1 ; y[NR]=$2 }
END {
    for(i=1;i<=m;i++) {
       if (i in P) printf "%0.4f %0.4f\n", i, P;
       else { printf "%0.4f %0.4f\n", i, \
           (((i-x[2]) * (i-x[3])) / ((x[1]-x[2]) * (x[1]-x[3]))) * y[1] + \
           (((i-x[1]) * (i-x[3])) / ((x[2]-x[1]) * (x[2]-x[3]))) * y[2] + \
           (((i-x[1]) * (i-x[2])) / ((x[3]-x[1]) * (x[3]-x[1]))) * y[3] ;
           printf "   (((%f-%f) * (%f-%f)) / ((%f-%f) * (%f-%f))) * %f +\n",
             i,x[2],i,x[3],x[1],x[2],x[1],x[3],y[1];
           printf "   (((%f-%f) * (%f-%f)) / ((%f-%f) * (%f-%f))) * %f +\n",
             i,x[1],i,x[3],x[2],x[1],x[2],x[3],y[2];
           printf "   (((%f-%f) * (%f-%f)) / ((%f-%f) * (%f-%f))) * %f\n",
             i,x[1],i,x[2],x[3],x[1],x[3],x[1],y[3];
           printf "  = %f + %f + %f\n",
             (((i-x[2]) * (i-x[3])) / ((x[1]-x[2]) * (x[1]-x[3]))) * y[1],
             (((i-x[1]) * (i-x[3])) / ((x[2]-x[1]) * (x[2]-x[3]))) * y[2],
             (((i-x[1]) * (i-x[2])) / ((x[3]-x[1]) * (x[3]-x[1]))) * y[3] ;
        }
    }
}' infile

Exampe output for X=7.0:

7.0000 16.2130
   (((7.000000-4.000000) * (7.000000-8.000000)) / ((1.000000-4.000000) * (1.000000-8.000000))) * 22.312500 +
   (((7.000000-1.000000) * (7.000000-8.000000)) / ((4.000000-1.000000) * (4.000000-8.000000))) * 22.500000 +
   (((7.000000-1.000000) * (7.000000-4.000000)) / ((8.000000-1.000000) * (8.000000-1.000000))) * 22.187500
  = -3.187500 + 11.250000 + 8.150510
1 Like
7

GREAT. It is working. But how to declare if x=0, then y=0. Because for y[3] and y[4], the x[1]=0, thus y[0]=0.

awk '
{ P[$1]=$2 ; m=$1>m?$1:m; x[NR]=$1 ; y[NR]=$2 }
END {
    for(i=1;i<=m;i++) {
       if (i in P) printf "%0.4f %0.4f\n", i, P;
       else { printf "%0.4f %0.4f\n", i, \
           (((i-x[2]) * (i-x[3])) / ((x[1]-x[2]) * (x[1]-x[3]))) * y[1] + \
           (((i-x[1]) * (i-x[3])) / ((x[2]-x[1]) * (x[2]-x[3]))) * y[2] + \
           (((i-x[1]) * (i-x[2])) / ((x[3]-x[1]) * (x[3]-x[2]))) * y[3] ;
           printf "   (((%f-%f) * (%f-%f)) / ((%f-%f) * (%f-%f))) * %f +\n",
             i,x[2],i,x[3],x[1],x[2],x[1],x[3],y[1];
           printf "   (((%f-%f) * (%f-%f)) / ((%f-%f) * (%f-%f))) * %f +\n",
             i,x[1],i,x[3],x[2],x[1],x[2],x[3],y[2];
           printf "   (((%f-%f) * (%f-%f)) / ((%f-%f) * (%f-%f))) * %f\n",
             i,x[1],i,x[2],x[3],x[1],x[3],x[2],y[3];
           printf "  = %f + %f + %f\n",
             (((i-x[2]) * (i-x[3])) / ((x[1]-x[2]) * (x[1]-x[3]))) * y[1],
             (((i-x[1]) * (i-x[3])) / ((x[2]-x[1]) * (x[2]-x[3]))) * y[2],
             (((i-x[1]) * (i-x[2])) / ((x[3]-x[1]) * (x[3]-x[2]))) * y[3] ;
        }
    }
}' infile
8

OK now I see how is should work, try this:

awk '
{ P[$1]=$2 ; m=m?m:$1; n=$1>n?$1:n; x[NR]=$1 ; y[NR]=$2 }
END {
    x[0]=y[0]=0
    p=-2
    for(i=m;i<=n;i++) {
       if (i in P) { p++; printf "%0.4f %0.4f\n", i, P; }
       else {
           printf "%0.4f %0.4f\n", i, \
             (((i-x[p+2]) * (i-x[p+3])) / ((x[p+1]-x[p+2]) * (x[p+1]-x[p+3]))) * y[p+1] + \
             (((i-x[p+1]) * (i-x[p+3])) / ((x[p+2]-x[p+1]) * (x[p+2]-x[p+3]))) * y[p+2] + \
             (((i-x[p+1]) * (i-x[p+2])) / ((x[p+3]-x[p+1]) * (x[p+3]-x[p+2]))) * y[p+3] ;
           printf "   (((%f-%f) * (%f-%f)) / ((%f-%f) * (%f-%f))) * %f +\n",
             i,x[p+2],i,x[p+3],x[p+1],x[p+2],x[p+1],x[p+3],y[p+1];
           printf "   (((%f-%f) * (%f-%f)) / ((%f-%f) * (%f-%f))) * %f +\n",
             i,x[p+1],i,x[p+3],x[p+2],x[p+1],x[p+2],x[p+3],y[p+2];
           printf "   (((%f-%f) * (%f-%f)) / ((%f-%f) * (%f-%f))) * %f\n",
             i,x[p+1],i,x[p+2],x[p+3],x[p+1],x[p+3],x[p+2],y[p+3];
           printf "  = %f + %f + %f\n",
             (((i-x[p+2]) * (i-x[p+3])) / ((x[p+1]-x[p+2]) * (x[p+1]-x[p+3]))) * y[p+1],
             (((i-x[p+1]) * (i-x[p+3])) / ((x[p+2]-x[p+1]) * (x[p+2]-x[p+3]))) * y[p+2],
             (((i-x[p+1]) * (i-x[p+2])) / ((x[p+3]-x[p+1]) * (x[p+3]-x[p+2]))) * y[p+3] ;
        }
    }
}' infile
1 Like
9

Thank you, man. You are GENIUS...