The UNIX and Linux Forums Script to solve second order (polynomial) interpolation
 Forums Search Forums Register Forum Rules Man Pages Albums FAQ Members Calendar Search Today's Posts Mark Forums Read

 Shell Programming and Scripting Post questions about KSH, CSH, SH, BASH, PERL, PHP, SED, AWK and OTHER shell scripts and shell scripting languages here.

#1
01-21-2013
 Registered User Join Date: Jan 2013 Posts: 9 Thanks: 5 Thanked 0 Times in 0 Posts
Script to solve second order (polynomial) interpolation

Currently I have awk command to do linear interpolation

Code:
```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[s]+(i-s)*(P[t]-P[s])/(t-s)
}
} ' m=1 n=8 infile```

FILE CONTENT INPUT# a.txt

Code:
```#X Y
1 22.3125
4 22.5
8 22.1875```

Any idea how to change the code to polynomial interpolation?

Last edited by Scrutinizer; 01-21-2013 at 03:54 AM.. Reason: icode tags => code tags
#2
01-21-2013
 Mead Rotor Join Date: Aug 2005 Location: Saskatchewan Posts: 16,371 Thanks: 490 Thanked 2,534 Times in 2,417 Posts
What do you mean by polynomial interpolation?
#3
01-21-2013
 Registered User Join Date: Jan 2013 Posts: 9 Thanks: 5 Thanked 0 Times in 0 Posts
Quote:
 Originally Posted by Corona688 What do you mean by polynomial interpolation?

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
01-21-2013
 Registered User Join Date: Oct 2010 Posts: 2,100 Thanks: 64 Thanked 608 Times in 586 Posts

Code:
```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[i];
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```

 The Following User Says Thank You to Chubler_XL For This Useful Post: Tzeronone (01-21-2013)
#5
01-21-2013
 Registered User Join Date: Jan 2013 Posts: 9 Thanks: 5 Thanked 0 Times in 0 Posts
Quote:
 Originally Posted by Chubler_XL How about this: Code: ```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[i]; 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```
It is working but shows the wrong result, thus it remains dissolved. Any how, thanks.
#6
01-21-2013
 Registered User Join Date: Oct 2010 Posts: 2,100 Thanks: 64 Thanked 608 Times in 586 Posts
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:

Code:
```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[i];
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:

Code:
```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```

Last edited by Chubler_XL; 01-21-2013 at 11:24 PM.. Reason: Fix indenting
 The Following User Says Thank You to Chubler_XL For This Useful Post: Tzeronone (01-22-2013)
#7
01-22-2013
 Registered User Join Date: Jan 2013 Posts: 9 Thanks: 5 Thanked 0 Times in 0 Posts
Quote:
 Originally Posted by Chubler_XL 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: Code: ```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[i]; 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: Code: ```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```
GREAT. It is working. But how to declare if x=0, then y[x]=0. Because for y[3] and y[4], the x[1]=0, thus y[0]=0.

Code:
```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[i];
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```

Last edited by Tzeronone; 01-22-2013 at 02:00 AM.. Reason: icode tags changed to code tags

 Tags awk, numbering, text processing

 More UNIX and Linux Forum Topics You Might Find Helpful Thread Thread Starter Forum Replies Last Post nex_asp Shell Programming and Scripting 12 12-05-2012 08:38 AM ardy_yana Shell Programming and Scripting 7 04-05-2012 03:04 AM z1dane Web Programming 0 01-27-2010 12:10 AM egkumpe UNIX for Dummies Questions & Answers 0 08-05-2004 06:05 PM

All times are GMT -4. The time now is 04:57 AM.