<<

. 9
( 11 .)



>>

end
input(™Druecke RETURN™)
hold off
clear
46 Matlab, ein mathematisches Labor




3. Lorenz-Gleichung
Die rechte Seite der DGL ist in ™lorenzf™ de¬niert.

function y = lorenzf(t,x)
y = zeros(3,1);
y(1) = 10.*(x(2)-x(1));
y(2) = -x(1).*x(3)+28.*x(1)-x(2);
y(3) = x(1).*x(2)-(8/3).*x(3);

[t,y] = ode45(™lorenzf™,[0 40], [0.1 0 0]);

Darstellung in Abh¨ngigkeit von t
a

plot(t,y(:,1),™k-™,t,y(:,2),™k:™,t,y(:,3),™k--™);
xlabel(™t™);
ylabel(™x_1(t) (-), x_2(t) (\cdot\cdot\cdot), x_3(t) (--)™);
input(™Druecke RETURN™)
clear
[t,y] = ode45(™lorenzf™,[0 100], [0.1 0 0]);

Darstellung als Kurve

figure
plot3(y(:,1),y(:,2),y(:,3),™k™);
xlabel(™x_1(t)™);
ylabel(™x_2(t)™);
zlabel(™x_3(t)™);
view(-110,-28);
input(™Druecke RETURN™)

Animierte Darstellung der Kurve

hold off
plot3(0, 0, 0)
axis([-20 20 -30 30 0 50])
view(-110,-28)
set(gca,™nextplot™,™replacechildren™)
s = size(t);
j = 1;
for i=1:400
while (t(j)<i/4)
j=j+1;
if (j==s(1))
Matlab, ein mathematisches Labor 47




break
end
end
plot3(y(1:j,1),y(1:j,2),y(1:j,3),y(j:j,1),y(j:j,2),y(j:j,3),™r.™)
F2(i) = getframe;
end
hold off

Nachfolgend illustrieren wir MATLABs Di¬erentialgleichungsl¨ser.
o
MATLAB stellt eine ganze Familie von L¨sern zur Verfugung:
o ¨
ode45 ode23 ode113 ode15s ode23s ode23tb ode23t
Wir werden hier die Loser ode45, ode23, ode15s und ode23s und ihre Anwen-
¨
dungsgebiete genauer untersuchen.
Als Standardl¨ser bietet sich fur allgemeine gew¨hnliche Di¬erentialgleichun-
o o
¨
gen die Routine ode45 an.

[t,x] = ode45(™dgl_f1™, [0 4], 1);
plot(t,x);

Dies ist ein eingebettetes Runge-Kutta Verfahren nach Dormand-Prince mit Kon-
sistenzordnung 4 und 5. Fur den Aufruf muss die rechte Seite der Di¬erentialglei-
¨
chung angegeben werden sowie das Losungs-Zeitintervall und der Anfangswert.
¨
Eine Schrittweite muss nicht angegeben werden, da diese bei allen Verfahren au-
tomatisch gesteurt wird.
¨
Auch wenn wir uns hier mit den L¨sungsroutinen besch¨ftigen wollen, soll kurz
o a
die gra¬sche Darstellung von Di¬erentialgleichungen als Richtungsfeld erwahnt¨
werden, die in MATLAB mit ™quiver™ erzeugt wird.

hold on
clear
[t,x] = meshgrid(0:0.2:4,-0.6:0.1:1);
z = dgl_f1(t,x);
quiver(t,x,ones(size(t)),z,™r-™);
axis([0 4 -0.6 1]);
hold off
input(™Druecke RETURN™)

Den DGL-L¨sern k¨nnen eine Reihe von Optionen ubergeben werden, die durch
o o ¨
die Verbundvariable ™options™ gesteurt werden, welche wiederum durch die An-
¨
weisung ™odeset™ gesetzt wird. Einige der wichtigen Optionen sind
RelTol: Gewunschte Relative Genauigkeit Voreinstellung 1e-3 (1 Prozent)
¨
AbsTol: Gewunschte Absolute Genauigkeit Voreinstellung 1e-6
¨
48 Matlab, ein mathematisches Labor




Die Gesamtgenauigkeit wird aus diesen Werten mittels der Formel: Fehler
¡ max(RelTol*abs(x(i)),AbsTol) berechnet, wobei ™Fehler™ den in jedem Schritt
hinzukommenden lokalen Fehler bezeichnet. Ist die Bedingung verletzt, wird uber
¨
eine Schrittweitensteuerung die Schrittweite h verkleinert (Details dazu in der
Numerik-Vorlesung am 15.7.)
Re¬ne: Feinheit der Ausgabe, steuert die Anzahl der Punkte, die ausgegeben
wird. Voreinstellung: 4 fur ™ode45™ 1 fur alle anderen L¨ser
o
¨ ¨
Stats: Ausgabe des numerischen Aufwands eines Verfahrens ™on™ oder ™o¬™
(Voreinstellung)
Vectorized: Angabe, ob das M-File der rechten Seite der DGL Vektoreingaben
korrekt verarbeiten kann ™on™ oder ™o¬™ (Voreinstellung)
MaxStep: Maximal erlaubte Schrittweite Voreinstellung: 1/10 der Rechenin-
tervallbreite
InitialStep: Vorschlag fur Anfangsschrittweite Voreinstellung: Automatische
¨
Auswahl
NormControl: ™on™: Die Genauigkeit der L¨sung wird uber die 2-Norm gemes-
o ¨
sen ™o¬™: Die Genauigkeit wird komponentenweise gemessen (Voreinstellung)
Im Folgenden werden wir eine Reihe von Optionen und ihre Auswirkungen
illustrieren. Hier verwenden wir neben ode45 auch ode23, was ein ¨hnliches Ver-
a
fahren mit (niedrigerer) Konsistenzordnung 2 und 3 ist.

clear
options = odeset(™Stats™, ™on™, ™AbsTol™, 1e-10);
fprintf(™\node45, Absolute Genauigkeit 1e-10\n™);
[t,x] = ode45(™dgl_f1™, [0 4], 1, options);
fprintf(™\node23, Absolute Genauigkeit 1e-10\n™);
[t,x] = ode23(™dgl_f1™, [0 4], 1, options);
options = odeset(™Stats™, ™on™, ™AbsTol™, 1e-4);
fprintf(™\node45, Absolute Genauigkeit 1e-4\n™);
[t,x] = ode45(™dgl_f1™, [0 4], 1, options);
fprintf(™\node23, Absolute Genauigkeit 1e-4\n™);
[t,x] = ode23(™dgl_f1™, [0 4], 1, options);
input(™Druecke RETURN™)

Bei dieser Gleichung ist ode45 o¬enbar klar e¬zienter, was an der hoheren Kon-
¨
sistenzordnung liegt. Es kann aber sinnvoll sein, Verfahren niedrigerer Konsi-
stenzordnung zu verwenden, wenn die rechte Seite der Di¬erentialgleichung nicht
hinreichend oft di¬erenzierbar ist.
Der ™re¬ne™ Parameter ist insbesondere fur eine sch¨nere gra¬sche Ausgabe
o
¨
nutzlich
¨
Matlab, ein mathematisches Labor 49




clear
[t,x] = ode45(™dgl_f1™, [0 4], 1);
plot(t,x);
input(™Druecke RETURN™)
figure
clear
options = odeset(™Refine™, 20);
[t,x] = ode45(™dgl_f1™, [0 4], 1, options);
plot(t,x);
input(™Druecke RETURN™)
Mit re¬ne=1 erh¨lt man genau die ™Rechenpunkte™ (ti , xi ), was auch fur die gra-
a ¨
¬sche Ausgabe sinnvoll verwendet werden kann.
hold on
clear
options = odeset(™Refine™, 1);
[t,x] = ode45(™dgl_f1™, [0 4], 1, options);
plot(t,x,™rx™);
hold off
input(™Druecke RETURN™)
Neben den ™options™ k¨nnen noch weitere Parameter an die DGL-L¨ser ubergeben
o o¨
werden. Diese bein¬‚ussen allerdings nicht die numerische Routine, sondern wer-
den direkt an die rechte Seite der Di¬erentialgleichung weitergegeben. Auf diese
Weise kann man Di¬erentialgleichungen fur verschiedene Parameter l¨sen, ohne
o
¨
das M-File der DGL zu andern.
¨
Wir veranschaulichen das an einer 3d Di¬erentialgleichung, die als ™R¨ssler-
o
System™ bekannt ist.
clear
options = []; % Ein Trick, weil hier eine ™options™-Variable
% uebergeben werden muss, auch wenn keine
% Optionen geaendert werden sollen
[t,x]=ode45(™dgl_f2™,[0,100],[1; 1; 1], options, 0.2, 0.2, 2.5);
plot3(x(:,1),x(:,2),x(:,3))
title(™c=2.5™)
grid
figure
[t,x]=ode45(™dgl_f2™,[0,100],[1; 1; 1], options, 0.2, 0.2, 5);
plot3(x(:,1),x(:,2),x(:,3))
title(™c=5™)
50 Matlab, ein mathematisches Labor




grid
input(™Druecke RETURN™)

Die bisher betrachteten Verfahren sind sogenannte ™explizite™ Einschrittverfahren,
bei denen die Losung x(i+1) direkt aus der vorhergehenden Losung x(i) berechnet
¨ ¨
wird. Fur gewisse Klassen von Di¬erentialgleichungen, die sogenannten ™steifen™
¨
Di¬erentialgleichungen sind diese Verfahren sehr ine¬zient, auch wenn sie nicht
unbedingt ein falsches Ergebnis liefern. Steife Di¬erentialgleichungen zeichnen
sich dadurch aus, dass die Konstante ™K(T)™ im Konvergenzsatz sehr gross ist,
was sich typischerweise durch eine grosse Lipschitz-Konstante fur die rechte Seite
¨
der Di¬erentialgleichung f ausdruckt. Geometrisch betrachtet haben wir es mit
¨
L¨sungen zu tun, die sich zun¨chst sehr schnell ver¨ndern k¨nnen, sich dann aber
o a a o
lange Zeit fast konstant verhalten.
Fur diese Gleichungen sind die sogenannten impliziten Verfahren deutlich bes-
¨
ser geeignet. Hierbei erh¨lt man keine explizite Formel fur x(i+1), sondern - et-
a ¨
was vereinfacht gesagt - ein nichtlineares Gleichungssystem der Form x(i+1) =
F(x(i),x(i+1)), dass in jedem Schritt (z.B. mit dem Newton-Verfahren) gel¨st o
werden muss.
MATLAB stellt mit ode15s und ode23s zwei solcher Verfahren zur Verfugung¨
(ode23s ist ein Einschrittverfahren, ode15s ein sogenantes Mehrschrittverfah-
ren, bei dem der Wert x(i+1) aus mehreren vorhergehenden Werten x(i-s), x(i-
s+1),...,x(i) berechnet wird). Wir vergleichen die Ergebnisse zun¨chst an einer
a
einfachen 1d DGL und dann an einer komplizierteren 3d DGL.

clear
options = odeset(™Stats™, ™on™, ™Refine™, 1);
tt = [0 0.5];
x0 = 1;
fprintf(™\node45:\n™)
[t1,x1] = ode45(™dgl_f3™, tt, x0, options);
plot(t1,x1,™rx™,t1,x1);
title(™ode45™)
axis([0 0.5 -1 1])
figure

fprintf(™\node15s:\n™)
[t2,x2] = ode15s(™dgl_f3™, tt, x0, options);
plot(t2,x2,™rx™,t2,x2);
title(™ode15s™)
axis([0 0.5 -1 1])

<<

. 9
( 11 .)



>>