> restart;
> with(linalg): with(plots):
Warning, new definition for norm
Warning, new definition for trace
> kubspline:=proc(a,h,x)
> local f;
> # option trace;
> if (x>=a) then if (x<a+h) then f:=evalf((x-a)^3/(6*h^3)) fi fi;
> if (x>=a+h) then if (x<a+2*h) then f:=evalf((h^3+3*h^2*(x-a-h)+3*h*(x-a-h)^2-3*(x-a-h)^3)/(6*h^3)) fi fi;
> if (x>=a+2*h) then if (x<a+3*h) then f:=evalf((h^3+3*h^2*(a+3*h-x)+3*h*(x-a-3*h)^2-3*(a+3*h-x)^3)/(6*h^3)) fi fi;
> if (x>=a+3*h) then if (x<a+4*h) then f:=evalf((a+4*h-x)^3/(6*h^3)) fi fi;
> if (x<a) then f:=0 fi;
> if (x>=a+4*h) then  f:=0 fi;
> RETURN(f)
> end:
> m:=100:
>  a:=0: h:=1:
> for j from 1 to m do z[j]:=evalf(a+(j-1)*(4*h)/(m-1)) od:
> for j from 1 to m do spl[j]:=kubspline(a,h,a+(j-1)*4*h/(m-1)) od:
> listplot([seq([z[t],spl[t]],t=1..m)],scaling =CONSTRAINED);

> # Interpolation mit Hermite-Endbedingungen
> kubInterpol:=proc(x,y,z,f0,f1)
> # x ist die Vektor der äquidistanten Stützstellen
> # y ist der Vektor der Stützwerte
> # f0 = f'(x1), f1= f'(x_n) (Hermite-Endbedingungen)
> # z ist der Vektor der Stützstellen, an denen der Spline-Interpolant ausgewertet wird
> # option trace;
> global s;
> # s ist der Vektor der Werte des Spline-Interpolants an Stützstellen z;
> local n, kubspline,j,k,m,a,anf,h,f,koeffmatrix;
> kubspline:=proc(a,h,x)
> local f;
> # option trace;
> if (x>=a) then if (x<a+h) then f:=evalf((x-a)^3/(6*h^3)) fi fi;
> if (x>=a+h) then if (x<a+2*h) then f:=evalf((h^3+3*h^2*(x-a-h)+3*h*(x-a-h)^2-3*(x-a-h)^3)/(6*h^3)) fi fi;
> if (x>=a+2*h) then if (x<a+3*h) then f:=evalf((h^3+3*h^2*(a+3*h-x)+3*h*(x-a-3*h)^2-3*(a+3*h-x)^3)/(6*h^3)) fi fi;
> if (x>=a+3*h) then if (x<a+4*h) then f:=evalf((a+4*h-x)^3/(6*h^3)) fi fi;
> if (x<a) then f:=0 fi;
> if (x>=a+4*h) then  f:=0 fi;
> RETURN(f)
> end;
> n:=vectdim(x);
> f:=vector(n+2);
> for j from 1 to n do f[j+1]:=y[j] od;  f[1]:=f0; f[n+2]:=f1;
> koeffmatrix:=matrix(n+2,n+2,0);
> for j from 2 to n+1 do koeffmatrix[j,j]:=2/3; 
>  koeffmatrix[j,j-1]:=1/6; koeffmatrix[j,j+1]:=1/6 
> od;
> h:=x[2]-x[1];
> koeffmatrix[1,1]:=-1/(2*h);    koeffmatrix[1,3]:=1/(2*h); 
> koeffmatrix[n+2,n+2]:=1/(2*h);  koeffmatrix[n+2,n]:=-1/(2*h); 
> a:=vector(n+2);
> # Berechnung der Lösung des LGS koeffmatrix * a=f
> a:=linsolve(koeffmatrix, f);
> #Berechnung der Funktionswerte des Interpolationssplines
> m:=vectdim(z);
> for j from 1 to m do 
> s[j]:=0; anf:=x[1];
> for k from 1 to n+2 do s[j]:=s[j]+a[k]*kubspline(anf+(k-4)*h ,h,z[j]) od
> od;
> RETURN(s)
> end:
> x:=vector(11,[-5,-4,-3,-2,-1,0,1,2,3,4,5]):
> y:=vector(11):
> for j from 1 to 11 do y[j]:=evalf(sin(2*x[j]+1)+cos(x[j]-3/2)) od:
> m:=101:
> z:= vector(m):
> a:=-5: b:=5: h:=(b-a)/(m-1):
> for j from 1 to m do z[j]:=a+ (j-1)*h od:
> f0:=evalf(2*cos(2*(-5)+1)-sin(-5-3/2)): 
> f1:=evalf(2*cos(2*(5)+1)-sin(5-3/2)):
> kubInterpol(x,y,z,f0,f1): with(plots):
> G1:=listplot([seq([z[t],s[t]],t=1..m)]):
> F:=plot(sin(2*x+1)+cos(x-3/2), x=-5..5):
> display({F,G1});

>