Die Testfunktion für die Interpolation ist eine abgewandelte Form der Runge-Funktion 1 / (1+x2), die ein Besipiel für das nach Carl Runge benannte Runge-Phänomen darstellt: Interpolation der Funktion durch einfache Polynominterpolation und Erhöhung des Polynomgrades kann zu einer Verschlechterung der Interpolationsgüte führen (siehe Wikipedia).

• function y = runge(x)
• y = 1./(1+25.*x.*x);

Zum Testen der Spline-Interpolation habe ich zwei Funktionen geschrieben. Die erst (test_spline) Führt zu einer Gegebenen Zerlegung eines Intervalls I eine Spline-Interpolation mit kubischen Splines durch. Die zweite Funktion (animation_spline) stellt Interpolationen zu verschiedenen Zerlegungen des Intervalls als einfach Animation dar.

• function erg = c_spline(z,x,M)
• % Zur Berechnung des/der Funktionswerte(s) des IP an der/den Stelle(n) {z_i}
• % s({z_i}) = c_spline({z_i},{x_i,f_i},{M_i})
•  
• if nargin < 3
• error('Zahl der Aufrufargumente stimmt nicht')
• end
•  
• m=length(z);
• n=length(x(:,1));
• if(n~=length(M))
• 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)
• error('Stuetzstellen muessen paarweise verschieden sein!');
• end
• end
•  
• erg = zeros(m);
• for k=1:m
• if(z(k)<x(n,1))
• 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))
• erg(k) = x(n,2);
• else
• error('Eine zu berechnende Stelle liegt außerhalb des Definitionsbereiches des Splines.'));
• end
• end
• end

• function M=momente(x)
• % Zur Berechnung der Momente M_k = s''(x_k) eines C^2-Splines
• % M_k = momente({x_i,f_i})
•  
• if nargin < 1
• error('Zahl der Aufrufargumente stimmt nicht')
• end
•  
• 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)
• 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); A(i,i-1) = m; end
• if(i<n-2); 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
•  
• if(n>2)
• N = Ab;
• M = [0;N(:,1);0];
• elseif n==2
• M = [0;0];
• else
• error('Zu wenig Stutzstellen gegeben. Eine Interplation ist nicht moeglich.');
• end
•  
• end

• function test_spline(f,I,l)
• % Zum Testen der SIP einer geg. Funktionen f im Intervall I, Feinheit l.
• % test_spline(@f,I,l)
•  
• x(:,1) = linspace(I(1),I(2),l);
• x(:,2) = f(x(:,1));
•  
•  
• z = linspace(min(x(:,1)),max(x(:,1)),(I(2)-I(1))*50);
•  
• M = momente(x);
• s = c_spline(z,x,M);
• y = f(z);
•  
• [ny nx] = size(y);
• if(nx~=1); y=y'; end
•  
• plot(x(:,1),x(:,2),'o',z,s,z,y)
• legend('Stuetzpunkte','C-Spline','Funktion')
• end

• function animation_spline(f,I,l)
• % Funktion zur animierten Anzeige einer Spline-IP
• % durch stetige incrementierung der Stützstellenzahl
• % an
• % aniimation_spline(@function_name,[xmin xmax],l)
• % l = Anzahl Stuetzstellen+1
•  
• if nargin < 3
• error('Zahl der Aufrufargumente stimmt nicht')
• end
•  
• for i=4:l
• test_spline(@f,I,i)
• pause(0.1)
• end
• end