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 ai
s from the zi
s 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
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.
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]}
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.
It seems to me that Sasha's approach can be mimicked by just Block
ing 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
If you love us? You can donate to us via Paypal or buy me a coffee so we can maintain and grow! Thank you!
Donate Us With