Lösung von Aufgabe 17

• Die Aufgabe lässt sich mit etwas Handarbeit erledigen, indem man den Integranden und die Berechnung der Simpson-Näherungen I₁ und I₂ jeweils als Matlab-Funktionen definiert und per Hand iteriert.
  • Man erhält folgende Schritte:

    •    a      b       I1      I2
      ===========================================
       0.400  4.000   2.7244  1.4803   halbieren
       0.400  2.200  -0.2736 -1.1846   halbieren
       0.400  1.300  -0.6431 -0.3510   halbieren
       0.400  0.850   0.0275  0.0273   ok
       0.850  1.300  -0.3786 -0.3995   halbieren
       0.850  1.075  -0.1988 -0.1996   ok
       1.075  1.300  -0.2007 -0.2001   ok
       1.300  2.200  -0.5415  0.1617   halbieren
       1.300  1.750   0.1498  0.0699   halbieren
       1.300  1.525   0.1156  0.1137   ok
       1.525  1.750  -0.0457 -0.0435   ok
       1.750  2.200   0.0120 -0.0064   halbieren
       1.750  1.975  -0.1296 -0.1302   ok
       1.975  2.200   0.1232  0.1230   ok
       2.200  4.000   1.7540  1.7764   halbieren
       2.200  3.100   0.8764  0.8839   ok
       3.100  4.000   0.9000  0.9000   ok

  • benutzte x-Werte sind also
    • 0.400 0.850 1.075 1.300 1.525 1.750 1.975 2.200 3.100 4.000
  • im Bild
    • 
  • Für das Integral ergibt sich damit
    • I ≈ 1.4750 (genauer Wert: 1.4738)
• Mit etwas mehr Matlab-Programmierung lässt sich das Problem sehr schnell und elegant durch Rekursion lösen:
  • Funktion simpson_adapt berechnet die beiden Integrale (mit einem Schritt bzw. zwei halben Schritten).
  • Ist die Genauigkeit erreicht, gibt es das (extrapolierte) Ergebnis zurück.
  • Sonst teilt es das Intervall in zwei Hälften, ruft sich selbst zur Berechnung der Teilstücke auf und addiert die beiden Teilergebnisse.
• Als Matlab-Code:

  • function I = simpson_adapt(f, a, b, tol)
    % berechnet das Integral mit adaptivem Simpsonverfahren
    
    I1 = (b-a)/6 * (f(a) + 4*f((a+b)/2) + f(b));
    I2 = (b-a)/12 * (f(a) + 4*f((3*a+b)/4) + 2*f((a+b)/2) ...
                      + 4*f((a+3*b)/4) + f(b));
    
    if  abs(I1 - I2) < tol
      I = (16*I2 - I1)/15;
    else
      c = (a+b)/2;
      Il = simpson_adapt(f, a, c, tol);
      Ir = simpson_adapt(f, c, b, tol);
      I = Il + Ir;
    end

• Achtung: Bei numerischen Problemen mit dem Integral (etwa einem Pol) kann die Funktion in eine Endlosschleife geraten!
• Mit einem solchen rekursiven Verfahren lassen sich auch die Endpunkte bequem sammeln sowie alle Zwischenergebnisse ausgeben.
• Das Matlab-Skript ex17.m zeigt die komplette Lösung.