Skip to content
Permalink
Browse files

Serie 4

  • Loading branch information
bauerda3 committed Mar 29, 2020
1 parent 553e4a4 commit 2964b024f78da657a2e274ebd2b21f65c0bf33a6
Binary file not shown.
@@ -0,0 +1,77 @@
% NMIT2
% Serie 4
% Aufgabe 1
% Daniel Bauer (bauerda3@students.zhaw.ch)
% 27.03.2020, IT18aWin, FS20
clear, clc, close all

t = 0:10:120;
h = [2, 286, 1268, 3009, 5375, 8220, 11505, 15407, 20127, 25593, ...
31672, 38257, 44931];
m = [2051113, 1935155, 1799290, 1681120, 1567611, 1475282, 1376301, ...
1277921, 1177704, 1075683, 991872, 913254, 880377];

R0 = 6378137;
M = 5.976e24;
G = 6.673e-11;

% Teilaufgabe a
v = ableitung_wertepaare(t, h);
a = ableitung_wertepaare(t, v);
dm_dt = ableitung_wertepaare(t, m);
dm_dh = ableitung_wertepaare(h, m);

figure('Name', 'Space Shuttle', 'NumberTitle', 'off')

subplot(4, 2, 1), plot(t, h, 'm')
title('Hoehe'), xlabel('t'), ylabel('h'), grid minor

subplot(4, 2, 2), plot(t, v, 'c')
title('Geschwindigkeit'), xlabel('t'), ylabel('v'), grid minor

subplot(4, 2, 3), plot(t, a, 'r')
title('Beschleunigung'), xlabel('t'), ylabel('a'), grid minor

subplot(4, 2, 4), plot(t, m, 'g')
title('Masse'), xlabel('t'), ylabel('m'), grid minor

subplot(4, 2, 5), plot(t, dm_dt, 'b')
title('Veraenderung Masse zur Zeit'), xlabel('t'), ylabel('dm/dt'), grid minor

subplot(4, 2, 6), plot(h, dm_dh, 'k')
title('Veraenderung Masse zur Hoehe'), xlabel('h'), ylabel('dm/dt'), grid minor

subplot(4, 2, [7, 8]), plot(t, dm_dt, t, dm_dh .* v)
title('Veraenderung Masse zur Zeit = Veraenderung Masse zur Hoehe * Geschwindigkeit'), xlabel('h'), ylabel('dm/dt'), grid minor

% Teilaufgabe b
R0_h = h + R0;
Ekin = kinetische_energie(R0_h, dm_dt, v, m, a, t);
Epot = potentielle_energie(R0_h, G, M, m);
E = Ekin + Epot;

figure('Name', 'Energie', 'NumberTitle', 'off')
plot(t, E, t, Ekin, t, Epot)
title('Energieverlauf'), xlabel('t'), ylabel('E'), grid minor
legend('E', 'Ekin', 'Epot')

% Teilaufgabe c
E_Haushalt = 10^10;
N_Haushalt = E(end) / E_Haushalt;
fprintf("Anzahl Haushalte die 1 Jahr lang versorgt werden koennten: %.1f\n", N_Haushalt)

function Ekin = kinetische_energie(R0_h, dm_dt, v, m, a, t)
Ekin = zeros(size(t));
for i = 1 : length(Ekin)
Ekin1 = summierte_trapezregel_nicht_aequidistant(R0_h(1:i), dm_dt(1:i) .* v(1:i));
Ekin2 = summierte_trapezregel_nicht_aequidistant(R0_h(1:i), m(1:i) .* a(1:i));
Ekin(i) = Ekin1 + Ekin2;
end
end

function Epot = potentielle_energie(R0_h, G, M, m)
Epot = zeros(size(m));
for i = 1 : length(Epot)
Epot(i) = summierte_trapezregel_nicht_aequidistant(R0_h(1:i), G * M * m(1:i) ./ R0_h(1:i).^2);
end
end
@@ -0,0 +1,40 @@
% NMIT2
% Serie 4
% Aufgabe 2
% Daniel Bauer (bauerda3@students.zhaw.ch)
% 27.03.2020, IT18aWin, FS20
clear, clc, close all

V = @(t) 3500 * sin(140 * pi * t) * exp(-63 *pi * t);
E = @(t) V(t)^2 / 50;

% Teilaufgabe a
t = 0 : 0.01 : 0.05;
E3 = zeros(size(t));
E6 = zeros(size(t));
for i = 1 : length(t)
[E3i, ~] = romberg_extrapolation(E, 0, t(i), 3);
[E6i, ~] = romberg_extrapolation(E, 0, t(i), 6);
E3(i) = E3i(1, end);
E6(i) = E6i(1, end);
end
figure('Name', 'Defibrillator', 'NumberTitle', 'off')
plot(t, E3, t, E6), title('Puls eines Defibrillators')
legend('Romberg m = 3', 'Romberg m = 6'), xlabel('t'), ylabel('E'), grid on

% Teilaufgabe b
f = @(t) E(t) - 250;
i = 1;
t(1) = 0.005;
tol = 10^-5;
while f(t(i) - tol) * f(t(i) + tol) < 0
t(i + 1) = t(i) - (f(t(i)) / f_diff(f, t(i)));
i = i + 1;
end
fprintf("fDiff: df = %e\n", f_diff(f, 0.005))
fprintf("Newton-Verfahren: t = %.3f\n", t(end))

function y = f_diff(f, x)
[D, ~] = h2_extrapolation(f, @vorwaertsdifferenz, x, 0.01, 3);
y = D(1, end);
end
@@ -0,0 +1,80 @@

%----------------------------------------------------------------------
% Anleitung zu Aufgabe 1 in Serie 4
%----------------------------------------------------------------------

% Gegebene Daten

% t = 0:10:120;
% h = [2, 286, 1268, 3009, 5375, 8220, 11505, 15407, 20127, 25593, ...
% 31672, 38257, 44931];
% m = [2051113, 1935155, 1799290, 1681120, 1567611, 1475282, 1376301, ...
% 1277921, 1177704, 1075683, 991872, 913254, 880377];
%
% R0 = 6378137;
% M = 5.976e24;
% G = 6.673e-11;

% Teilaufgabe a)

% In a) sind aus den gegebenen Listen t, h und m entsprechende Listen fuer
% die Geschwindigkeit v = dh/dt, die Beschleunigung a = dv/dt sowie
% die Massensenderungen dm_dt = dm/dt und dm_dh = dm/dh zu berechnen.
% Dazu koennen Sie die in Aufgabe 3 a) von Serie 1 erstellte Funktion
% verwenden. Fuer zum Beispiel v ergibt der Aufruf
%
% v = Name_Vorname_Klasse_S1_Aufg3a(t,h)
%
% das Ergebnis
%
% v = [28.4000, 63.3000, 136.1500, 205.3500, 260.5500, 306.5000, ...
% 359.3500, 431.1000, 09.3000, 577.2500, 633.2000, 662.9500, 667.4000]

% Teilaufgabe b)

% Zur Berechnung der Integrale in b) muss wegen den Integrationsgrenzen
% noch eine Liste erstellt mit den um den Erdradius erhoehten Hoehen:
%
% R0_h = h + R0;
%
% Fuer die Integrale koennen Sie die in Aufgabe 4 a) von Serie 3 erstellte
% Funktion verwenden. Dazu zwei Beispiele:
%
% - Um Ekin(t) zum Endzeitpunkt t = 120 auszurechnen, ergeben die Aufrufe
%
% Ekin1 = Name_Vorname_Klasse_S3_Aufg4a(R0_h, dm_dt.*v);
% Ekin2 = Name_Vorname_Klasse_S3_Aufg4a(R0_h, m.*a);
% Ekin = Ekin1 + Ekin2
%
% das Ergebnis
%
% Ekin = 1.0593e+11
%
% - Um Ekin(t) zu einem Zwischenzeitpunkt t zu berechnen, verwenden Sie nur
% den Teil der Daten bis zu t. Fuer Ekin(t) zum Zeitpunkt t = 30 geht das
% wie folgt:
%
% Ekin1 = Name_Vorname_Klasse_S3_Aufg4a(R0_h(1:4), dm_dt(1:4).*v(1:4));
% Ekin2 = Name_Vorname_Klasse_S3_Aufg4a(R0_h(1:4), m(1:4).*a(1:4));
% Ekin = Ekin1 + Ekin2
%
% Das Ergenis lautet
%
% Ekin = 2.9141e+10
%
% UM Epot(t) ZU BERECHNEN, IST ZU BEACHTEN, DASS IM INTEGRAND NICHT WIE
% AUF DEM AUFGABENBLATT ANGEGEBEN DURCH h^2 SONDERN DURCH (R0 + h)^2 zu
% ZU TEILEN IST.
%
% Um Epot(t) zum Endzeitpunkt t = 120 zu berechnen, lautet der Aufruf
%
% Epot = Name_Vorname_Klasse_S3_Aufg4a(R0_h, G*M*m./(R0_h).^2)
%
% und ergibt
%
% Epot = 5.2572e+11

% Teilaufgabe c)

% Kein Hinweis noetig

@@ -0,0 +1,66 @@

%----------------------------------------------------------------------
% Anleitung zu Aufgabe 2 in Serie 4
%----------------------------------------------------------------------

% Gegebene Daten

V = @(t) 3500*sin(140*pi*t)*exp(-63*pi*t);
R = 50;

% Teilaufgabe a)

% DIE ANGABE t AUS [0,50] IST UNVOLSTAENDIG. GEMEINT IST t AUS
% [0 MS, 50 MS]. FUER DIE INTEGRATION IST t DEMENTSPRECHEND AUS
% [0, 0.05] ZU WAEHLEN.
%
% In Aufgabe 3 von Serie 3 haben Sie ein Funktion erstellt zur
% des Integrals von a bis b der Funktion f(x) m-stufiger Rhomberg-
% Extrapolation. Fuer diev Integration koennen Sie darauf zurueck-
% greifen. Um E(t) mit dieser Funktion zum Zeitpunkt t = 0.01
% mit 3-stufiger Extrapolation zu berechnen, lautet der Aufruf
%
% E = Name_Vorname_Klasse_S3_Aufg3(@(t) V(t)^2/R, 0, 0.01, 3)
%
% Er ergibt
%
% E = 248.3496
%
% Der Grund fuer das eigenartige Verhalten von E(t) ist subtil und
% wohl nicht einfach zu ermitteln. Schuld ist die Anfangsschrittweite h0
% der Extrapolation. In der Funktion ist sie mit h0 = b - a implementiert.
% Es ist als hier h0 = t. Mit zunehmendem t wird h0 gross und ergibt
% schlechte Ergebnisse.
%
% Man kann die Genauigkeit auch fuer grosse t sicherstellen, indem man
% auf mehr Stufen m extrapoliert. Ausreichend ist hier m = 6.

% Teilaufgabe b)

% Berechnen Sie E(t) auf die genaue Art mit 6-stufiger Extrapolation.
% Sie koennen E(t) wie folgt kompakt definieren:
%
% E = @(t) Name_Vorname_Klasse_S3_Aufg3(@(t) V(t)^2/R, 0, t, 6);
%
% Der Beispielaufruf E(0.02) ergibt 257.2870.
%
% Fuer die Newtoniteration ist die Funktion f(t) = E(t) - 250 zu verwenden.
%
% Fuer die Bestimmung der in der Newton-Iteration auftrenden
% Ableitung f'(t) verwenden Sie die in Aufgabe 2 von Serie 2 erstellte
% Funktion. Diese bestimmt die Ableitung einer Funktion f an einer Stelle
% x0 durch m-stifige h^2-Extrapolation mit Startschrittweite h0. Waehlen
% Sie h0 = 0.01 und m = 3.
%
% Wenn Sie die Iteration mit t0 = 0.005 starten ist im ersten Schritt
% die Ableitung f'(0.005) zu bilden. Das geht via
%
% df = Name_Vorname_Klasse_S2_Aufg2(f, 0.005, 0.01, 3)
%
% und ergibt
%
% df = 2.2451e+04




No changes.
@@ -0,0 +1,31 @@
% Ableitung von Wertepaaren mit D1f, D2f, D3f

% Beispiel: Serie 1, Aufgabe 3a
% n = 1000;
% i = 1:n;
% t = 0 + i * 0.1;
% x = 10 * exp(-0.05 * t) .* cos(0.2 * pi * t);
% dx = ableitung_wertepaare(t, x)

function dx = ableitung_wertepaare(x, y)
if length(x) ~= length(y)
error 'Vektoren x und y haben nicht gleiche Länge';
end

n = length(x);
dx = zeros(1, n);

% Vorwärtsdifferenz (D1f)
h1 = x(2) - x(1);
dx(1) = (y(1 + 1) - y(1)) / h1;

for i = 2 : n - 1
% Zentrale Differenz (D2f)
h2 = x(i + 1) - x(i - 1);
dx(i) = (y(i + 1) - y(i - 1)) / h2;
end

% Rückwärtsdifferenz (D3f)
h3 = x(n) - x(n - 1);
dx(n) = (y(n) - y(n - 1)) / h3;
end
@@ -0,0 +1,29 @@
% h2-Extrapolation mit zugehörigen Diskretisierungsfehlern
% Ergebnis der Extrapolation in D(1,m+1)

% Beispiel: Kapitel 6, Aufgabe 6.5
% [D, E] = h2_extrapolation(@(x) sin(x), @zentrale_differenz, 1, 0.1, 3)

function [D, E] = h2_extrapolation(f, Dxf, x0, h, m)
D = zeros(m+1);
E = zeros(m+1);

% Ableitung benötigt für Diskretisierungsfehler
fDiff = matlabFunction(diff(sym(f)));

% Erste Spalte mit D2f füllen
for row = 1:m+1
i = row - 1;
D(row,1) = Dxf(f, x0, h/2^i);
E(row,1) = abs(D(row,1) - fDiff(x0));
end

% Rekursion
for col = 2:m+1
k = col - 1;
for row = 1:m+1-k
D(row,col) = (4^k * D(row+1, col-1) - D(row, col-1)) / (4^k - 1);
E(row,col) = abs(D(row,col) - fDiff(x0));
end
end
end
@@ -0,0 +1,32 @@
% Romberg-Extrapolation mit zugehörigen Fehlern
% Ergebnis der Extrapolation in T(1,m+1)

% Beispiel: Kapitel 6, Beispiel 6.4
% [T, E] = romberg_extrapolation(@(x) 1/x, 2, 4, 3)

% Beispiel: Kapitel 6, Aufgabe 6.10
% [T] = romberg_extrapolation(@(x) exp(-x^2), 0, 0.5, 3)

function [T, E] = romberg_extrapolation(f, a, b, m)
T = zeros(m+1);
E = zeros(m+1);

% Integral benötigt für Fehler
fInt = matlabFunction(int(sym(f)));

% Erste Spalte mit summierter Trapezregel füllen
for row = 1:m+1
i = row - 1;
T(row,1) = summierte_trapezregel(f, a, b, 2^i);
E(row,1) = abs(T(row,1) - (fInt(b) - fInt(a)));
end

% Rekursion
for col = 2:m+1
k = col - 1;
for row = 1:m+1-k
T(row,col) = (4^k * T(row+1, col-1) - T(row, col-1)) / (4^k - 1);
E(row,col) = abs(T(row,col) - (fInt(b) - fInt(a)));
end
end
end
@@ -0,0 +1,20 @@
% Summierte Trapezregel fuer nicht aequidistante x-Werte

function Tf_neq = summierte_trapezregel_nicht_aequidistant(x, y)
if size(x) ~= size(y)
error("x und y muessen dieselbe Dimension haben");
end

Tf_neq = 0;

if (size(x, 2) == 1)
return
end

y_diff = (y(1:end-1) + y(2:end)) ./ 2;
x_diff = x(2:end) - x(1:end-1);

flaeche_intervalle = y_diff .* x_diff;
Tf_neq = sum(flaeche_intervalle);

end

0 comments on commit 2964b02

Please sign in to comment.
You can’t perform that action at this time.