6. Numerische Differentiation und Integration

6. Numerische Differentiation und Integration

In dieser Übung geht es um das numerische Berechnen der Ableitung bzw. des Integrals. Dabei wird auch auf vergangene Übungen zurückgegriffen. So müssen zum Beispiel Messwertlücken zunächst interpoliert werden, bevor eine Änderungsrate berechnet werden kann.

Situation

Die numerische Berechnung von Ableitungen und Integralen spielt für die Berechnung und Simulation technischer Systeme eine entscheidende Rolle. Genauso entscheidend ist es aber auch, Einschränkungen und Grenzen der eingesetzten Hard- und Software zu kennen.

Daher ist das Thema dieser Übung numerische Differentiation und Quadratur.

Differentiation 1

Berechnen Sie nach dem Ansatz in der Vorlesung (auch bekannt als „h-Methode“) numerisch die Ableitungen der u.g. Funktionen an den genannten Stellen.

Dazu soll eine Matlab-Function angelegt werden, welche als Eingang ein beliebiges Function-Handle sowie den Ort \(x0\) erhält und als Ausgang den Wert der Ableitung bei \(x0\) zurückgibt.

Um den Einfluss der Schrittweite h zu prüfen, sollen auch dafür verschiedene Werte verwendet werden.

\[x_0=0,1,10,100,1000\]
\[h=10^{-4},10^{-8},10^{-12},10^{-16}\]

Berechnen Sie zur Bewertung der Ergebnisse die analytischen Ableitungen und stellen Sie die Werte gegenüber.

\[f_1 (x)=e^{x^2}\]
\[f_2 (x)=\ln⁡(x)\]
\[f_3 (x)= \sin⁡(x)\]

Tipp

Nutzen Sie zum Abarbeiten der verschiedenen Werte eine Schleife.

% your code here

Differentiation 2

Bei der Durchführung verschiedener Messreihen sind jeweils folgende Werte ermittelt worden:

t

-4

6

5

1

3

y

1

1

5

8

4

t

0

1.5708

3.1416

4.7124

6.2832

y

0

1

0

-1

0

t

0

2.5

5.0

7.5

10

y

1

1

1

0

0

Für eine simple Änderungsauswertung kann zunächst der Befehl diff getestet werden. Er gibt den Wert einer Änderung von \(y_i\) zu \(y_{i+1}\) aus.

Da die Funktionswerte hier aber nur an bestimmten Stellen bekannt sind, muss ein Interpolationspolynom zur Approximation einer Ableitung verwendet werden.

Schreiben Sie eine function, welche als Eingang \(x\)- und \(y\)-Wertepaare sowie einen Abfrageort \(x0\) enthält. Die function soll dann ein Interpolationspolynom bestimmen und dessen Ableitung an der Stelle \(x0\) zurückgeben.

Tipp

Die Koeffizienten des abgeleiteten Polynomes, z.B. mit den Koeffizienten \(p = [a, b, c]\), können mittels polyder(p) bestimmt werden.

Schauen Sie auch nochmals in die Unterlagen in 3. Polynomfitting.

% your code here
function fdot = derivate(f,x0)
end

Quadratur

Berechnen Sie die Fläche welche zwischen den beiden Funktionen \(f_1\) und \(f_2\) im Intervall \([0,3]\) aufgespannt wird.

\[f_1(x) = e^{\frac{x}{2}}\]
\[f_2(x) = 4\sin(x)\]

Tipp

Ermitteln Sie zunächst die Schnittpunkte mit einem Newtonverfahren (z.B. fsolve) und verwenden Sie dann die Routine integral.

% your code here

Lösung

- Hinweis: Dieses Video ist von YouTube aus eingebunden und nicht Teil des frei lizenzierten Materials! -

Hier finden Sie außerdem den Code der Lösung.

Differentiation 1:

%% Hauptprogramm
function main
    % Definition der Schrittweite h als globale Variable
    global h
    % Verschiedene Testpunkte x0
    x0 = [0;1;10;100;1000];
    % Berechnung der Ableitung mit verschiedenen Schrittweiten h
    for i = [4,8,12,16]
        h = 10^(-i);
        % Ausgabe von jeweils numerischer und analytischer Ableitung
        [ableit(@f1,x0),df1(x0)]
    end
end

%% Numerische Berechnung der Ableitung
function dy = ableit(f,x0)
    global h
    dy = (f(x0+h)-f(x0))/h;
end

%% Testfunktion: Exponentialfunktion
% Funktion f
function y = f1(x)
    y = exp(x.^2);
end
% Ableitung df
function dy = df1(x)
    dy = 2*x.*exp(x.^2);
end

Differentiation 2:

%% Hauptprogramm
% Tipp: Schauen Sie sich hierzu auch nochmals Übung 3 (Polynomfitting) an.
function main
    %-- Messreihe t/y (siehe Aufgabenblatt)
    t = [-4,6,5,1,3];
    y = [1,-1,5,8,4];
    %-- Erzeugung äquidistanter Zeitwerte zwischen tmin und tmax
    xp = linspace(min(t),max(t),100);
    %-- Gesuchter Ort (kann frei gewählt werden)
    x0 = 2;
    %-- Plot der Messreihe
    hold off
    plot(t,y,'rx','LineWidth',2)
    grid on
    hold on
    %-- Berechnung Interpolationspolynomkoeffizienten 4. Grades
    %  (weil 5 Messwerte verfügbar sind)
    p = polyfit(t,y,4);
    yp = polyval(p,xp);
    %-- Plot des Polynoms
    plot(xp,yp)
    %-- Berechnung der Ableitungskoeffizienten des Polynoms
    dp = polyder(p);
    dy = polyval(dp,xp);
    %-- Plot der Polynomableitung
    plot(xp,dy)
    %-- Berechnung des Ableitungswertes bei x0 mit eigener Funktion (s.u.)
    dy = ableit2(t,y,x0)
end

%% Berechnung der Ableitung des Interpolationspolynoms
function dy = ableit2(x,y,x0)
    % Interpolationspolynom
    N = length(x);
    p = polyfit(x,y,N-1);
    % Ableitungskoeffizienten
    dp = polyder(p);
    % Wert der Ableitung bei x0
    dy = polyval(dp,x0);
end

Quadratur:

%% Hauptprogramm
function main
    %-- Definition eines Plot-Intervalls (sodass eingeschlossene Fläche
    %   sichtbar ist)
    xp = linspace(0,3,100);
    %-- Plot von f1 und f2
    hold off
    plot(xp,f1(xp));
    hold on
    plot(xp,f2(xp));
    %-- Berechnung der Schnittpunkte zwischen f1 und f2 mittels
    %   Newtonverfahren fsolve.
    %   Die Startwerte 0.2 und 0.3 sind durch Ablesen am obigen Plot
    %   gewählt worden.
    xa = fsolve(@df,0.3);
    xb = fsolve(@df,2.2);
    %-- Berechnung des Integrals zwischen den gerade berechneten Grenzen xa
    %   und xb mithilfe der Routine "integral".
    A = integral(@df,xa,xb)
end

%% Differenzfunktion f1-f2
function y = df(x)
    y = f2(x)-f1(x);
end

%% Funktion f1
function y = f1(x)
    y = exp(x./2);
end

%% Funktion f2
function y = f2(x)
    y = 4*sin(x);
end