Sie sind hier: Numerik > Interpolationsaufgaben > Interpolation durch kubischen Spline > MatLab Variante zur Interpolation ...

MatLab Variante zur Interpolation mit kubischen Splines

Dienstag 14. November 2006 von
Simon Praetorius
Eine abgewandelte Form der Runge-FunktionEine abgewandelte Form der Runge-Funktion
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).
matlab Code
  • 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.
matlab Code
  • 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)
  • 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 auerhalb des Definitionsbereiches des Splines.'));
  • end
  • end
  • end
matlab Code
  • 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(iend
  •  
  • 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
matlab Code
  • 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
matlab Code
  • function animation_spline(f,I,l)
  • % Funktion zur animierten Anzeige einer Spline-IP
  • % durch stetige incrementierung der Sttzstellenzahl
  • % 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
Besucher: 53590 | Permalink | Kategorie: Numerik
Tags: , ,

Kommentar hinzufügen

Dieses Feld bitten nicht ausfüllen: