> #rationale Interpolation
> restart;
> with(linalg):
Warning, new definition for norm
Warning, new definition for trace
> inverseDiff:=proc(x,f)
> #x ist der Vektor der Stuetzstellen
> #f ist der Vektor der Stuetzwerte
> local y,n,i,k,h;
> n:=vectdim(x);
> y:=vector(n);
> for i from 1 to n do 
>    h:=f[i];   
>    for k from 2 to i do
>       h:=evalf((x[i]-x[k-1])/(h-y[k-1]))
>    od;
>    y[i]:=h
> od;
> evalm(y)
> end:
> 
> x:=vector(4,[0,1,2,3]):
> y:=vector(4,[0,-1,-2/3,9]):
> w:=inverseDiff(x,y);

               w := [0, -1., -.5000000000, .5000000000]

> rationalinterpol:=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 inverseDiff,y,r,m,n,j,k;
> inverseDiff:=proc(x,f)
> #x ist der Vektor der Stuetzstellen
> #f ist der Vektor der Stuetzwerte
> local y,n,i,k,h;
> n:=vectdim(x);
> y:=vector(n);
> for i from 1 to n do 
>    h:=f[i];   
>    for k from 2 to i do
>       h:=evalf((x[i]-x[k-1])/(h-y[k-1]))
>    od;
>    y[i]:=h
> od;
> evalm(y)
> end:;
> # Berechnung der Funktionswerte des Interpolationspolynoms
> y:=inverseDiff(x,f);
> n:=vectdim(y);
> m:=vectdim(z);
> r:=vector(m);
> for j from 1 to m do
>     r[j]:=0;
>     for k from n-1 by (-1) to 1 do r[j]:=evalf((z[j]-x[k])/(y[k+1]+r[j])) od;
>     r[j]:=r[j]+y[1]
> od;
> evalm(r)
> end:
> 
> n:=100: z:=vector(n):
> a:=0: b:=3:
> for j from 1 to n do z[j]:=evalf(a+(j-1)*(b-a)/(n-1)) od:
> p:=rationalinterpol(x,y,z):
> with(plots):
> listplot([seq([z[t],p[t]],t=1..n)]);

>