Sie sind hier: Numerik > Interpolationsaufgaben > Interpolation durch kubischen Spline

Interpolation durch kubischen Spline

Dienstag 14. November 2006 von Simon Praetorius
Spline der Runge-Funktion mit 8 StützstellenSpline der Runge-Funktion mit 8 Stützstellen
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

Spline der Runge-Funktion mit 16 StützstellenSpline der Runge-Funktion mit 16 Stützstellen
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.

scilab Code
  • // 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.");

Kommentar hinzufügen

Dieses Feld bitten nicht ausfüllen: