kill(all)$

program : Lagrange interpolation polynomial (version 1.0)
peter.vlasschaert@gmail.com,31/05/2019

load :packages.  (we load other packages later , "see  below")
==============

load(lrats)$

How to use : reverse command, 'Maxima'.
===================================

x^3+2*x+1=10*x;

reverse(%);

%/10;

info : "Lagrange interpolation polynomial".
====
----------------------------------------------------------------------------------------
see :  https://www.youtube.com/watch?v=_zK_KhHW6og
see :  https://en.wikipedia.org/wiki/Lagrange_polynomial
----------------------------------------------------------------------------------------

(n)umber (p)oints : np
==================
every  point : (x[i],y[i])    , i=1,2,....np

np:3;

reason : combination 'sum & product '
================================

why :  som         ,   init value = 0     =>  suminit  = 0.
             product  ,   init value = 1     =>  prodinit = 1.

suminit:0$
for i:1 thru  np do
     (productinit:1,
          for j:1 thru np do
               if j≠ i then
                   productinit : productinit*(x-x[j])/(x[i]-x[j]),
   (suminit :suminit + productinit*y[i])
       );

rem  :  only result when  i≠j otherwise division by zero.

output :
=======

suminit;

procductinit;

to make function
==============

1e) method     'define'
==========

p1:define(F(x),suminit);

2e) method   ' ''   '
==========

p2:F(x):=''suminit;

function , mylagrange.
===================

kill(suminit);

suminit;

mylagrange(nn):= block([],
    suminit:0,
for i:1 thru  nn do
     (productinit:1,
          for j:1 thru nn do
               if j≠ i then
                   productinit : productinit*(x-x[j])/(x[i]-x[j]),
   (suminit :suminit + productinit*y[i])
       ),return(suminit)
                                        )$
    

mylagrange(4);

p3a:suminit;

p3b:productinit;

list of 'np' points ' x[i],y[i] ', i=1,2,..,k,....np
==================================

xylist:makelist([x[k],y[k]],k,1,np);

How to get same result with a package.
=================================

load(interpol);

p3:lagrange(xylist);

rem: 'p3'  = rhs 'p2'
-----------------------------

example: p1:[-2,2],p2:[1,4],p3:[2,5]
=============================

p4a1:[-2,2];
p4a2:[1,4];
p4a3:[2,5];

p4:[p4a1,p4a2,p4a3];

F(x) : ?
======

p5:lagrange(p4),expand;

wxplot2d  , 'discrete'
=================

wxplot2d([p5,['discrete,p4]],[x,-2.5,2.5],[y,-10,10]);

wxplot2d([p5,['discrete,p4]],[x,-2.5,2.5],[y,-10,10],
    [style,[lines,1,2],[points,3,1,1]],[legend,"F(x)","p1,p2,p3"]);