Hello,
I suppose I'm late to the party, but I was reading some notes about
polynomial interpolation and remembered having seen this thread. My
contribution is a version using the new fold operator:
V =: {{ ((}.-}:) y) % ((-x){.m) - x{.m }} NB. iterate divided differences
A =: {{ ({.y) , y {.F::(x V) i.&.<:#y }} NB. newton coefficients
N =: {{ (({."1 u)A({:"1 u)) +/ . * 1,*/\ y-}:{."1 u }} NB. evaluate
polynomial at y
({:"1 -: (zs N"0)@{."1) zs =: (,. 3&o.) _3r2 _3r4 0 3r4 3r2
pd 'reset; type line'
pd (; zs N"0) 1.51 * (%~ i:) 1000
pd (; 3&o.) 1.51 * (%~ i:) 1000
pd 'show'
Joseph
On 10/19/20, 'Mike Day' via Programming <[email protected]> wrote:
> Funny, I think I’ve usually done the quicker one.
>
> Thanks for reminding me,
>
> Mike
>
> Sent from my iPad
>
>> On 19 Oct 2020, at 22:03, Raul Miller <[email protected]> wrote:
>>
>> A minor side note:
>>
>> (+/\@:= -: >:/~) i.2000
>> 1
>> timespacex '+/\=i.2000'
>> 0.010131 3.78202e7
>> timespacex '>:/~i.2000'
>> 0.002597 4.21203e6
>>
>> Thanks,
>>
>> --
>> Raul
>>
>> On Thu, Oct 15, 2020 at 2:21 PM 'Michael Day' via Programming
>> <[email protected]> wrote:
>>>
>>> It took me a while to see what was required using the matrix, M and the
>>> derived coefficients, a.
>>>
>>> But this works ok. I haven't turned the explicit code into a function,
>>> but that should be quite
>>> straightforward.
>>>
>>> NB. fixed point font best past here....
>>> a =. y %. M =. 1,. }:"1 */\"1 (-/~ * +/\@:=@i.@#) x
>>> M
>>> 1 0 0 0 0
>>> 1 3r4 0 0 0
>>> 1 3r2 9r8 0 0
>>> 1 9r4 27r8 81r32 0
>>> 1 3 27r4 81r8 243r32
>>> a
>>> _14.10142 17.559765 _10.878424 4.8348551 0
>>>
>>> 6{."1 [n =: 1,}:*/\ x-~/steps _1.5 1.5 30 NB. n is Newton basis
>>> form of new x
>>> 1 1 1 1 1 1
>>> 0 0.1 0.2 0.3 0.4 0.5
>>> 0 _0.065 _0.11 _0.135 _0.14 _0.125
>>> 0 0.091 0.143 0.162 0.154 0.125
>>> 0 _0.19565 _0.29315 _0.3159 _0.2849 _0.21875
>>>
>>> plot a +/ . * n
>>>
>>> Earlier, I'd overlooked the (rather obvious!) fact that you need to
>>> apply a to the Newton
>>> basis form, n here, of your required prediction points, not the
>>> untransformed points.
>>>
>>> Mike
>>>
>>>
>>>> On 15/10/2020 14:50, Stefan Baumann wrote:
>>>> Dear all.
>>>> Recently I stumbled upon the Newton polynomial and took it as a practice
>>>> to
>>>> implement it without using loops, but failed. I first didn't get a grab
>>>> on
>>>> how to create the matrix used in
>>>> https://en.wikipedia.org/wiki/Newton_polynomial#Main_idea, and
>>>> eventually
>>>> came up with code following
>>>> https://en.wikipedia.org/wiki/Newton_polynomial#Application:
>>>>
>>>> NB. Newton polynomial
>>>>
>>>> np=: 4 : 0
>>>>
>>>> a=. {. y
>>>>
>>>> for_i. }.>:i.#x do.
>>>>
>>>> y=. (2 (-~)/\ y) % i ({:-{.)\ x NB. Divided differences
>>>>
>>>> a=. a, {. y NB. Coefficients are the topmost entries
>>>>
>>>> end.
>>>>
>>>> NB. Convert the summands aᵢ(x-x₀)…(x-xᵢ₋₁) of the polynomial
>>>>
>>>> NB. from multiplier-and-roots to coefficients form and add them up
>>>>
>>>> +/@,:/ p."1 (;/a) ,. (<''), }:<\ x
>>>>
>>>> )
>>>>
>>>> x=: _3r2 _3r4 0 3r4 3r2
>>>>
>>>> y=: 3&o. x
>>>>
>>>> load'plot'
>>>>
>>>> load'stats'
>>>>
>>>> plot (];(x np y)&p.) steps _1.5 1.5 30
>>>>
>>>> I also tried replacing the loop with fold F:. but again was not able to
>>>> do
>>>> so. Anyone out there who can enlighten me?
>>>>
>>>> Many thanks. Stefan.
>>>> ----------------------------------------------------------------------
>>>> For information about J forums see http://www.jsoftware.com/forums.htm
>>>
>>>
>>> --
>>> This email has been checked for viruses by Avast antivirus software.
>>> https://www.avast.com/antivirus
>>>
>>> ----------------------------------------------------------------------
>>> For information about J forums see http://www.jsoftware.com/forums.htm
>> ----------------------------------------------------------------------
>> For information about J forums see http://www.jsoftware.com/forums.htm
> ----------------------------------------------------------------------
> For information about J forums see http://www.jsoftware.com/forums.htm
>
----------------------------------------------------------------------
For information about J forums see http://www.jsoftware.com/forums.htm
