________________________________________________________________________________
Inhalt....: Das Euler-Cauchy Verfahren zur Lösung von Differenzialgleichung
Kategorie.: Arbeitsblatt
Mathematik: Analysis, Numerik
MuPAD.....: 3.0.0
Datum.....: 2002-03-06
Autoren...: Agnes Havasi <havasiagnes@yahoo.de>
Autoren...: Kai Gehrs <acrowley@mupad.de>
Funktionen: diff, plotfunc2d, solve, simplify, subs, plot, plot::Line2d,
Funktionen: plot::Point2d, plot::VectorField2d, plot::Ode2d, PointSize
Funktionen: numeric::ode2vectorfield, ode
________________________________________________________________________________
Das Euler-Cauchy Verfahren zur Lösung von
Differenzialgleichungen
In diesem Notebook stellen wir das Euler-Cauchy-Polygonzug Verfahren zur approximativen
grafischen Lösung von Differenzialgleichungen vor. Der mathematische Inhalt geht stellenweise
recht entscheidend über das hinaus, was standardmässig in Schulen behandelt werden kann.
Die Motivation bedient sich Beispielen der Physik, die ebenfalls recht komplex sind.
Es ist jedoch nicht unbedingt notwendig, ein physikalisch bzw. mathematisch
fundiertes Wissen auf dem Gebiet der gewühnlichen Differenzialgleichungen zu besitzen,
um die wesentlichen mathematischen Ideen, die hier vorgestellt werden, nachvollziehen
zu können. Vielmehr sollte man das Notebook als Appetithappen oder als eventuelle
Anregung für Facharbeiten bzw. mathematisch technisch interessierte Schülerinnen und
Schüler bzw. Studentinnen und Studenten verstehen.
Motivation:
Differenzialgleichungen treten in vielen Bereichen der Natur und Technik auf.
So wird zum Beispiel die Bewegung einer an einer elastischen Feder
aufgehängten Kugel der Masse m beschrieben durch die Gleichung
wobei y(t) den Ort der Kugel zur Zeit t beschreibt, y''(t) ihre Beschleunigung
und k eine Konstante ist. Diese entsteht durch Gleichsetzen der entsprechenden
Ausdrücke für die beteiligten Kräfte, nämlich
für die durch die Bewegung bedingte Kraft und
für die sogenannte Rückstellkraft der Feder mit Federkonstante c.
Ein weiteres Beispiel ist das Wachstum einer Bakterienkultur.
Beschreibt man die Populationsgröße der Kultur zur Zeit t mit y(t),
so ist das Wachstum in natürlicher Weise gegeben durch y'(t).
Dieses Wachstum ist proportional zur Größe der Population, also
Wir betrachten und behandeln diese Gleichung nun für den speziellen Fall
p=1:
Diese Gleichung ist eine sogenannte Differenzialgleichung. In Differenzial-
gleichungen tauchen Funktionen und ihre Ableitungen als Unbekannte auf.
Alle Funktionen, die die obige Gleichung erfüllen, müssen gleich ihrer
ersten Ableitung sein. Wir kennen alle die e-Funktion aus dem Schul-
unterricht. Leitet man diese Funktion ab, so erhält man wieder die gleiche
Funktion:
y := t -> exp(t);
y'
Diesen Sachverhalt kann man sich leicht anhand der folgenden beiden
Witze merken:
Witze in der Mathematik
Zwei Funktionen beschimpfen sich.
Nach einem längeren Austausch von "Komplimenten",
schließlich die eine: "Ich differenziere und integriere dich,
bis du nicht mehr weißt, wer du eigentlich bist!"
Darauf antwortet die andere: "Ätsch, ich bin e hoch x !"
Im Raum der stetigen Funktionen findet ein Tanzball statt. Auf der
Tanzfläche tanzen Cosinus und Sinus auf und ab und die Polynome
bilden einen Ring. Nur die Exponentialfunktion steht den ganzen Abend
alleine herum. Aus Mitleid geht die Identität irgendwann zu ihr hin und
sagt: "Mensch, integrier dich doch einfach mal!" "Schon versucht!",
antwortet die Exponentialfunktion, "Das hat aber auch nichts geändert!"
Wir kennzeichnen im folgenden mit y1 immer die erste Ableitung der
Funktion y.
Ist die Funktion exp(x) die einzige Lösung der obigen Differenzialgleichung?
Die Antwort ist nein, denn für eine beliebige reelle Zahl c gilt:
c * y(t) = diff(c * y(t), t)
d.h. alle Funktionen der Form y(x) = c * exp(x) erfüllen die obige Differenzial-
gleichung.
Wie kann man die Bedingungen an eine Lösung nun so weit verschärfen,
dass die Lösung eindeutig wird?
Zunächst schauen wir uns einige der Kurven der Form c * exp(x) an:
y_c:= t -> c * exp(t)
plotfunc2d(y_c(t) $ c = -4..4, t = -1..1)
Beobachtung: Jede Kurve hat mit der y-Achse einen anderen Schnittpunkt.
Nun liegt es nahe, dass wir die Lösungskurve durch eine Bedingung der
Form y(0) = d für eine reelle Zahl d "eindeutig machen", denn zu vorgegebenen
d können wir aus y(x) = c * exp(x) den entsprechenden Wert für c berechnen.
solve(d = c * y(0), c)
d.h. wir müssen c = d wählen.
BEISPIEL: d = 3
solve(3 = c * y(0), c)
plotfunc2d(y_c|c=3, t = -1..1)
Jetzt stellt sich die Frage: Wie löst man kompliziertere Gleichungen, wo man
die Lösung nicht wie oben "erraten" kann?
Ein mögliches Verfahren ist das sogenannte Euler-Cauchy-Polygonzugverfahren.
HISTORISCHES:
Euler (1707 - 1783) war einer der größten Mathematiker aller Zeiten.
Er entwickelte die Grundlagen der modernen Zahlentheorie und Algebra,
der Topologie, der Wahrscheinlichkeitsrechnung und Kombinatorik, der
Integralrechnung, der Theorie der Diffenrentialgleichungen und der
Differenzialgeometrie, der Variationsrechnung, entdeckte den Zusammenhang
zwischen trigonometrischen Funktionen und Exponentialfunktionen. Leonard
Euler entwickelte die Hydrodynamik und Strömungslehre, schuf die Grundlagen
für die Theorie des Kreisels. Er war ein genialer Naturwissenschaftler, ein
hervorragender Lehrer und Mentor.
Cauchy, geboren 21. 08. 1789 Paris, gestorben 22. 05. 1857 Sceaux (bei Paris).
Cauchy, dessen Vater hohe administrative Ämter bekleidete, erhielt eine
ausgezeichnete Privatausbildung. Er wollte Ingenieur werden und war auch
einige Jahre nach Beendigung seines Studiums als solcher tätig. Im
Selbststudium eignete er sich Werke von LAGRANGE und P.S. LAPLACE an.
Ab 1813 war Cauchy wieder in Paris. Seine mathematische Karriere begann aber
schon 1811, als es ihm gelang ein Problem von LAGRANGE zu lösen. Für die
Lösung eines Problems der Hydromechanik bekam Cauchy 1816 den Preis der
Pariser Akademie.
Mit seinen Arbeiten zur Elastizitätstheorie gehörte er spätestens ab 1822 zu
den herausragendsten Mathematikern seiner Zeit. In den Jahren 1815/16 war
Cauchy Professor an der Ecole Polytechnique und wurde 1816 Mitglied der
Französischen Akademie. Viele mathematische Begriffe und Sätze sind mit
seinem Namen verbunde, zum Beispiel der Satz von Cauchy, Cauchy-Folge,
Cauchy-Riemannsche Differenzialgleichungen, Cauchysche Integralformel,
Cauchyscher Integralsatz, Cauchysches Konvergenzkriterium. Cauchy gehört
zu den produktivsten Mathematikern aller Zeiten.
Zurück zur Mathematik an sich:
Die Idee des Verfahrens ist am leichtesten zu verstehen, wenn man sich
ein Beispiel anschaut. Wir wollen im folgenden die Differenzialgleichung
Wir beginnen mit einem Anfangswert y(0) = 2. Dann können wir den
Wert y'(0) ausrechnen, denn y(t) = y(0) = 2 und t = 0 (Rechnung siehe unten) 
delete y, y1:
y1:= t-> t * y - y^2
y1(0) | y=2
Wenn die Steigung bekannt ist, so können wir die Gerade durch (0|2) mit
Steigung -4 berechnen.
delete b:
g:= x -> -4 * x + b
solve( 2 = g(0), b)
Die Gerade ist also
g:= g | b=2
Die Idee ist nun: Berechne den Funktionswert von g an der Stelle 1/4 (d.h. wir gehen
einen viertel Schritt nach rechts):
g(1/4)
Das bedeutet: y(1/4)= 1.
l1 := plot::Line2d([0, 2], [1/4, 1]):
p1 := plot::Point2d([0, 2], PointSize = 3*unit::mm):
p2 := plot::Point2d([1/4, 1], Color = RGB::Blue,
PointSize = 3*unit::mm):
plot(l1,p1,p2):
Setzen wir das Verfahren fort. Wir gehen von dem neuen Punkt (1/4|1)
aus wieder um 1/4 nach rechts und berechnen die entsprechende Steigung
und den Funktionswert wie oben:
y1(1/4) | y=1
Die Steigung ist jetzt -3/4, die Gerade läuft durch (1/4|1). Wir bestimmen
jetzt eine weitere Gerade h mit Steidung -3/4, die durch den Punkt (1/4|1)
verläuft.
h:= x -> -3/4 * x + c
solve(1 = h(1/4), c)
Die Gerade sieht also folgendermaßen aus:
h:= h | c=19/16
Berechnen wir jetzt den Funktionswert von h an der Stelle 1/2, d.h. wir gehen
noch einen viertel Schritt nach rechts:
h(1/2)
Also gilt näherungsweise: y(1/2) = 13/16
Zeichnen wir sie mit den bisherigen Ergebnissen zusammen in ein
Koordinatensystem:
p3:= plot::Point2d([1/2, 13/16], Color = RGB::Green,
PointSize = 3*unit::mm):
l2:= plot::Line2d([1/4, 1], [1/2, 13/16]):
plot(p1,p2,p3,l1,l2):
Wir setzen das Verfahren fort:
y1(1/2) | y=13/16
Die Steigung der Gerade ist -65/256, die Gerade läuft durch den Punkt
(1/2 | 13/16). Diese Gerade bezeichnen wir mit k:
delete e:
k:= x -> -65/256 * x + e
solve( 13/16 = k(1/2), e)
Die Geradengleichung is daher gegeben durch:
k:= k | e=481/512
Berechnen wir jetzt den Funktionswert von k an der Stelle 3/4, d.h. wir
gehen einen weiteren viertel Schritt nach rechts:
k(3/4)
d.h. y(3/4) = 767/1024.
Zeichnen wir auch diese zusammen mit den vorherigen Ergebnissen:
p4:= plot::Point2d([3/4, 767/1024], Color = RGB::Pink,
PointSize = 3*unit::mm):
l3:= plot::Line2d([1/2, 13/16], [3/4, 767/1024]):
plot(p1,p2,p3,p4,l1,l2,l3):
Wenn man dieses Verfahren immer weiter fortsetzt, so erhält man einen
Polygonzug, der bei hinreichend kleiner Schrittweite, eine Kurve annähert.
Dieser Polygonzug nähert dann die Lösung der Differenzialgleichung
zu der festen Anfangsbedingung y(0) = 2 an.
Damit wir das Verfahren nicht in allen Schritten per Hand durchführen
müssen, bietet MuPAD eine Möglichkeit, zunächst die Steigung von
Differenzialgleichungen in Form von Richtungsfeldern in der Ebene zu
visualisieren.
Sehen wir uns eine geometrische Veranschaulichung des Richtungsfelds
der Differenzialgleichung y'(t) = t * y(t) - y(t) ^ 2 mit MuPAD an:
f:= (t, y) -> t*y - y^2
V1 := plot::VectorField2d([1, f(t, y)], t = 0..3, y = 0..3,
Mesh = [10, 10],
Color = RGB::Black):
plot(V1):
Wie interpretiert man das Bild? Man wählt auf der y-Achse einen Anfangswert,
z.B. y(0) = 2. Nun folgt man von diesem Punkt aus den Pfeilen in der Ebene.
Dies liefert eine Näherung der Lösung der Differenzialgleichung zum gegebenen
Anfangswert. Im obigen Beispiel, wo wir mit den Geraden gerechnet haben,
haben wir die Pfeile mathematisch über die Steigung der berechneten Geraden
angenähert.
Wir zeichnen unsere Näherungslösung mit dem Richtungsfeld in ein
Koordinatensystem:
plot(p1,p2,p3,p4,l1,l2,l3,V1):
Offensichtlich folgen wir mit dem obigen Verfahren zumindest im
Anfangswert y(0) = 2 unmittelbar dem eingezeichneten Pfeil.
Natürlich können wir mit MuPAD auch eine Kurve zu einem gegebenen
Anfangswert mit in die Grafik einzeichnen. Man sieht sehr schön, dass
die von MuPAD berechnete Näherung genauer ist, denn unsere Näherung
weicht deutlich davon ab.
f2 := plot::Ode2d(
(t,Y) -> [f(t, Y[1])], [i/3 $ i=0..9], [2],
[(t, Y) -> [t, Y[1]], Style = Points, Color = RGB::Red],
[(t, Y) -> [t, Y[1]], Style = Splines, Color = RGB::Blue]):
f2 := plot::modify(f2, PointSize = 3*unit::mm):
plot(V1, f2, p1,p2,p3,p4,l1,l2,l3,
Scaling = Constrained, GridVisible,
XAxisTitle = "t", YAxisTitle = "y",
XTicksDistance = 1.0, YTicksDistance = 0.5
):
Unsere Näherung verbessert sich dann, wenn wir die Schrittweite von
1/4 auf kleinere Werte 1/n für große n verkleinern.
Wie man sieht, interpoliert MuPAD zwischen den rot eingefärbten Punkten
NICHT mit linearen Funktionen, so wie wir es oben gemacht haben. Das
liegt daran, dass wir die Option Style = Splines verwendet haben.
Sie bewirkt, dass die Näherungsfunktion wirklich wie eine Kurve aussieht
und liefert daher offensichtlich viel genauere Werte als unser obiges
Verfahren.
Wir zeichnen noch einige Lösungskurven der betrachteten Differenzial-
gleichung zum Abschluß:
[f, t0, Y0] := [numeric::ode2vectorfield(
{y'(t) = t*y(t) - y(t)^2, y(0) = 2}, [y(t)])]:
G := (t, Y) -> [t, Y[1]]:
p1 := plot::Ode2d(f, [-0.5, i $ i = 0..4], Y0,
[G, Style = Points, Color = RGB::Blue],
[G, Style = Splines, Color = RGB::Red]):
Wählen wir noch 3 und 4 als weitere Anfangswerte.
[f, t0, Y0] := [numeric::ode2vectorfield(
{y'(t) = t*y(t) - y(t)^2, y(0) = 3}, [y(t)])]:
G := (t, Y) -> [t, Y[1]]:
p2 := plot::Ode2d(f, [-0.5, i $ i = 0..5], Y0,
[G, Style = Points, Color = RGB::Blue],
[G, Style = Splines, Color = RGB::Green]):
[f, t0, Y0] := [numeric::ode2vectorfield(
{y'(t) = t*y(t) - y(t)^2, y(0) = 4}, [y(t)])]:
G := (t, Y) -> [t, Y[1]]:
p3 := plot::Ode2d(f, [-0.5, i $ i = 0..5], Y0,
[G, Style = Points, Color = RGB::Blue],
[G, Style = Splines, Color = RGB::Black]):
plot(p1,p2,p3):
Wir betrachten noch einige weitere Beispiele aus dem Bereich der Differenzial-
gleichungen:
Wir zeichnen das Richtungssfeld und eine Näherung zur Anfangsbedingung
y(0) = 1/2:
f:= (t, y) -> t*sin(t + y^2) - y
p1:= plot::VectorField2d([1, f(t, y)], t = 0..4, y = -1..1,
Mesh = [10, 5], Color = RGB::Black):
p2:= plot::Ode2d(
(t,Y) -> [f(t, Y[1])], [i/3 $ i=0..12], [0.5],
[(t, Y) -> [t, Y[1]], Style = Points, Color = RGB::Red],
[(t, Y) -> [t, Y[1]], Style = Splines, Color = RGB::Blue]):
p2:= plot::modify(p2, PointSize = 3*unit::mm):
plot(p1, p2, Scaling = Constrained,
XAxisTitle = "t", YAxisTitle = "y",
GridVisible = TRUE,
XTicksDistance = 1.0,
YTicksDistance = 0.5):
Ein letztes Beispiel betrachten wir die Differenzialgleichung
f := (t, y) -> t*1/y
p1 := plot::VectorField2d([1, f(t, y)], t = 0..3, y = 0.1..3,
Mesh = [15, 15], Color = RGB::Black):
p2 := plot::Ode2d(
(t,Y) -> [f(t, Y[1])], [i/3 $ i=0..9], [0.5],
[(t, Y) -> [t, Y[1]], Style = Points, Color = RGB::Red],
[(t, Y) -> [t, Y[1]], Style = Splines, Color = RGB::Blue]):
p2 := plot::modify(p2, PointSize = 3*unit::mm):
plot(p1, p2, Scaling = Constrained,
XAxisTitle = "t", YAxisTitle = "y",
GridVisible = TRUE,
XTicksDistance = 1.0,
YTicksDistance = 0.5):
Nun stellt sich die Frage: Ist die Differenzialgleichung, die wir betrachten,
auch exakt lösbar?
Eine geschlossene Lösung können wir mit MuPAD wie folgt berechnen:
eq := ode(y'(x) = t*1/(y(x)), y(x))
solve(eq)
________________________________________________________________________________
Anmerkungen:
1. Weitere Anregungen finden Sie unter: http://schule.mupad.de bzw. http://studium.mupad.de
_______________________________________________________________________________