HJ Blog > Code Snippet #316: Pang 2.2

2009-06-28 01:58:28 Pang 2.2

Going thru some of my computational physics books, in this case, Tao Pang's intro on the subject. Here's my implementation of problem 2.2 on Lagrangian interpolation done up in F-Script. The formula being programmed is:

$Click on image to see the raw LaTeX code$ f(x) = \sum_{k=1}^{n+1} f_k P_k^{(n)}(x) where $Click on image to see the raw LaTeX code$ P_k^{(n)}(x) = \prod_{j \ne k}^{n+1} \frac{x-x_j}{x_k-x_j}

for a fixed set of known values x for f. The code takes advantage of F-Script's (x,y) point type, which has the setting syntax "x<>y", in which one can later read off, say, the y value, via "myPoint y".

( fscript ) ✂

#!/usr/local/bin/fscript

xValues := { 0.3, 0.5 }.
trueFxValues := { 0.328627, 0.520500 }.

knownDataPoints := {
    0.0 <> 0.0,       "f1"
    0.4 <> 0.428392,
    0.8 <> 0.742101,
    1.2 <> 0.910314,
    1.6 <> 0.970348   "f5 = f(n+1)"
}.

weightFunctionP := [ :k :x |
    weight := 1.0.
    xk := (knownDataPoints at:k) x.
    0 to:(knownDataPoints count -1) do:[ :j |
        (j=k) ifFalse:[
            xj := (knownDataPoints at:j) x.
            weight := weight * (x — xj) / (xk — xj).
        ].
    ].
    weight.
].

f := [ :x |
    value := 0.0.
    0 to:(knownDataPoints count -1) do:[ :k |
        fk := (knownDataPoints at:k) y.
        value := value + (fk * (weightFunctionP value:k value:x)).
    ].
    value.
].

[ :index |
    fx := f value:(xValues at:index).
    trueFx := trueFxValues at:index.
    accuracy := ((10000 * (trueFx — fx) / trueFx) intValue abs) / 100. "to 2 places"

    xString := (xValues at:index) stringValue.
    sys out println:
        '\nf('++xString++') = '++(fx stringValue)++'\n'++
        'fTrue('++xString++') = '++(trueFx stringValue)++'\n'++
        '  accuracy: '++(accuracy stringValue)++'%'.

] value:@(xValues count iota).

  1  #!/usr/local/bin/fscript
  2  
  3  xValues := { 0.3, 0.5 }.
  4  trueFxValues := { 0.328627, 0.520500 }.
  5  
  6  knownDataPoints := {
  7      0.0 <> 0.0,       "f1"
  8      0.4 <> 0.428392,
  9      0.8 <> 0.742101,
 10      1.2 <> 0.910314,
 11      1.6 <> 0.970348   "f5 = f(n+1)"
 12  }.
 13  
 14  weightFunctionP := [ :k :x |
 15      weight := 1.0.
 16      xk := (knownDataPoints at:k) x.
 17      0 to:(knownDataPoints count -1) do:[ :j |
 18          (j=k) ifFalse:[
 19              xj := (knownDataPoints at:j) x.
 20              weight := weight * (x — xj) / (xk — xj).
 21          ].
 22      ].
 23      weight.
 24  ].
 25  
 26  f := [ :x |
 27      value := 0.0.
 28      0 to:(knownDataPoints count -1) do:[ :k |
 29          fk := (knownDataPoints at:k) y.
 30          value := value + (fk * (weightFunctionP value:k value:x)).
 31      ].
 32      value.
 33  ].
 34  
 35  [ :index |
 36      fx := f value:(xValues at:index).
 37      trueFx := trueFxValues at:index.
 38      accuracy := ((10000 * (trueFx — fx) / trueFx) intValue abs) / 100. "to 2 places"
 39  
 40      xString := (xValues at:index) stringValue.
 41      sys out println:
 42          '\nf('++xString++') = '++(fx stringValue)++'\n'++
 43          'fTrue('++xString++') = '++(trueFx stringValue)++'\n'++
 44          '  accuracy: '++(accuracy stringValue)++'%'.
 45  
 46  ] value:@(xValues count iota).

which gives the following output:
( terminal ) ✂

f(0.3) = 0.3293449013671874
fTrue(0.3) = 0.328627
  accuracy: 0.21%

f(0.5) = 0.5199387451171874
fTrue(0.5) = 0.5205
  accuracy: 0.1%

  1  f(0.3) = 0.3293449013671874
  2  fTrue(0.3) = 0.328627
  3    accuracy: 0.21%
  4  
  5  f(0.5) = 0.5199387451171874
  6  fTrue(0.5) = 0.5205
  7    accuracy: 0.1%

RSS Feed for this page

Jeff (Sat, January 7th, 2012, 3:27pm UTC)

Interesting follow-up on this. Now that C has blocks, it looked straight-forward to transliterate this F-Scirpt code into C code. Not only was easy but it's actually better-looking code, despite my love of the former. (Also, it would be trivial in the C version to produce 64-bit results.) The numerical results are slightly different, so I include them at the end.

There is just one gotcha on referencing a C array in a block: you can't, you have to reference a pointer to it. So, one extra line of code to clear that up.

Apart from the few lines of set up at the top, it's almost the same length of code, and you'll note several places where the C code is cleaner, thanks to bracket subscripting, += and *= syntax. The complication, of course, was having to set all variable types and creating a dataPair struct.

See what you think:
( c ) ✂

#include <stdio.h>
#include <stdlib.h>
#include <math.h>

#define BOUNDS(a) ((sizeof (a))/(sizeof ((a)[0])))

int main(int argc, const char *argv[]) {
    
    double xValues[] = { .3, .5 };
    double trueFxValues[] = { .328627, .520500 };
    struct dataPair { double x,y; };
    struct dataPair knownDataPointsRaw[] = {
        { 0, 0 },
        {  .4, .428392 },
        {  .8, .742101 },
        { 1.2, .910314 },
        { 1.6, .970348 }
    };
    struct dataPair *knownDataPoints = knownDataPointsRaw;
    int knownDataPairCount = BOUNDS( knownDataPointsRaw );
    
    double (^weightFunctionP)( int k, double x ) = ^( int k, double x ) {
        double weight = 1.0;
        double xk = knownDataPoints[k].x;
        for( int j=0; j<knownDataPairCount; j++)
            if (j!=k) {
                double xj = knownDataPoints[j].x;
                weight *= (x — xj) / (xk — xj);
            }
        return weight;
    };
    
    double (^f)( double x ) = ^( double x ) {
        double value = 0.0;
        for( int k=0; k<knownDataPairCount; k++ ) {
            double fk = knownDataPoints[k].y;
            value += fk * weightFunctionP( k, x );
        }
        return value;
    };
    
    for (int i=0; i<BOUNDS( xValues); i++ ) {
        double x = xValues[ i ];
        double fx = f( xValues[ i ] );
        double trueFx = trueFxValues[ i ];
        double accuracy = fabs(10000 * (trueFx — fx) / trueFx) / 100;
        printf(
            "\nf(%3.1f) = %18.16f\n" \
            "fTrue(%3.1f) = %f\n" \
            "  accuracy: %5.2f%%\n\n",
            x, fx, x, trueFx, accuracy
        );
    }
    
    return EXIT_SUCCESS;
}

  1  #include <stdio.h>
  2  #include <stdlib.h>
  3  #include <math.h>
  4  
  5  #define BOUNDS(a) ((sizeof (a))/(sizeof ((a)[0])))
  6  
  7  int main(int argc, const char *argv[]) {
  8      
  9      double xValues[] = { .3, .5 };
 10      double trueFxValues[] = { .328627, .520500 };
 11      struct dataPair { double x,y; };
 12      struct dataPair knownDataPointsRaw[] = {
 13          { 0, 0 },
 14          {  .4, .428392 },
 15          {  .8, .742101 },
 16          { 1.2, .910314 },
 17          { 1.6, .970348 }
 18      };
 19      struct dataPair *knownDataPoints = knownDataPointsRaw;
 20      int knownDataPairCount = BOUNDS( knownDataPointsRaw );
 21      
 22      double (^weightFunctionP)( int k, double x ) = ^( int k, double x ) {
 23          double weight = 1.0;
 24          double xk = knownDataPoints[k].x;
 25          for( int j=0; j<knownDataPairCount; j++)
 26              if (j!=k) {
 27                  double xj = knownDataPoints[j].x;
 28                  weight *= (x — xj) / (xk — xj);
 29              }
 30          return weight;
 31      };
 32      
 33      double (^f)( double x ) = ^( double x ) {
 34          double value = 0.0;
 35          for( int k=0; k<knownDataPairCount; k++ ) {
 36              double fk = knownDataPoints[k].y;
 37              value += fk * weightFunctionP( k, x );
 38          }
 39          return value;
 40      };
 41      
 42      for (int i=0; i<BOUNDS( xValues); i++ ) {
 43          double x = xValues[ i ];
 44          double fx = f( xValues[ i ] );
 45          double trueFx = trueFxValues[ i ];
 46          double accuracy = fabs(10000 * (trueFx — fx) / trueFx) / 100;
 47          printf(
 48              "\nf(%3.1f) [color= %18.16f\n" \
 49              "fTrue(%3.1f) [color= %f\n" \
 50              "  accuracy: %5.2f%%\n\n",
 51              x, fx, x, trueFx, accuracy
 52          );
 53      }
 54      
 55      return EXIT_SUCCESS;
 56  }

( terminal ) ✂

f(0.3) = 0.3293449013671874
fTrue(0.3) = 0.328627
  accuracy:  0.22%


f(0.5) = 0.5199387451171875
fTrue(0.5) = 0.520500
  accuracy:  0.11%

  1  f(0.3) = 0.3293449013671874
  2  fTrue(0.3) = 0.328627
  3    accuracy:  0.22%
  4  
  5  
  6  f(0.5) = 0.5199387451171875
  7  fTrue(0.5) = 0.520500
  8    accuracy:  0.11%

Jeff (Sat, January 7th, 2012, 4:25pm UTC)

Here's a side-by-side image of the two bits of code for quick comparison: