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.
Berechnen Sie zur Bewertung der Ergebnisse die analytischen Ableitungen und stellen Sie die Werte gegenüber.
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.
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#
Hier finden Sie außerdem den Code der Lösung.
Differentiation 1:
Show code cell source
%% 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:
Show code cell source
%% 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:
Show code cell source
%% 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