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:
f(x) = \sum_{k=1}^{n+1} f_k P_k^{(n)}(x) where 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%