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.

• 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.

• 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:

• [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.