<<

. 8
( 11 .)



>>

P7 = polyfit(xx,ff,2)
hold on;
y = polyval(P7,t);
plot(t,y,™g-™)
hold off;
input(™Druecke RETURN™)
Matlab, ein mathematisches Labor 41




MATLAB hat auch Routinen zur Erzeugung und Auswertung von Funktionen,
die stuckweise aus Polynomen bestehen. Wir werden in der Vorlesung noch sehen,
¨
dass diese zur L¨sung von Interpolationsproblemen sehr nutzlich sein k¨nnen.
o o
¨
Nehmen wir an, dass wir eine Funktion de¬nieren wollen, die durch die fol-
genden Formel beschrieben wird:

f (x) = (x + 10)2 /81 ’ 2 falls x >= ’10 und x <= ’1
f (x) = (x + 1) ’ 1 falls x >= ’1 und x <= 1
f (x) = (x ’ 1)3 + (x ’ 1)2 + 1 falls x >= 1 und x <= 10
Die MATLAB Konvention sieht hierbei vor, dass jede Teilfunktion ein Poly-
nom in der Variablen (x-xi) ist, wobei xi die linke Intervallgrenze des betre¬enden
De¬nitionsabschnittes ist.
Um diese Funktion in MATLAB zu realisieren, de¬nieren wir zunachst einen
¨
Vektor mit den ™Nahtstellen™ (inklusive der Randpunkte):
xx = [-10 -1 1 10]
Dann de¬nieren wir die drei Polynome, die die Funktion auf den entsprechenden
Abschnitten de¬nieren. Beachte: Im Allgemeinen ist es egal, ob die Polynom- Ko-
e¬zienten in Zeilen oder Spaltenvektoren gespeichert werden. Fur diese Anwen-
¨
dung mussen es aber zwingend Zeilenvektoren sein. Ausserdem mussen alle Vek-
¨ ¨
toren die gleiche Dimension haben, weswegen wir PP2 mit einer Null ™au¬ullen™
¨
mussen.
¨
PP1 = [1/81 0 -2]
PP2 = [0 1 -1]
PP3 = [1 1 1]
Mit der Anweisung ™mkpp™ erzeugen wir nun die gewunschte Funktion
¨
PP = mkpp(xx,[PP1; PP2; PP3])
Beachte: PP ist eine sogenannte ™Verbundvariable™ (engl. ™struct™), in der alle
n¨tigen Werte gespeichert sind.
o
Mit ™ppval™ k¨nnen wir diese Funktion nun an beliebigen Stellen auswerten.
o
ppval(PP,0)
ppval(PP,-5)
und auch plotten
t = [-3:0.01:3];
y = ppval(PP,t);
plot(t,y,™-™)
42 Matlab, ein mathematisches Labor




2.5 Integration
In diesem Abschnitt wollen wir MATLABs Integrationsroutinen kurz besprechen.
ACHTUNG: Diese Version ist fur MATLAB 6, fur MATLAB 5.3 gibt es leichte
¨ ¨
Ab¨nderungen.
a
Wir wollen das Integral
b
f (x) dx
a
numerisch berechnen. Als Beispielintegranden betrachten wir die Funktion

f (x) = 1/((x ’ 0.3)2 + 0.01) + 1/((x ’ 0.9)2 + 0.04) ’ 6,

die in MATLAB unter dam Namen ™humps™ vorde¬niert ist.

fplot(™humps™,[0,1])

Das Integral von 0 bis 1 uber diese Funktion ist etwa 29.858325395
¨
MATLAB hat zwei Newton-Cotes Formeln eingebaut, die aber unterschiedlich
aufgerufen werden.
Zun¨chst gibt es die Trapez-Regel, die mittels ™trapz™ zur Verfugung steht. Hier
a ¨
wird allerdings keine Funktion sondern einfach eine Liste von Datenpunkten f(i)
zu ¨quidistanten Stutzstellen ubergeben. Um das richtige Integrationsergebnis zu
a ¨ ¨
erhalten, muss das Ergebnis durch die Anzahl der Stutzstellen -1 geteilt werden
¨

n = 100;
t = [0:1/n:1];
f = humps(t);
trapz(f)/n
¨
Wegen der expliziten Ubergabe von Daten eigent sich diese Routine insbeson-
dere zur Integration von Funktionen, die nur als Messwerte oder Wertetabellen
vorliegen.
Zur Integration von ™echten™ Funktionen dienen die Routinen ™quad™ und ™quadl™.
™quad™ benutzt die Simpson-Regel w¨hrend ™quadl™ eine sogenannte ™(4,7) Gauss-
a
Lobatto-Kronrod Regel™ erwendet. Diese ist bei hinreichend oft di¬erenzierbaren
Funktionen i.A. e¬zienter.
Der Aufruf dieser Funktionen lautet ™quad(™fun™,a,b,tol,trace)™, mit:
fun : zu integrierende Funktion
a,b : Integrationsgrenzen
tol : gewunschte Genauigkeit (Voreinstellung 1e-6)
¨
trace: =0: Keine Ausgabe von Zwischenergebnissen (Voreinstellung), =1: Ausga-
be von Zwischenergebnissen
Matlab, ein mathematisches Labor 43




tol und trace konnen weggelassen werden, dann wird die Voreinstellung gewahlt.
¨ ¨
Bei ™trace=1™ werden als Zwischenergebnisse jeweils die Anzahl der Funktions-
auswertungen, der Anfangspunkt und die L¨nge der verwendeten Teilintervalle
a
und der Integralwert auf diesem Teilinterval ausgegeben.

quad(™humps™,0,1)
quad(™humps™,0,1,1e-15)
quad(™humps™,0,1,1e-6,1)
input(™Druecke RETURN™)

Beide Routinen verwenden eine sogenannte adaptive Stutzstellenwahl, d.h. die
¨
Stutzstellen werden abh¨ngig von der vorgegebenen Genauigkeit iterativ aus-
a
¨
gew¨hlt und sind im Allgemeinen nicht ¨quidistant.
a a
Mit dem Aufruf

[I,p] = quad(™humps™,0,1)
input(™Druecke RETURN™)

wird in der Variablen ™p™ die Anzahl der Funktionsauswertungen zuruckgegeben.
¨
Leider liefert ™quad™ keine explizite Liste der verwendeten Stutzstellen. Mit
¨
einem Trick k¨nnen wir diese trotzdem gra¬sch darstellen: Wir de¬nieren eine
o
Funktion ™myhumps™

function y = myhumps(x)
y = humps(x);
plot(x,y,™ro™)
plot(x,zeros(size(x)),™k+™)

die bei jedem Aufruf den Wert von ™humps™ zuruckgibt und zugleich einen kleinen
¨
Kreis auf dem Graphen und ein ™+™ auf der x-Achse an der entsprechenden Stelle
zeichnet.

hold on
quad(™myhumps™,0,1,1e-4)
hold off
input(™Druecke RETURN™)

MATLAB kann auch Doppelintegrale fur Funktionen mit zwei Variablen berech-
¨
nen, also
b1 b2

f (x, y) dydx
a1 a2

Die Anweisung hierfur lautet ™dblquad(™fun™,a1,b1,a2,b2,tol,method)™ mit
¨
44 Matlab, ein mathematisches Labor




fun : Zu integrierende Funktion in x und y
a1,b1 : Integrationsgrenzen in x
a2,b2 : Integrationsgrenzen in y
tol : gewunschte Genauigkeit (Voreinstellung 1e-6)
¨
method: @quad (Voreinstellung) oder @quadl
Als Methode kann auch der Name einer selbstgeschriebenen Integrationsroutine
mit fuhrendem ™@™ eingegeben werden, sofern diese den Parameterkonventionen
¨
der eingabuten Routinen entspricht.

dblquad(™sin(x).*y™,0,pi,0,1)


2.6 Di¬erentialgleichungen
Wir betrachten einige Beispieldi¬erentialgleichungen
1. Radioaktiver Zerfall. Hier kennen wir die exakte L¨sung, die wir einfach
o
plotten konnen
¨

t = 0:0.1:10;
plot(t,exp(-t.*0.5),™k-™);
xlabel(™t™);
ylabel(™x_1(t)™);
input(™Druecke RETURN™)
clear

2. Restringiertes Drei-K¨rper-Problem
o
Hierbei mussen wir die voreingestellte Genauigkeit etwas erh¨hen, um eine
o
¨
sch¨ne L¨sung zu erhalten
o o
Die rechte Seite der DGL ist in ™satellitf™ de¬niert.

function y = satellitf(t,x)
y = zeros(4,1);
mu1 = 0.012277471;
mu2 = 1 - mu1;
D1 = ((x(1) + mu1).^2 + x(3).^2).^(3/2);
D2 = ((x(1) - mu2).^2 + x(3).^2).^(3/2);
y(1) = x(2);
y(2) = x(1) + 2.*x(4) - mu2.*(x(1) + mu1)./ D1 ...
- mu1.*(x(1) - mu2)./ D2;
y(3) = x(4);
y(4) = x(3) - 2.*x(2) - mu2.*x(3) ./ D1 - mu1.*x(3)./ D2;
Matlab, ein mathematisches Labor 45




options = odeset(™RelTol™,1e-10);
[t,y] = ode45(™satellitf™,[0,17.066], [0.994,0,0,-2.0015851063790825],...
options);

Darstellung in Abh¨ngigkeit von t.
a

plot(t,y(:,1),™k-™,t,y(:,3),™k:™);
xlabel(™t™);
ylabel(™x_1(t) (-), x_3(t) (\cdot\cdot\cdot)™);
input(™Druecke RETURN™)

Darstellung als Kurve

figure
plot(y(:,1),y(:,3),™k™);
axis([-1.5,1.5,-1.5,1.5]);
xlabel(™x_1(t)™);
ylabel(™x_3(t)™);
input(™Druecke RETURN™)

Animierte Darstellung der Kurve

hold off
plot(-0.0123, 0, ™ko™, 1-0.0123, 0, ™k.™)
axis([-1.5,1.5,-1.5,1.5]);
set(gca,™nextplot™,™replacechildren™)
s = size(t);
j = 1;
for i=1:69
while (t(j)<i/4)
j=j+1;
if (j==s(1))
break
end
end
plot(0, 0, ™ko™, 1, 0, ™k.™,y(1:j,1),y(1:j,3))
F1(i) = getframe;

<<

. 8
( 11 .)



>>