Hello,


I needed to make some interpolations in my project, which led me to this
macro that I'm sharing here.


// ---------------------------------------------------------------------
// Macro LagrangeInterpolation()
//
// Input :
// - TabularValues : a two dimensional array than contain the x and the
//                   y coordinates.
// - xValue : the given abscissae.
// Output :
// - The interpolated value for the given xValue.
//
//
// Note 1 : A remarkable feature of Lagrange interpolation is that
//          the values entered initially do not have to be in order,
//          or evenly spaced. Accuracy is usually better with uniform
//          spacing, however.
//
// Note 2 : It's remarkable that, even for values of x outside the given
//          range, the Lagrange interpolation formula still produces
//          fairly good values.
// ---------------------------------------------------------------------
#macro LagrangeInterpolation(TabularValues, xValue)
 #local n = dimension_size(TabularValues,1);
 #local yValue = 0.0;
 #local i = 1;
 #while ( i < n )
  #local p = 1;
  #local j = 1;
  #while ( j < n )	
   #if(i != j )
     #local p = p*(xValue - TabularValues[j][0])/(TabularValues[i][0] -
TabularValues[j][0]);
   #end
   #local j = j + 1;
  #end	
  #local yValue = yValue + p*TabularValues[i][1];
  #local i = i +1;
 #end
 yValue
#end


With a small example.

// Six values
#declare N = 6;
// First table column : angle in degrees
// Second column : nothing for the moment
#declare SinusData = array[N][2] {
 {29.43,0},
 {30.97,0},
 {27.69,0},
 {28.11,0},
 {31.58,0},
 {33.05,0}
 }
// Filling the second column with sinuses	
#declare i = 0;
#while( i < N)
 #declare SinusData[i][1]=sin(radians(SinusData[i][0]));
 // Display array
 #debug concat(str(i,2,0)," : ",str(SinusData[i][0],0,2),
 " | ",str(SinusData[i][1],0,6),"\n")
 #declare i = i + 1;
#end

// A value within the range or not
#declare xp = 30.0;
// Interpolation for this value
#declare Val = LagrangeInterpolation(SinusData, xp);
// Because real value is know, calculate the difference
#declare e = sin(radians(xp))-Val;
// And displaying the results
#debug concat("For x = ",str(xp,0,2)," interpolated value is
",str(Val,0,6)," with a gap of ",str(e,0,6),".\n")

Which gives :

 0 : 29.43 | 0.491360
 1 : 30.97 | 0.514589
 2 : 27.69 | 0.464688
 3 : 28.11 | 0.471166
 4 : 31.58 | 0.523689
 5 : 33.05 | 0.545371

For x = 30.00 interpolated value is 0.500000 with a gap of -0.000000.




If that helps...
-- 
Kurtz le pirate
Compagnie de la Banquise

Post a reply to this message