interpolacni.html

procedura Newton()

> Newton := proc(k,body)
        global N;
        local A,B,a,i,j,n,m,y,f;

        with(linalg):
        A:=array(1..k,1..k):
        for i from 1 by 1 to k do
            A[i,1]:=1:
        end do:
        for i from 1 by 1 to k do
            for j from 2 by 1 to k do
                A[i,j]:=0:
            end do:
        end do:
        f(y):=1:
        for j from 2 by 1 to k do
            f(y):=f(y)*(y-body[j-1][1]):
            for i from j by 1 to k do
                A[i,j]:=eval(f(y),y=body[i][1]):
            end do:
        end do:
        B:=array(1..k):
        for i from 1 by 1 to k do
            B[i]:=body[i][2]:
        end do:
        a:=linsolve(A,B):
        N(x):=sum(a[n]*product((x-body[m][1]),m=1..n-1), n=1..k);
        N(x):=simplify(N(x));
    end: