# Newton's divided differences iteration

I am trying to write a program that forms the interpolation polynomial for a given function on a given interval for any number of data points n. I wish to write a formula that will compute all of the necessary divided differences. Here is part of my code,

g[x_] := 1/(1 + x^2)
f = Table[g[x], {x, -5.0, 5.0, 1}]
x = Table[a, {a, -5, 5, 1}]

Table[(f[[i]] - f[[i - 1]])/(x[[i + j]] - x[[i - 1]]), {j, 0, 9}, {i, 2, 11 - j}]


For j=0 I get a list of the correct outputs I need, {0.020362, 0.0411765, 0.1, 0.3, 0.5, -0.5, -0.3, -0.1, -0.0411765,-0.020362} but for j=1, I must find a way to have the terms in the numerator replaced with these output values over the loop, and so forth for each j. So for j=1,i=2 I wish to write a formula that will calculate (0.0411765-0.020362)/(x[] - x[]) for example, if this is possible.

I am a Mathematica novice, so bear with me and I hope I have been as clear as possible. I would appreciate any help, tips, tricks or guidance if it seems my approach is not a good one.

• I don't if this is just a (homework) exercise, but if you want to use Mathematica, either InterpolatingPolynomial[Transpose@{xx, ff}, x] or Interpolation[Transpose@{xx, ff}, InterpolationOrder -> All][x] constructs the interpolation you seek. – Michael E2 Sep 28 at 2:12
• Yeah, this is part of an assignment. I am aware of the built-in interpolation but my professor explicitly mentioned not to use them and to "get our hands dirty a bit." I am able to write something but it is highly impractical and essentially impossible for large amounts of data points. I essentially wrote an entirely new function for each j, calling the results from the previous outputs. – Sandwiches29 Sep 28 at 2:21
• Is this what you need? – Αλέξανδρος Ζεγγ Sep 28 at 6:07
• I've seen this post, and it obviously produces the result that I am looking for. I was just interested in understanding whether my approach was a valid one or not, or whether my work could be adjusted to produce the same results. – Sandwiches29 Sep 29 at 16:37

Here is one possibility (although a bit a tricky one):

We first define a function that returns the Newton polynomial as a function :

getNewton[xs_, ys_] := Module[{dd, i, newtoncoef},
dd[{i_Integer}] := ys[[i]] ;(*dd implements Newton recursion*)
dd[i : {__} ] := (dd[ Rest@i] - dd[Most@i])/(xs[[Last@i]] -
xs[[First@i]]);

newtoncoef =
dd[Range[#]] & /@ Range[Length[xs]];(*Newton coefficients*)
Function[x, newtoncoef.FoldList[#1 (x - #2) &, 1, Most@xs]]
]


Now we try it out. We first create some data and then apply our function to the x values:

g[x_] := 1/(1 + x^2);
ys = Table[g[x], {x, -5.0, 5.0, 1}];
xs = Table[x, {x, -5, 5, 1}];
Print["xs=", xs]; Print["ys=", ys];

Print["Newton poly at data points= ", getNewton[xs, ys] /@ xs]

• Thank you for your help, it is much appreciated. I tried running this code, but rather than printing me the Newton polynomial, it printed a list of ys values. – Sandwiches29 Sep 30 at 0:11
• @Sandwiches29 That is correct. getNewton does everything in one go. The coefficients are calculated by: "newtoncoef". getNewton1 will return the coefficients: getNewton1[xs_, ys_] := Module[{dd, i, newtoncoef}, dd[{i_Integer}] := ys[[i]];(dd implements Newton recursion) dd[i : {__}] := (dd[Rest@i] - dd[Most@i])/(xs[[Last@i]] - xs[[First@i]]); newtoncoef = dd[Range[#]] & /@ Range[Length[xs]] ] – Daniel Huber Sep 30 at 9:05