Blog

H02 – Iterative Integralberechnung

Zu berechnen sind die Integrale {I_n} = \int_0^1 {\frac{{{x^n}}}{{x+20}}dx} ,\quad n \geq 0.

  1. Zeigen Sie \frac{1}{{21}} \cdot \frac{1}{{n+1}} \leq {I_n} \leq \frac{1}{{20}} \cdot \frac{1}{{n+1}}
  2. Zeigen Sie, dass die Integrale {I_n} der linearen Rekursion {I_n} = \frac{1}{n}-20{I_{n-1}} genügen
  3. Wir gehen im Folgenden von exakter Gleitpunktarithmetik aus. Der Startwert {I_0} = \log 21-\log 20 lässt sich im Rechner nicht exakt darstellen. Stattdessen wird ein gerundeter Wert {\tilde I_0} = {I_0}+\Delta {I_0} verwendet. Wie pflanzt sich der Fehler \Delta {I_0} fort? Welchen absoluten Fehler \Delta {I_n} hat der berechnete Wert {I_n}?
  4. Zu berechnen seien {I_1}, \ldots {I_{10}}. Wäre {I_{10+k}} für ein gewisses k bekannt, könnte man die Rückwärtsrekursion {I_{n-1}} = \frac{1}{{20}}\left( {\frac{1}{n}-{I_n}} \right) verwenden. Wie pflanzt sich der Startfehler \Delta {I_{10+k}} fort? Wie groß muss k sein, um die Integrale {I_{10}} bis {I_0} in Matlab mit voller Genauigkeit zu berechnen, wenn man als Startwert {\tilde I_{10+k}} = \frac{1}{{20,5}} \cdot \frac{1}{{11+k}} wählt?
  5. Schreiben Sie zwei Matlab-Funktionen zur Berechnung der Vorwärts- und Rückwärtsrekursion. Berechnen Sie die Werte für {I_1} bis {I_{10}} mit beiden Funktionen und geben Sie das Ergebnis in Form einer Matrix mit den drei Spalten n, vorwaerts, rueckwaerts aus. Vergleichen Sie anschließend mit den exakten Werten.

Lösung

a )

Da in dem Integral x \in \left[ {0,1} \right] ist, schätzen wir nach unten bzw. oben ab mit x = 1 bzw. x = 0.

{I_n} = \int_0^1 {\frac{{{x^n}}}{{x+20}}dx} \leq \int_0^1 {\frac{{{x^n}}}{{0+20}}dx = \frac{1}{{20}}\int_0^1 {{x^n}dx} } = \frac{1}{{20}} \cdot \left[ {\frac{{{x^n}}}{{n+1}}} \right]_0^1 = \frac{1}{{20}} \cdot \frac{1}{{n+1}}

{I_n} = \int_0^1 {\frac{{{x^n}}}{{x+20}}dx} \geq \int_0^1 {\frac{{{x^n}}}{{1+20}}dx = \frac{1}{{21}}\int_0^1 {{x^n}dx} } = \frac{1}{{21}} \cdot \left[ {\frac{{{x^n}}}{{n+1}}} \right]_0^1 = \frac{1}{{21}} \cdot \frac{1}{{n+1}}

\quad \Rightarrow \quad \frac{1}{{21}} \cdot \frac{1}{{n+1}} \leq {I_n} \leq \frac{1}{{20}} \cdot \frac{1}{{n+1}}

b )

{I_n} = \frac{1}{n}-20{I_{n-1}},\quad {I_n} = \int_0^1 {\frac{{{x^n}}}{{x+20}}dx} ,\quad n \geq 0

\quad \Rightarrow \quad \int_0^1 {\frac{{{x^n}}}{{x+20}}dx} = \frac{1}{n}-20\int_0^1 {\frac{{{x^{n-1}}}}{{x+20}}dx}

\quad \Rightarrow \quad \int_0^1 {\frac{{{x^n}}}{{x+20}}dx} +\int_0^1 {\frac{{20{x^{n-1}}}}{{x+20}}dx} = \frac{1}{n}

\quad \Rightarrow \quad \int_0^1 {\frac{{{x^n}+20{x^{n-1}}}}{{x+20}}dx} = \frac{1}{n}

\quad \Rightarrow \quad \left[ {\frac{1}{n} \cdot {x^n}} \right]_0^1 = \frac{1}{n}

\quad \Rightarrow \quad \frac{1}{n} = \frac{1}{n}

c )

{{\tilde I}_0} = {I_0}+\Delta {I_0}\quad \Rightarrow \quad \Delta {I_0} = {{\tilde I}_0}-{I_0}

{{\tilde I}_1} = 1-20{{\tilde I}_0} = 1-20\left( {{I_0}+\Delta {I_0}} \right) = \underbrace {1-20{I_0}}_{ = {I_1}}-20\Delta {I_0} = {I_1}-20{I_0}

\Rightarrow \quad \Delta {I_1} = -20{I_0}

{{\tilde I}_2} = \frac{1}{2}-20{I_1} = \frac{1}{2}-20\left( {{I_1}+\Delta {I_1}} \right) = \underbrace {\frac{1}{2}-20{I_1}}_{ = {I_2}}-20\left( {-20{I_0}} \right)

\Rightarrow \quad \Delta {I_2} = {\left( {-20} \right)^2}{I_0}

Allgemein:

\Delta {I_n} = {\left( {-20} \right)^n}\Delta {I_0}

Beweis mit vollständiger Induktion:

\Delta {I_n} = {\left( {-20} \right)^n}\Delta {I_0},\quad n \to n+1

\Delta {I_{n+1}} = {{\tilde I}_{n+1}}-{I_{n+1}} = \frac{1}{{n+1}}-20\left( {{I_n}+\Delta {I_n}} \right)-{I_{n+1}}

\Delta {I_{n+1}} = \underbrace {\frac{1}{{n+1}}-20{I_n}}_{ = {I_{n+1}}}-20\Delta {I_n}-{I_{n+1}} = -20\Delta {I_n}

Einsetzen der Induktionsannahme:

\Delta {I_n} = {\left( {-20} \right)^n}\Delta {I_0},\quad n \to n+1

\Delta {I_{n+1}} = {{\tilde I}_{n+1}}-{I_{n+1}} = \frac{1}{{n+1}}-20\left( {{I_n}+\Delta {I_n}} \right)-{I_{n+1}}

\Delta {I_{n+1}} = -20{\left( {-20} \right)^n}\Delta {I_0} = {\left( {-20} \right)^{n+1}}\Delta {I_0}\quad \quad \quad \quad \square

d )

{I_{n-1}} = \frac{1}{{20}}\left( {\frac{1}{n}-{I_n}} \right)

Der Fehler bei der Rückwärtsrekursion ist:

\Delta {I_{n-1}} = {\tilde I_{n-1}}-{I_{n-1}} = \frac{1}{{20}}\left( {\frac{1}{n}-{{\tilde I}_n}} \right)-{I_{n-1}} = \underbrace {\frac{1}{{20}}\left( {\frac{1}{n}-{I_n}} \right)}_{ = {I_{n-1}}}-\frac{1}{{20}}\Delta {I_n}-{I_{n-1}} = -\frac{1}{{20}}\Delta {I_n}

Mit Induktion folgt:

\Delta {I_0} = {\left( {-\frac{1}{{20}}} \right)^n}\Delta {I_n}\quad \Rightarrow \quadstarke Dämpfung des Fehlers

Es sollen die ersten 10 Werte berechnet werden, wobei der Fehler bei {I_{10}} noch am größten ist:

\Delta {I_{10}} = {\left( {-\frac{1}{{20}}} \right)^{k-10}}\Delta {I_k} = {\left( {-\frac{1}{{20}}} \right)^k}\Delta {I_{k+10}}

Dabei ist der Fehler des Startwertes höchstens halb so groß wie das Intervall, in dem der Startwert liegt. Wir schätzen daher nach oben ab:

\frac{1}{{21}} \cdot \frac{1}{{n+1}} \leq {I_n} \leq \frac{1}{{20}} \cdot \frac{1}{{n+1}}\quad \Rightarrow \quad \Delta {I_n} = \frac{1}{2}\left( {\frac{1}{{20}}-\frac{1}{{21}}} \right)\frac{1}{{n+1}} = \frac{1}{{840}} \cdot \frac{1}{{n+1}}

Dieser Fehler darf für den Wert {I_{10}} nicht zu groß sein:

\left| {\Delta {I_{10}}} \right| = {\left( {\frac{1}{{20}}} \right)^k}\Delta {I_{k+10}} = {\left( {\frac{1}{{20}}} \right)^k} \cdot \frac{1}{{840}} \cdot \frac{1}{{k+11}}

Dieser Wert soll nun kleiner als {10^{-16}} sein. Es ergibt sich:

\left| {\Delta {I_{10}}} \right| \leq {10^{-16}}\quad \Rightarrow \quad {\left( {\frac{1}{{20}}} \right)^k} \cdot \frac{1}{{840}} \cdot \frac{1}{{k+11}} \leq {10^{-16}}\quad \Rightarrow \quad k \geq 10

e )

vorwaerts.m

function int = vorwaerts(I0, n)
% int: vector I1, ..., I10
% I0: start value I0
% n: index of last integral

int = zeros(1, n);  % pre-allocating for speed
int(1) = 1 - 20 * I0;

for i = 2 : n
    int(i) = (1 / i) - 20 * int(i - 1);
end

rueckwaerts.m

function int = rueckwaerts(Ink, n, k)
% int: vector I1, ..., Ik
% Ink: start value I(n + k)
% n  : index of last integral
% k  : number of precalculated integrals

intAll = zeros(1, n + k);
intAll(n + k) = Ink;

for i = n + k - 1 : -1 : 1
    intAll(i) = 1 / 20 * (1 / (i + 1) - intAll(i + 1));
end

int = zeros(1, n);

for i = 1 : n
    int(i) = intAll(i);
end

Aufruf in u14.m

format long

% exact values:
I = zeros(1, 10);
I(1) = 0.024196716611359938;
I(2) = 0.016065667772801226;
I(3) = 0.012019977877308810;
I(4) = 0.0096004424538237932;
I(5) = 0.0079911509235241352;
I(6) = 0.0068436481961839618;
I(7) = 0.0059841789334636197;
I(8) = 0.0053164213307276052;
I(9) = 0.0047826844965590059;
I(10) = 0.0043463100688198801;

% calculate forward integrals:
n = 10;
intFw1 = vorwaerts(log(21) - log(20), n);
intFw2 = vorwaerts(log(21/20), n);

% calculate backward integrals:
k = 10;
startValue = 1 / (20.5 * (11 + k)); 
intBw1 = rueckwaerts(startValue, n, k);
k = 0;
startValue = 1 / (20.5 * (11 + k)); 
intBw2 = rueckwaerts(startValue, n, k);

% calculating forward errors:
devFw1 = I - intFw1;
devFw2 = I - intFw2;

% calculating backward errors:
devBw1 = I - intBw1;
devBw2 = I - intBw2;

% gathering results:
result1 = zeros(n, 3);
result1(:, 1) = 1 : n;
result1(:, 2) = intFw2;
result1(:, 3) = intBw2;
result1

Ausgabe:

n fw log(20)-log(21) deviation fw log(21/20) deviation
1 0.024196716611357 0.000000000000003 0.024196716611359 0.000000000000001
2 0.016065667772864 -0.000000000000063 0.016065667772819 -0.000000000000018
3 0.012019977876056 0.000000000001253 0.012019977876944 0.000000000000365
4 0.009600442478890 -0.000000000025066 0.009600442461126 -0.000000000007302
5 0.007991150422208 0.000000000501316 0.007991150777480 0.000000000146044
6 0.006843658222499 -0.000000010026315 0.006843651117072 -0.000000002920888
7 0.005983978407162 0.000000200526301 0.005984120515710 0.000000058417754
8 0.005320431856752 -0.000004010526024 0.005317589685809 -0.000001168355081
9 0.004702473976077 0.000080210520482 0.004759317394938 0.000023367101621
10 0.005950520478451 -0.001604210409632 0.004813652101235 -0.000467342032415
n bw k = 10 dev bw k = 0 deviation
1 0.024196716611360 0 0.024196716611360 0.000000000000000
2 0.016065667772801 0 0.016065667772805 -0.000000000000003
3 0.012019977877309 0 0.012019977877240 0.000000000000069
4 0.009600442453824 0 0.009600442455203 -0.000000000001379
5 0.007991150923524 0 0.007991150895937 0.000000000027587
6 0.006843648196184 0 0.006843648747932 -0.000000000551748
7 0.005984178933464 0 0.005984167898497 0.000000011034966
8 0.005316421330728 0 0.005316642030057 -0.000000220699329
9 0.004782684496559 0 0.004778270509978 0.000004413986581
10 0.004346310068820 0 0.004434589800443 -0.000088279731624

Wie erwartet wird der Fehler bei der Vorwärtsrekursion größer, bei der Rückwärtsrekursion liegt er wie berechnet für k=10 unterhalb der Maschinengenauigkeit. Selbst für k = 0 ergeben sich bessere Werte als mit der Vorwärtsrekursion.

Leave a Reply

Close Menu