> #aufgabe 23
> restart;
> with(linalg): with(plots):
Warning, new definition for norm
Warning, new definition for trace
> # kubische Spline-Interpolation mit natürlichen Endbedingungen
> kubInterpol:=proc(x,y,z)
> # x ist die Vektor der äquidistanten Stützstellen
> # y ist der Vektor der Stützwerte
> # 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]:=0; f[n+2]:=0;
> 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/(h^2);   koeffmatrix[1,2]:=-2/(h^2);  koeffmatrix[1,3]:=1/(h^2); 
> koeffmatrix[n+2,n+2]:=1/(h^2);  koeffmatrix[n+2,n+1]:=-2/(h^2); koeffmatrix[n+2,n]:=1/(h^2); 
> a:=vector(n+2);
> a:=multiply(inverse(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(1/(25+x[j]^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:
> kubInterpol(x,y,z): 
> G1:=listplot([seq([z[t],s[t]],t=1..m)]):
> F:=plot(1/(25+x^2), x=-5..5):
> display({F,G1});

>