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 Sep 28 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$ – Sandwiches29 Sep 28 at 2:21
  • 1
    $\begingroup$ Is this what you need? $\endgroup$ – Αλέξανδρος Ζεγγ Sep 28 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$ – 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]] - 
  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]
| improve this answer | |
  • $\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$ – Sandwiches29 Sep 30 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$ – Daniel Huber Sep 30 at 9:05

Your Answer

By clicking “Post Your Answer”, you agree to our terms of service, privacy policy and cookie policy

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