Splines sind Funktionen, die stückweise aus Teilfunktionen, insbesondere aus Polynomen, zusammengesetzt sind. Ein kubischer C2-Spline ist eine Art Standardspline. Dieser besteht aus einer Verbindung von Polynomen 3.Grades deren Verbindungsstellen zwei mal stetig differenzierbar sind.
Will man nun eine Funktion f durch einen solchen Spline s interpolieren, legt man in jedes Intervall [xk, xk+1], k=0,...,n-1 ein Polynom pk (wie oben). Durch die Interpolationsbedingungen f(xk)=pk(xk) und f(xk+1)=pk(xk+1) und die Glattheitsbedingungen p'k(xk+1)=p'k+1(xk+1) bzw. p''k(xk+1)=p''k+1(xk+1) ist der Spline noch nicht eindeutig bestimmt. Man hat noch 2 Freiheitsgrade. Diese werden durch zusätzliche Randbedingungen, z.B. den natürlichen Randbedingungen s''(x0) = 0 = s''(xn) erfüllt.
Bei der Implementierung der Splineinterpolation nutzt man die Momente Mk := s''(xk) der Funktion und stellt damit ein tridiagonales Gleichungssystem auf, das dann mit entsprechenden Algorithmen gelöst werden kann.
Eine nähere Erläuterung zu kubischen Splines findet man z.B. auf www.arndt-bruenner.de/mathe/scripts/kubspline.htm

Es soll zusätzlich eine Funktion geschrieben werden, die den interpolierten Wert an einer beliebigen Stelle z berechnet. Dann soll die Spline-Interpolation grafisch dargestellt werden, durch den Vergleich einer Funktion mit dem interpolierten Polynomen.

• // Programm zur Berechnung der Moment M_k = s''(x_k) des kubischen C^2-Splines s mit natuerlichen Randbedingungen
•  
• function erg=spline(z,x,M)
• m=length(z);
• n=length(x(:,1));
• if(n<>length(M)) then; error("Falsche Werte uebergeben - Dimensionen von Momenten- und Stuetzstellenvektor stimmen nicht ueberein."); end
•  
• h = zeros(n);
• for i=1:n-1
• h(i+1)=x(i+1,1)-x(i,1);
• if(h(i+1)==0) then; error("Stuetzstellen muessen paarweise verschieden sein!"); end
• end
•  
• erg = zeros(m);
• for k=1:m
• if(z(k)<x(n,1)) then
• j = max(find(x(:,1)<=z(k)));
• erg(k) = M(j)*(x(j+1,1)-z(k))^3 / (6*h(j+1)) + M(j+1)*(z(k)-x(j,1))^3 / (6*h(j+1)) + ((x(j+1,2)-x(j,2))/h(j+1) - (M(j+1)-M(j))*h(j+1)/6)*(z(k)-x(j,1)) + x(j,2)-M(j)*(h(j+1)^2)/6;
• elseif(z(k)==x(n,1)) then
• erg(k) = x(n,2);
• else
• error("Eine zu berechnende Stelle liegt ausserhalb des Definitionsbereiches des Splines.");
• end
• end
• endfunction
•  
• function M=momente(x)
• n = length(x(:,1));
•  
• A = eye(n-2,n-2)*2;
• b = zeros(n-2);
• h = zeros(n-1);
•  
• for i=1:n-1
• h(i+1)=x(i+1,1)-x(i,1);
• if(h(i+1)==0) then; error("Stuetzstellen muessen paarweise verschieden sein!"); end
• end
•  
• for i=1:n-2
• m = h(i+1) / (h(i+1)+h(i+2));
• if(i>1) then; A(i,i-1) = m; end
• if(i<n-2) then; A(i,i+1) = 1-m; end
•  
• b(i) = 6/(h(i+1)+h(i+2))*( (x(i+2,2)-x(i+1,2))/h(i+2) - (x(i+1,2)-x(i,2))/h(i+1) );
• end
•  
• M = [0;Ab;0];
• endfunction
•  
• function y=runge(x)
• for j=1:length(x)
• y(j)=1/(1+25*x(j)^2);
• end
• endfunction
•  
• function test_spline(f,I,l)
• x=zeros(l,2);
• for i=0:l
• a=I(1)+i*((I(2)-I(1))/l);
• b=f(a);
• x(i+1,:) = [a b];
• end
•  
• z = linspace(min(x(:,1)),max(x(:,1)),100);
•  
• M = momente(x);
• s = spline(z,x,M);
• y = f(z);
•  
• [ny nx] = size(y);
• if(nx<>1) then; y=y'; end
•  
• ymin = min(y);
• ymax = max(y);
• dy = (ymax-ymin)*0.05;
• rect = [I(1),ymin-dy,I(2),ymax+dy];
•  
• xbasc()
• plot2d(x(:,1),x(:,2),-9,"011","",rect)
• plot2d(z,[s y], [2 3],"100", leg="Spline@Runge-Funktion@y=exp(x)")
• xset("wshow")
• endfunction
•  
• function sleep(s)
• x=0;
• for i=1:0.05:s
• for j=1:0.05:s
• x=x+1;
• end
• end
• endfunction
•  
• function animation_spline(f,I,l)
• driver("X11")
• xset("pixmap",1)
• xtitle("Spline-Interpolation")
• for i=1:l
• test_spline(f,I,i)
• sleep(15)
• end
• xset("pixmap",0)
• driver("Rec")
• endfunction
•  
• disp("Programm enthaelt folgende Funktionen:");
• disp(" M_k = momente({x_i,f_i}) ................... Zur Berechnung der Momente M_k = s""(x_k) eines C^2-Splines");
• disp(" s({z_i}) = spline({z_i},{x_i,f_i},{M_i}) ... Zur Berechnung des/der Funktionswerte(s) des IP an der/den Stelle(n) {z_i}");
• disp(" y=runge({x_i}) ............................. Implementierung der Runge-Funktion");
• disp(" test_spline(f,I,l).......................... Zum Testen der SIP einer geg. Funktionen f im Intervall I, Feinheit l.");