Einführung in die Numerik Mathematik, FS 2009 Prof. D. Cohen / Dr. T. Mitkova, M. Utzinger Universität Basel
Musterlösungen zur Serie 11
Aufgabe 1: (P)*
Die M-Datei NewtonIter_System.m:
function [x, iter] = NewtonIter_System(f1, f2, x0, max_iter, tol) h = 1e-8; % Schrittweite zum Differenzenquotienten
x = [x0(1); x0(2)]; % Startwert iter = 0; % Anzahl der Iterationen
while (max(abs(f(x))) >= tol & iter < max_iter)
% Jacobi-Matrix
Df(1, 1) = (f1(x(1)+h, x(2)) - f1(x(1), x(2)))/h;
Df(1, 2) = (f1(x(1), x(2)+h) - f1(x(1), x(2)))/h;
Df(2, 1) = (f2(x(1)+h, x(2)) - f2(x(1), x(2)))/h;
Df(2, 2) = (f2(x(1), x(2)+h) - f2(x(1), x(2)))/h;
x = x - Df\[f1(x(1), x(2)); f2(x(1), x(2))]; % die Loesung iter = iter + 1;
end
Die M-Datei Test_NewtonIter_System˙m:
f1 = inline(’sinh(x1*x2) + x1^2 + x2^2 + x1 - 1;’, ’x1’, ’x2’);
f2 = inline(’x1^3 - x2^2 + x2 + 1’, ’x1’, ’x2’);
x0 = [-1, 1];
max_iter = 20;
tol = 10^(-12);
[loesung, iter] = NewtonIter_System(f1, f2, x0, max_iter, tol)
Mit dem Startwert x(0) = (−1,1)T erhalten wir nach 7Schritten die Lösung mit der gewünschten Genauigkeit: x≃ −0.606290595782463,y≃1.513476530934236.
Aufgabe 2: (P)*
Die M-Datei Euler.m:
function [x, y] = Euler(f, y0, h, x_0, x_max) x = x_0 : h : x_max;
y = zeros(length(x), 1);
y(1) = y0;
for i = 1 : (length(x)-1)
y(i+1) = y(i) + h*f(x(i), y(i));
end
Die M-Datei Test_Euler.m
f = inline(’(x+2)*y/(x+1)’,’x’,’y’);
loesung = inline(’(1+x).*exp(1+x)’, ’x’);
[x, y] = Euler(f, exp(1), 0.01 , 0, 5);
plot(x, loesung(x), x, y)
legend(’Exakte Loesung’, ’Numer. Loesung’, ’Location’,’NorthWest’) axis([0 5.5 -0.001 2500])
Aufgabe 3: (P)
Die M-Datei Euler_System.m:
function [x, y1, y2] = Euler_System(f1, f2, y0, h, x_0, x_max) x = x_0 : h : x_max;
y1 = zeros(length(x), 1);
y2 = zeros(length(x), 1);
y1(1) = y0(1);
y2(1) = y0(2);
for i = 1 : (length(x)-1)
y1(i+1) = y1(i) + h*f1(x(i), y1(i), y2(i));
y2(i+1) = y2(i) + h*f2(x(i), y1(i), y2(i));
end
Die M-Datei Mod_Euler_System.m:
function [x, y1, y2] = ModEuler_System(f1, f2, y0, h, x_0, x_max) x = x_0 : h : x_max;
y1 = zeros(length(x), 1);
y2 = zeros(length(x), 1);
y1(1) = y0(1);
y2(1) = y0(2);
for i = 1 : (length(x)-1)
y1(i+1) = y1(i) + h*f1(x(i)+h/2, ...
y1(i)+(h/2)*f1(x(i), y1(i), y2(i)), ...
y2(i)+(h/2)*f1(x(i), y1(i), y2(i)));
y2(i+1) = y2(i) + h*f2(x(i)+h/2, ...
y1(i)+(h/2)*f2(x(i), y1(i), y2(i)), ...
y2(i)+(h/2)*f2(x(i), y1(i), y2(i)));
end
Die M-Datei Test_Euler_System.m:
y0 = [0, 1];
h = 0.025;
x_0 = 0; x_max = 5;
% Test mit lambda = 0
f1 = inline(’y2’,’x’,’y1’, ’y2’);
f2 = inline(’-y1’,’x’,’y1’, ’y2’);
loesung_y1 = inline(’sin(x)’, ’x’);
loesung_y2 = inline(’cos(x)’, ’x’);
[x, y1, y2] = Euler_System(f1, f2, y0, h, x_0, x_max);
figure(1)
subplot(1, 2, 1)
plot(x, y1, x, loesung_y1(x))
legend(’numerisch’, ’exakt’, ’Location’, ’SouthWest’) xlabel(’x’); ylabel(’y_1(x)’)
axis([-0.5 5.5 -1.2 1.2]) subplot(1, 2, 2)
plot(x, y2, x, loesung_y2(x)) legend(’numerisch’, ’exakt’) xlabel(’x’); ylabel(’y_2(x)’) axis([-0.5 5.5 -1.2 1.2])
% Test mit lambda = 12
f1 = inline(’y2’,’x’,’y1’, ’y2’);
f2 = inline(’12*(1-y1^2)*y2-y1’,’x’,’y1’, ’y2’);
[x, y1, y2] = Euler_System(f1, f2, y0, h, x_0, x_max);
figure(2)
subplot(1, 2, 1) plot(x, y1)
xlabel(’x’); ylabel(’y_1(x)’) axis([-0.5 5.5 -0.2 2])
subplot(1, 2, 2) plot(x, y2)
xlabel(’x’); ylabel(’y_2(x)’) axis([-0.5 5.5 -0.2 11])
Die M-DateiTest_ModEuler_System.mist ähnlich zur DateiTest_Euler_System.m.
Aufgabe 4: (A+B)
Existenz: Seix0 ∈Ωbeliebig. Wir betrachten die Fixpunktiterationxn+1 =Φ(xn), n= 0,1,2, . . .. DaΦ eine Kontraktion ist, gilt
kxn+1−xnk=kΦ(xn)−Φ(xn−1)k ≤Lkxn−xn−1k ≤ · · · ≤Lnkx1−x0k, für alle n ∈ N und ein L <1. Dabei ist die k · k eine beliebige Norm auf Rn. Damit zeigen wir, dass {xn}n≥0 eine Cauchy-Folge ist. Für m > n(O. B. d. A.) gilt
kxm−xnk=
m−n
X
k=1
(xn+k−xn+k−1)
≤
m−n
X
k=1
kxn+k−xn+k−1k
≤ kx1−x0k
m−n
X
k=1
Ln+k−1 ≤ kx1−x0kLn
∞
X
k=0
Lk = Ln
1−Lkx1−x0k. Da Ln → 0, n → ∞, existiert also tatsächlich zu jedem ǫ > 0 ein N ∈ N, so dass kxm−xnk< ǫfür alle n, m≥N, und damit ist {xn}n≥0 eine Cauchy-Folge.
DaΩ⊂Rn abgeschlossen und beschränkt, also kompakt, ist, konvergiert die Cauchy- Folge gegen einen Grenzwert x∗ ∈Ω. Wegen der Stetigkeit von Φ gilt
x∗ = lim
n→∞xn+1 = lim
n→∞Φ(xn) =Φ( lim
n→∞xn) =Φ(x∗). Also ist x∗ ein Fixpunkt von Φ.
Eindeutigkeit: Sei x′ ∈Ωein weiterer Fixpunkt von Φ. Dann gilt kx∗−x′k=kΦ(x∗)−Φ(x′)k ≤Lkx∗−x′k.
Damit gilt kx∗ −x′k = 0, da L < 1 und 0 ≤ (1−L)kx∗−x′k ≤ 0. Da k · k eine Norm ist, folgt x∗ =x′.
Aufgabe 5: (A+B)
Da Φ stetig differenzierbar ist, ist Φ′ stetig auf Ω. Also existiert ein ε > 0 so, dass kΦ′(x)k<1, x∈Bε(x∗). Deshalb gilt für x,y∈Bε:
kΦ(x)−Φ(y)k=
1
Z
0
Φ′(tx+ (1−t)y)dt(x−y)
≤ max
ξ∈{tx+(1−t)y:t∈[0,1]}kΦ′(ξ)k·kx−yk.
Da Bε konvex ist, gilt {tx+ (1−t)y:t ∈[0,1]} ⊂Bε und damit
ξ∈{tx+(1−t)y:t∈[0,1]}max kΦ′(ξ)k=:L <1.
Also gilt
kΦ(x)−Φ(y)k ≤Lkx−yk,
und Φ ist eine Kontraktion auf Bε. Um zu zeigen, dass Bε in sich abgebildet wird, nehmen wir ein x∈Bε und rechnen
kΦ(x)−x∗k=kΦ(x)−Φ(x∗)k ≤Lkx−x∗k<kx−x∗k ≤ε .
Also istΦ(x)∈Bε. Dies gilt für jedes x∈Bε, also wirdBε vonΦin sich abgebildet.
Allgemeine Informationen zur Vorlesung und Übungsblätter befinden sich auf der Webseite http://www.math.unibas.ch/˜cohen
