Sie sind hier: Numerik > Numerische Integration > Adaptive Integration mit Trapezregel

Adaptive Integration mit Trapezregel

Donnerstag 07. Dezember 2006 von
Simon Praetorius
Anzeige einer Funktion mit GitterstellenAnzeige einer Funktion mit Gitterstellen
Um ein bestimmtes Integral in den Grenzen a und b (a<b) näherungsweise zu berechnen, gibt es mehrere Verfahren. Ich möchte hier die mit der sogenannten Trapezregel implementieren. Die Integration soll adaptiv gelöst werden, d.h. die Intervalle, auf denen der (mittels Vergleich mit der Simpson-Regel) geschätzte Fehler größer ist als die vorgegebene Toleranz, werden halbiert und dann die Rechnung rekursiv auf den beiden Teilintervallen wiederholt.

matlab Code
  • function [Int I]=trapez(f,I,tol)
  • % [Integral Intervall] = trapez(@func_name,Intervall,toleranz)
  • % Funktion zur Berechnung des Integrals
  • % einer uebergebenen Funktion f
  • % in den Grenzen I(1) bis I(2)
  • % mittels Trapezregel und einer
  • % Toleranz von tol=(10^(-k)) / (I(2)-I(1))
  • % Die Funktion gibt den berechneten Integralwert
  • % und den Vektor der Gitterstruktur zurueck
  • Int = (I(2)-I(1))/2 * (feval(f,I(1)) + feval(f,I(2)));
  • S = (I(2)-I(1))/6 * (feval(f,I(1)) + 4*feval(f,(I(2)+I(1))/2) + feval(f,I(2)));
  • if(abs(Int-S)>=tol*(I(2)-I(1)))
  • I = [I(1);(I(2)+I(1))/2;I(2)];
  • [Int1 I1] = trapez(f,I(1:2),tol);
  • [Int2 I2] = trapez(f,I(2:3),tol);
  • Int = Int1 + Int2;
  • I = [I1;I2(2:end)];
  • end
  • end
Ich habe das Verfahren so implementiert, dass man die Intervall-Struktur im Nachhinein ablesen kann.

Zum Testen der Implementierung habe ich eine weitere Funktion geschrieben, so, dass auch x als Vektor eingegeben werden kann. Dies wird benötigt, um die Funktion zu plotten.
matlab Code
  • function y=f1(x)
  • % Exponential-Funktion y=f(x)=e^(-8x^2)
  • y=exp(-8.*(x.^2));
  • end
Um das ganze zu testen, kann man folgendermaßen vorgehen:
matlab Code
  • [Int I] = trapez(@f1,[0 2],(10^(-3))/2);
  • Int
  • x=linspace(0,2,100);
  • y=f1(x);
  • y2=f1(I);
  • plot(x,y,I,y2,'o')
Dies erzeugt das oben eingefügte Bild.
Besucher: 15101 | Permalink | Kategorie: Numerik
Tags: , , ,

Kommentar hinzufügen

Dieses Feld bitten nicht ausfüllen: