6

I am finally working on my n-point Pade code, again, and I am running across an error that was not occurring previously. The heart of the matter revolves around this code:

zi = {0.1, 0.2, 0.3}
ai = {0.904837, 1.05171, -0.499584}

Quiet[ RecurrenceTable[ {A[0] == 0, A[1] == ai[[1]],
          A[n+1]==A[n] + (z - zi[[n]]) ai[[n+1]] A[n-1]},
          A, {n, Length@ai -1 } ],
   {Part::pspec}]

(The use of Quiet is necessary as Part complains about zi[[n]] and ai[[n+1]] when n is purely symbolic.) The code itself is part of a function that I want a symbolic result from, hence z is a Symbol. But, when I run the above code I get the error:

RecurrenceTable::nlnum1: 
  The function value {0.904837,0.904837+0. z} is not a list of numbers with 
  dimensions {2} when the arguments are {0,0.,0.904837}.

Note the term {0.904837,0.904837+0. z} where 0. z is not reduced to zero. What do I need to do to force it to evaluate to zero, as it seems to be the source of the problem? Are there alternatives?

Additionally, as a general complaint about the help system for the Wolfram Research personnel who haunt stackoverflow: in v.7 RecurrenceTable::nlnum1 is not searchable! Nor, does the >> link at the end of the error take you to the error definition, but takes you to the definition of RecurrenceTable, instead, where the common errors are not cross-referenced.

Edit: After reviewing my code, the solution I came up with was to evaluate the RecurrenceTable completely symbolically, including the initial conditions. The working code is as follows:

Clear[NPointPade, NPointPadeFcn]
NPointPade[pts : {{_, _} ..}] := NPointPade @@ Transpose[pts]

NPointPade[zi_List, fi_List] /; Length[zi] == Length[fi] :=
 Module[{ap, fcn, rec},
  ap = {fi[[1]]};
  fcn = Module[{gp = #, zp, res},
     zp =  zi[[-Length@gp ;;]];
     res = (gp[[1]] - #)/((#2 - zp[[1]]) #) &[Rest@gp, Rest@zp];
     AppendTo[ap, res[[1]]];
     res
  ] &;

  NestWhile[fcn, fi, (Length[#] > 1 &)];

 (*
  The recurrence relation is used twice, with different initial conditions, so
  pre-evaluate it to pass along to NPointPadeFcn
 *)
 rec[aif_, zif_, a_, b_][z_] := 
  Evaluate[RecurrenceTable[
     {A[n + 1] == A[n] + (z - zif[n])*aif[n + 1]*A[n - 1], 
      A[0] == a, A[1] == b}, 
     A, {n, {Length@ap - 1}}][[1]]];

  NPointPadeFcn[{zi, ap, rec }]
 ]

NPointPadeFcn[{zi_List, ai_List, rec_}][z_] /; Length[zi] == Length[ai] :=
 Module[{aif, zif},
  zif[n_Integer] /; 1 <= n <= Length[zi] := zi[[n]];
  aif[n_Integer] /; 1 <= n <= Length[zi] := ai[[n]];
  rec[aif, zif, 0, ai[[1]]][z]/rec[aif, zif, 1, 1][z]
 ]

Format[NPointPadeFcn[x_List]] := NPointPadeFcn[Shallow[x, 1]];

Like the built-in interpolation functions, NPointPade does some pre-processing, and returns a function that can be evaluated, NPointPadeFcn. The pre-processing done by NPointPade generates the list of ais from the zis and the function values at those points, in addition to pre-evaluating the recurrence relations. When NPointPadeFcn is supplied with a z value, it evaluates two linear recurrence relations by supplying it with the appropriate values.

Edit: for the curious, here's NPointPade in operation

NPointPade in action

In the first plot, it is difficult to tell the difference between the two functions, but the second plot shows the absolute (blue) and relative (red) errors. As written, it takes a very long time to create a Pade for 20 points, so I need to work on speeding it up. But, for now, it works.

Community
  • 1
  • 1
rcollyer
  • 10,475
  • 4
  • 48
  • 75
  • @rcollyer Using `WorkingPrecision->Infinity` in both calls to `RecurrenceTable` fixes the problem. The issue was that `tmp` contained `aif[1.]` and the likes, which did not match against your patterns following `RecurrenceTable` calls. – Sasha Apr 26 '11 at 18:42
  • @Sasha, I fixed it another way entirely, see above. I'm still giving you credit, as your idea of delaying the evaluation of the `ai`s was the correct one. As a side not, changing the `WorkingPrecision` doesn't seem to work for me. – rcollyer Apr 26 '11 at 18:52
  • For my curiosity, would you give an example using your updated function? – Mr.Wizard Apr 26 '11 at 19:54
  • Thanks. Why does your syntax highlighting show pink/magenta `#` slot symbols? – Mr.Wizard Apr 26 '11 at 21:16
  • @Mr. Wizard: probably because he changed them? You can change colors in preferences. No need to edit internal noteboos. –  Apr 26 '11 at 22:36
  • @d'o-o'b I suppose that could be, but by default a similar color is used to indicate a Slot without a matching Function. – Mr.Wizard Apr 26 '11 at 22:46
  • @Mr. Wizard: Isn't it red in that case and not pink/magenta? –  Apr 26 '11 at 22:48
  • @d'o-o'b not on my system; maybe I am the one using strange syntax highlighting, or maybe it varies from platform to platform? If you type `# + #2` on a new line without anything else, what color are `#` and `#2`? – Mr.Wizard Apr 26 '11 at 22:54
  • I think I have a couple of improvements to your `NPointPade` function: `NPointPade[pts : {{_, _} ..}] := NPointPade @@ (pts\[Transpose])` And an inner function: `res = (gp[[1]] - #)/((#2 - zp[[1]]) #) &[Rest@gp, Rest@zp];` These test faster for me. Also, you seem to have a couple of orphan symbols in the `Module` in `NPointPadeFcn`: `A` and `B` – Mr.Wizard Apr 26 '11 at 23:21
  • @Mr.Wizard, as to the colors: the color I use for a non-matching `Slot` is very red, as the default matching and non-matching colors were very difficult to readily differentiate on my screen. Orphan symbols, argh! I will fix those, and apply your fixes. – rcollyer Apr 26 '11 at 23:57
  • @rcollyer Could you give me a reference to the logic of this code, i.e. the method used to constuct n-Pade approximation ? Why do you call it Pade ? – Sasha Apr 27 '11 at 03:24
  • @Sasha, the standard is [Essentials of Pade Approximants](http://www.amazon.com/gp/search?index=books&linkCode=qs&keywords=012074855X) by Baker, and the method is called n-point Pade approximation. This is more correctly referred to as *rational interpolation*, but since it uses the ideas in forming Pade approximants, either is correct. Through your linke-in profile, I have access to your email, and I can send you a summary of the method tomorrow, if you'd like. The comment field is not the best place for this. – rcollyer Apr 27 '11 at 03:46

2 Answers2

6

You can hide part extraction behind a function:

In[122]:= zi = {0.1, 0.2, 0.3};
ai = {0.904837, 1.05171, -0.499584};

In[124]:= zif[n_Integer] /; 1 <= n <= Length[zi] := zi[[n]]
aif[n_Integer] /; 1 <= n <= Length[ai] := ai[[n]]

In[127]:= RecurrenceTable[{A[0] == 0, A[1] == aif[1], 
  A[n + 1] == 
   A[n] + (z - zif[n]) aif[n + 1] A[n - 1]}, A, {n, (Length@ai) - 1}]

Out[127]= {0.904837, 0.904837, 
 0.904837 - 0.271451 aif[4] + 0.904837 z aif[4]}


EDIT

Here is the work-around for the problem:

In[4]:= zi = {0.1, 0.2, 0.3};
ai = {0.904837, 1.05171, -0.499584};

In[6]:= zif[n_Integer] /; 1 <= n <= Length[zi] := zi[[n]]
aif[n_Integer] /; 1 <= n <= Length[ai] := ai[[n]]

In[8]:= Block[{aif, zif}, 
 RecurrenceTable[{A[0] == 0, A[1] == aif[1], 
   A[n + 1] == A[n] + (z - zif[n]) aif[n + 1] A[n - 1]}, 
  A, {n, 0, (Length@ai) - 1}]]

Out[8]= {0, 0.904837, 0.904837}

Block serves to temporarily remove definitions of aif and zif while RecurrenceTable is executed. Then, when Block exits, the values are restored, and the output of RecurrenceTable evaluates.

Sasha
  • 5,935
  • 1
  • 25
  • 33
  • 1
    @Sasha: this is probably stupid, but shouldn't it be aif[3] rather than aif[4] above? and where is A[0]? – Phil Apr 26 '11 at 17:00
  • 1
    @Phi, you are right. Clear definitions of `zif` and `aif` and rerun `RecurrenceTable`. This is the output I would expect. But if we redefine `aif` and `zif`, and rerun the `RecurrenceTable` the same odd output is obtains. Looks like a bug. Good catch! – Sasha Apr 26 '11 at 17:14
  • @Sasha, I'm still getting the error. Nice attempt, but it doesn't seem to like me. – rcollyer Apr 26 '11 at 17:17
  • @rcollyer, I am getting the same result as @Sasha with no error message on Mathematica v.8 but it's only that the result is odd as noted above. @sasha clearing definitions does not change the output for me – Phil Apr 26 '11 at 17:24
  • @Phi You also need to replace `A[1]==ai[[1]]` with `A[1]==aif[1]`, clear definitions and then rerun `RecurrenceTable` – Sasha Apr 26 '11 at 17:29
  • Correction, if you change the order of evaluation where the `tmp = RecurrenceTable[...]` is evaluated before `aif` and `zif` are defined, you get something. But, the error pops back up when I use it in `Plot`, i.e. it is only masked – rcollyer Apr 26 '11 at 17:32
  • @sasha: cleared everything, no change. It seems to me RecurrenceTable has to be further debugged on the mathematica development side – Phil Apr 26 '11 at 17:35
  • @Phil, at one point, this worked for me. But, then I changed the code to allow for another scenario - which never worked, anyway. However, when I attempted to return the code to the working form, it no longer worked. I don't know what I changed. I'll post the full code above and see if more eyes can tell me what went wrong. – rcollyer Apr 26 '11 at 18:01
3

It seems to me that Sasha's approach can be mimicked by just Blocking Part.

zi = {0.1, 0.2, 0.3};
ai = {0.904837, 1.05171, -0.499584};

Block[{Part},
 RecurrenceTable[{A[0] == 0, A[1] == ai[[1]], 
   A[n + 1] == A[n] + (z - zi[[n]]) ai[[n + 1]] A[n - 1]}, 
  A, {n, Length@ai - 1}]
]
   {0, 0.904837, 0.904837} 

Addressing Sasha's criticism, here are two other ways one might approach this:

With[{Part = $z}, 
  RecurrenceTable[{A[0] == 0, A[1] == ai[[1]], 
    A[n + 1] == A[n] + (z - zi[[n]]) ai[[n + 1]] A[n - 1]}, 
   A, {n, Length@ai - 1}]
] /. $z -> Part

-

With[{Part = Hold[Part]}, 
  RecurrenceTable[{A[0] == 0, A[1] == ai[[1]], 
    A[n + 1] == A[n] + (z - zi[[n]]) ai[[n + 1]] A[n - 1]}, 
   A, {n, Length@ai - 1}]
] // ReleaseHold
Mr.Wizard
  • 24,179
  • 5
  • 44
  • 125
  • +1, it does, and it's simpler. I'll consider changing the accepted answer. – rcollyer Apr 26 '11 at 19:57
  • 2
    I would be cautious here. You do not know whether `Part` is ever used inside `RecurrenceTable` code. Doing what you suggest might work here, but is a good generic method to invite troubles. I do not mind rcollyer to revert the accept answer, just warning of possible dire consequences of disabling `Part` while executing system commands. – Sasha Apr 26 '11 at 21:07
  • @Sasha, good point. If RecurrenceTable were a kernel function, could this be presumed safe? – Mr.Wizard Apr 26 '11 at 21:12
  • I would say no, but I do not know for sure. – Sasha Apr 26 '11 at 21:23