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[[3]] - x[[1]]) 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.

  • $\begingroup$ 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. $\endgroup$
    – Michael E2
    Commented Sep 28, 2020 at 2:12
  • $\begingroup$ 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. $\endgroup$ Commented Sep 28, 2020 at 2:21
  • 1
    $\begingroup$ Is this what you need? $\endgroup$ Commented Sep 28, 2020 at 6:07
  • $\begingroup$ 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. $\endgroup$ Commented Sep 29, 2020 at 16:37

1 Answer 1


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]] - 
  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]
  • $\begingroup$ 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. $\endgroup$ Commented Sep 30, 2020 at 0:11
  • $\begingroup$ @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]] ] $\endgroup$ Commented Sep 30, 2020 at 9:05

Not the answer you're looking for? Browse other questions tagged or ask your own question.