Schreiben Sie eine Funktion Ih = quadratur(f,a,b,N,regel) (Matlab)

Question

Aufgabe:

Schreibe die Funktion in Matlab.

a) Schreiben Sie eine Funktion Ih = quadratur(f,a,b,N,regel), die die folgenden Argumente (in obiger Reihenfolge) erwartet: Eine Funktion f : [a, b] → R, reelle Zahlen a, b eine ganze Zahl N und einen String regel, welche den Wert ’rechteck’,
’trapez’ oder ’simpson’ annehmen kann. Je nachdem, welchen Wert regel hat, soll das Integral ∫ (a,b) f(x)dx mit der Rechtecksregel, der Trapezregel und der Simpsonregel approximiert werden. N gibt die Anzahl der (äquidistanten) Teilintervalle an.

b) Bestimmen Sie den exakten Wert des Integrals ( 0 und pi/2) cos(x)e^sin(x)dx. Schreiben Sie ein Skript quad_plot, welche folgende Aufgabe erfüllt: Berechnen Sie die Approximation des Integrals fur alle in (a) genannten Verfahren und mit N = 2, 4, 8, 16, 32, 64.

Sei h die Länge der jeweiligen Teilintervalle. Generieren Sie mit dem Matlab-Befehl loglog für jede Quadraturformel ein Schaubild, welches den Logarithmus des Fehlers als Funktion von log(h) aufträgt. Zeichnen Sie in den je selben Plot die Funktion h bei der Rechtecksregel, h^2 bei der Trapezregel und h^4 bei der Simpsonregel. Was beobachten Sie? Können Sie das Verhalten erklären?

Hinweise: Der Befehl loglog funktioniert wie der Befehl plot, jedoch wird log(f(log(x))) über log(x) aufgetragen. Mit figure(k) können Sie Matlab anweisen, den folgenden Plot in das k-te Schaubild zu zeichnen.

Answer 1

Antwort mit KI erzeugt

Teil A: Die quadratur Funktion in Matlab

Um die Funktion quadratur in Matlab zu implementieren, folgen Sie diesen Schritten:

1. Die Funktion quadratur nimmt als Input die Funktion f, die untere und obere Grenze des Integrals a und b, die Anzahl der Teilintervalle N und den Namen der Regel regel an.
2. Abhängig von dem Wert von regel, wird die entsprechende numerische Integration-Methode gewählt: Rechteckregel, Trapezregel oder Simpsonregel.
3. Die Funktion gibt die Approximation des Integrals Ih zurück.

matlab
function Ih = quadratur(f, a, b, N, regel)
    % Bestimme die Länge eines Teilintervalls
    h = (b - a) / N;
    
    % Initialisiere x Werte der Teilintervall Knoten
    x = linspace(a, b, N+1);
    
    % Evaluieren der Funktion an jedem Knotenpunkt
    y = f(x);
    
    switch regel
        case 'rechteck'
            % Rechteckregel: Verwende den Wert am linken Rand des Teilintervalls
            Ih = h * sum(f(a + (0:N-1)*h));
            
        case 'trapez'
            % Trapezregel: Mittelwert aus linkem und rechtem Rand, dann summiert
            Ih = h * (sum(y(1:end-1)) + sum(y(2:end))) / 2;
            
        case 'simpson'
            % Simpsonregel: Spezielle Gewichtung der Knotenpunkte
            if mod(N,2) ~= 0
                error('N muss für die Simpsonregel gerade sein');
            end
            x_mid = a + (0.5:h:b-h);
            Ih = (h/6) * (y(1) + 2*sum(y(3:2:end-2)) + 4*sum(f(x_mid)) + y(end));
    end
end

Teil B: Exakter Wert des Integrals und Skript quad_plot

Um den exakten Wert des Integrals \(\int_{0}^{\frac{\pi}{2}} \cos(x) e^{\sin(x)} \, dx\) zu berechnen, benötigen wir eine analytische Lösung oder eine sehr genaue numerische Approximation. Wir können bemerken, dass \(\frac{d}{dx}e^{\sin(x)}=\cos(x)e^{\sin(x)}\), also ist das Integral von 0 bis \(\pi/2\) einfach die Differenz \(e^{\sin(x)}\) bewertet bei \(\pi/2\) und 0, also \(e^1 - e^0 = e - 1\).

Für das Skript quad_plot folgen Sie diesen Schritten:

1. Definieren und berechnen Sie die Approximation des Integrals für jede Methode und verschiedene Werte von N.
2. Verwenden Sie loglog, um den Fehler relativ zum exakten Wert grafisch darzustellen.
3. Zeichnen Sie zusätzlich die Referenzlinien für die erwartete Ordnung der Fehlerreduktion.

Hier ist ein grundlegendes Skript, das dies tut:

matlab
% Definiere die Funktion
f = @(x) cos(x).*exp(sin(x));
a = 0;
b = pi/2;
exactValue = exp(1) - 1;

Ns = [2, 4, 8, 16, 32, 64];
errors_rechteck = zeros(1, length(Ns));
errors_trapez = zeros(1, length(Ns));
errors_simpson = zeros(1, length(Ns));
hs = zeros(1, length(Ns));

for k = 1:length(Ns)
    N = Ns(k);
    hs(k) = (b-a)/N;
    
    % Berechnungen für jede Regel
    approx_rechteck = quadratur(f, a, b, N, 'rechteck');
    approx_trapez = quadratur(f, a, b, N, 'trapez');
    approx_simpson = quadratur(f, a, b, N, 'simpson');
    
    % Fehler speichern
    errors_rechteck(k) = abs(exactValue - approx_rechteck);
    errors_trapez(k) = abs(exactValue - approx_trapez);
    errors_simpson(k) = abs(exactValue - approx_simpson);
end

figure(1); % Für Rechteckregel
loglog(hs, errors_rechteck, 'r-o', hs, hs, 'b--');
legend('Fehler der Rechteckregel', 'h', 'Location', 'northwest');
title('Fehler der Rechteckregel');

figure(2); % Für Trapezregel
loglog(hs, errors_trapez, 'g-o', hs, hs.^2, 'b--');
legend('Fehler der Trapezregel', 'h^2', 'Location', 'northwest');
title('Fehler der Trapezregel');

figure(3); % Für Simpsonregel
loglog(hs, errors_simpson, 'b-o', hs, hs.^4, 'r--');
legend('Fehler der Simpsonregel', 'h^4', 'Location', 'northwest');
title('Fehler der Simpsonregel');

Was Sie beobachten werden, ist, dass die Fehlerkurven für jede Regel ungefähr der Rate folgen, die durch die hinzugefügten Referenzlinien vorgegeben wird: \(h\) für die Rechteckregel, \(h^2\) für die Trapezregel und \(h^4\) für die Simpsonregel. Dies bestätigt die theoretische Fehlerordnung für jede dieser numerischen Integrationsmethoden: Rechteckregel ist \(O(h)\), Trapezregel ist \(O(h^2)\), und Simpsonregel ist \(O(h^4)\), was bedeutet, dass, wenn die Länge des Teilintervalls halbiert wird, der Fehler für die Trapezregel um einen Faktor von vier reduziert wird und für die Simpsonregel sogar um einen Faktor von 16.