Take the 2-minute tour ×
Mathematica Stack Exchange is a question and answer site for users of Mathematica. It's 100% free, no registration required.

Given some data pairs $(x_i,y_i)$, with $i=0,...,m$, and a degree $r$, I wish to build a piecewise polynomial function to interpolate these data. That interpolation should be continuous, and, on every interval $[x_k,x_{k+r}]$, with $k=0, r, 2r, ...$, should be a polynomial of degree $r$. This can be useful for example to represent the solution of a PDE obtained with finite element method of degree $r$.

Because with $r=1$ there are no problem I'll refer to $r=2$. For example, given the following data pairs:

$$ \{ (0,0), (1,1), (2,0), (3,1), (4,0) \} $$

I wish to get this result:

The resulting interpolation function need not to have continuous derivative at $x=2$.

I tried with Interpolation and various options:

Interpolation[{{0, 0}, {1, 1}, {2, 0}, {3, 1}, {4, 0}}, 
 InterpolationOrder -> 2, Method -> "Hermite"]
Plot[%[x], {x, 0, 4}]

Interpolation[{{0, 0}, {1, 1}, {2, 0}, {3, 1}, {4, 0}}, 
 InterpolationOrder -> 2, Method -> "Spline"]
Plot[%[x], {x, 0, 4}]

As a reference, under MATLAB, I can build a piecewise polynomial interpolation of arbitrary degree, in a some involved way, with mkpp, and later consume the interpolation with ppval. For piecewise linear interpolation there is a more simple and direct interp1 function.

Under MATLAB I give to mkpp the values of the polinomials and their derivatives at $x_0, x_r, x_{2r}, ...$ and I get the expected result. Under Mathematica this approach doesn't work:

Interpolation[{{{0}, 0, 2, -2}, {{2}, 0, 2, -2}, {{4}, 0, 2, -2}}, 
 InterpolationOrder -> 2, Method -> "Hermite"]
Plot[%[x], {x, 0, 4}]

I considered using Piecewise and constructiong an Interpolation or a polynomial pure function for every interval $[x_k,x_{k+r}]$ but I suspect this become unmanageably complex when there are hundred or thousand of intervals.

There is some builtin way, reasonably simple and fast, to get this result? Naturally I search a general way, for general data and general $r$.

UPDATE

@kguler answer is interesting but I need a way to generalize for every $r$.

share|improve this question

2 Answers 2

Based on What's inside InterpolatingFunction[{{1., 4.}}, <>]?, I would guess that a built-in way is not possible. However, one can take advantage of InterpolatingFunction to construct a Piecewise function. Here, split, does an overlapping partition starting a new list at position p, is a modification of Mr.Wizard's dPcore in this answer.

split[L_, pos_] := 
 Inner[L[[# - 1 ;; #2]] &, Prepend[pos + 1, 2], Append[pos, Length[L]], Head@L]

pwpolyifn[pts_, breaks_] := Function[x,
  Evaluate@
   Piecewise[{Interpolation[#, InterpolationOrder -> All][
        x], #[[1, 1]] <= x <= #[[-1, 1]]} & /@ split[pts, breaks]
    ]]

Example: split can do a ragged split.

split[{{0, 0}, {1, 1}, {2, 0}, {3, 1}, {4, 0}, {5, 1}, {6, 0}}, {3, 6}]
(*
  {{{0, 0}, {1, 1}, {2, 0}},
   {{2, 0}, {3, 1}, {4, 0}, {5, 1}},
   {{5, 1}, {6, 0}}}
*)

In this one, the OP's example, split[.., {3, 5}] is the same as Partition[.., 3, 2]:

split[{{0, 0}, {1, 1}, {2, 0}, {3, 1}, {4, 0}, {5, 1}, {6, 0}}, {3, 5}]
(*
  {{{0, 0}, {1, 1}, {2, 0}},
   {{2, 0}, {3, 1}, {4, 0}},
   {{4, 0}, {5, 1}, {6, 0}}}
*)

pwpolyifn[{{0, 0}, {1, 1}, {2, 0}, {3, 1}, {4, 0}, {5, 1}, {6, 0}}, {3, 5}][x]

Mathematica graphics

Plot[pwpolyifn[{{0, 0}, {1, 1}, {2, 0}, {3, 1}, {4, 0}, {5, 1}, {6, 0}}, {3, 5}][x],
    {x, 0, 6}]

Mathematica graphics

Another example:

SeedRandom[1];
pts = Table[{i, RandomReal[{0, 10}]}, {i, 0, 20}];
breaks = {3, 8, 10, 16};
Plot[Evaluate@pwpolyifn[pts, breaks][x], {x, 0, 20}, 
 GridLines -> {pts[[breaks]][[All, 1]], None}]

Mathematica graphics

If bullet-proofing the definition is desired, then one can check that pts is a list of pairs of numbers and that the breaks are increasing and lie between 2 and Length[pts] - 1.

share|improve this answer
    
Instead of split can't you just do Partition[list, r+1, r]? –  Rahul Narain Jul 28 at 17:25
    
@RahulNarain Maybe I should show it for the last example. It's similar to Internal`PartitionRagged, except with the overlap like your Partition example. –  Michael E2 Jul 28 at 17:26
    
Ah, right, one can't do the last example with Partition. In the original question the degree is constant, though, so Partition should be sufficient. +1 for the Piecewise and InterpolationOrder -> All. –  Rahul Narain Jul 28 at 17:28
    
What exactly All means as a value for the option InterpolationOrder of the Interpolation function? –  unlikely Jul 28 at 18:00
1  
@unlikely It means the degree should be as high as possible, the length of the list minus 1. (Originally I had Length[pts] - 1.) –  Michael E2 Jul 28 at 18:04

Perhaps too specific to OP's example case:

Interpolation[{{{0}, 0, Automatic}, {{1}, 1, 0}, {{2}, 0,  Automatic}, 
               {{3}, 1, 0}, {{4}, 0. Automatic}}, 
  InterpolationOrder -> 2]
  Plot[%[x], {x, 0, 4}]

enter image description here

and

 Interpolation[{{{0}, 0, Automatic}, {{1}, 1, 0}, {{2}, 0,  Automatic}, 
                {{3}, 1, 0}, {{4}, 0. Automatic}}, 
    InterpolationOrder -> 2, PeriodicInterpolation -> True]
 Plot[%[x], {x, 0, 8}]

enter image description here

share|improve this answer
    
You're probably right -- the OP does say he wants the polynomials all the same degree. –  Michael E2 Jul 28 at 17:37
1  
This apperars interesting... What exactly Automatic implies? The documentation is a bit lachonic. And for general $r$ and general case how should I constuct the data to pass to Interpolation? –  unlikely Jul 28 at 17:51
    
@unlikely, the only info in the documentation is in the Details section: "Partial derivatives not specified explicitly can be given as Automatic". I will post an update, if i can figure out a generalization to arbitrary r. –  kguler Jul 28 at 18:17

Your Answer

 
discard

By posting your answer, you agree to the privacy policy and terms of service.

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