> # Rungebeispiel
> 
> restart;
> with(linalg):
> Newtoninterpol:=proc(x,f,z)
> # f ist der Vektor der Stuetzwerte
> # x ist der Vektor der Stuetzstellen
> # z ist der Vektor der Stuetzstellen, fuer den ein Interpolationswert gesucht ist
> local Newtonkoeff,y,p,m,n,j,k;
> Newtonkoeff:=proc(x,f)
> #x ist der Vektor der Stuetzstellen
> #f ist der Vektor der Stuetzwerte
> local y,n,i,k;
> n:=vectdim(x);
> y:=vector(n);
> for i from 1 to n do y[i]:=f[i] od;
> for k from 1 to n-1 do 
>    for i from n by -1 to k+1 do
>       y[i]:=evalf((y[i]-y[i-1])/(x[i]-x[i-k]))
>    od
> od;
> evalm(y)
> end;
> # Berechnung der Funktionswerte des Interpolationspolynoms
> y:=Newtonkoeff(x,f);
> n:=vectdim(y);
> m:=vectdim(z);
> p:=vector(m);
> for j from 1 to m do
>     p[j]:=y[n];
>     for k from n-1 by (-1) to 1 do p[j]:=evalf(y[k]+(z[j]-x[k])*p[j]) od
> od;
> evalm(p)
> end:
> 
Warning, new definition for norm
Warning, new definition for trace
> m:=201:
> z:=vector(m):
> a:=-5: b:=5:
> for j from 1 to m do z[j]:=evalf(a+(j-1)*(b-a)/(m-1)) od:
> n:=11; x:=vector(n): f:=vector(n): x1:=vector(n):
> for j from 1 to n do x[j]:=evalf(a+(j-1)*(b-a)/(n-1)) od:
> for j from 1 to n do x1[j]:=evalf(-(cos((2*j-1)*Pi/(2*n))*(b-a)+a+b)/2) od:
> print(`äquidistante Knoten`,x);
> print(`Tschebyscheff-Knoten`,x1);
> for j from 1 to n do f[j]:=evalf(1/(x[j]^2+1)) od:
> for j from 1 to n do f1[j]:=evalf(1/(x1[j]^2+1)) od:
> 

                               n := 11


äquidistante Knoten, [-5., -4., -3., -2., -1., 0, 1., 2., 3., 4., 5.]


Tschebyscheff-Knoten, [-4.949107210, -4.548159977, -3.778747872,

    -2.703204087, -1.408662783, 0, 1.408662783, 2.703204087,

    3.778747872, 4.548159977, 4.949107210]

> p:=Newtoninterpol(x,f,z):
> p1:=Newtoninterpol(x1,f1,z):
> with(plots):
> F:=plot(1/(x^2+1), x=-5..5, y=-0.5..2, thickness=3): 
> G1:=listplot([seq([z[t],p[t]],t=1..m)],thickness=0):
> G2:=listplot([seq([z[t],p1[t]],t=1..m)],thickness=2):
> display({F,G1,G2});

>