Differentialgleichungen Demonstrationsübung I

Visualisierungen zum expliziten Euler-Verfahren / Taylor / Runge-Kutta
Mathematica und Differentialgleichungen

Numerische Analysis & Computer-Algebra mit Mathematica	
Prof. Dr. Bernhard Zgraggen, HSR

Mini-Einführung Mathematica

Dieses Dokument ist ein "selbsterklärendes" Mathematica-Notebook, das interaktiv bearbeitet werden kann ("Hands-On"). Input-Zellen können mittels der Tastenkombination ⇧⌤ (in Worten: Shift-Enter) ausgeführt werden. Falls Textzellen geschlossen sind, so können sie durch Klick auf das Icon ▼ oder auf die blaue Zellklammer rechterhand geöffnet werden.

Eine allgemeine Hands-On-Einführung in Mathematica mit Übungen ist im englischen Notebook 01Introduction.nb vorhanden. 

Wenn Sie in Mathematica auf Befehle oder Begriffe stossen, die Sie nicht verstehen, so steht Ihnen die mächtige Hilfe/Doc Center zur Verfügung; diese ist interaktiv.

Allfällige Aufgaben/Übungen befinden sich am Ende des Dokuments.

Mathematica und Differentialgleichungen -> DSolve und NDSolve

Das Lösen von Differentialgleichungen gehört zu einer der Stärken von Mathematica. 
Die grundlegenden Befehle für die analytische bzw. numerische Auflösung von Differentialgleichungen lauten DSolve und NDSolve. Die Hilfe (F1) liefert Erklärungen und Beispiele zu diesen mächtigen Befehlen und ihren zahlreichen Optionen. 

Im folgenden ein oder mehrere Beispiele:

Differentialgleichung der Winkelgeschwindigkeit der Erde:  w'==((2π a b)/( T (p^2)))((( 1-ϵ Cos w ))^2) mit w(0) = 0

Grunddaten zur Umlaufbahn der Erde um die Sonne (Entfernungen in m, Zeiten in sec etc.) als Liste {...} :

{
perihel=147.1*10^6*10^3 , (* Sonnennächster Punkt *)
aphel=152.1*10^6*10^3,  (* Sonnenfernster Punkt *)
ϵ=0.0167,  (* Exzentrizität *)
T=365.2*24*60*60,  (* Umlaufzeit anomalistisch *)

a=((perihel+aphel)/2)  ,(* Grosse Halbachse *)
b=(Sqrt[(a^2)-(ϵ^2)(a^2)]) , (* Kleine Halbachse *)
p=((b^2)/a)  (* Ellipsenparameter *)
}

Die Differentialgleichung der Winkels w als Funktion der Zeit t lautet w'(t)==((2π a b)/( T (p^2)))(((1-ϵ Cos(w(t))))^2) mit Anfangsbedingung w(0) = 0. 
Diese Gleichung lässt sich nicht elementar auflösen, und auch Mathematicas sonst mächtiger Befehl DSolve zum Lösen von Differentialgleichungen liefert keine Lösung. Allfällige Warnmeldungen werden hier mittels der Funktion Quiet (postfix) unterbunden.

DSolve[{w'[t]==((2π a b)/( T (p^2)))(((1-ϵ Cos[w[t]]))^2),w[0]==0},w,t]

Deshalb wird die Gleichung numerisch gelöst (im Zeitintervall zweier Jahre):

sol=NDSolve[{w'[t]==((2π a b)/( T (p^2)))(((1-ϵ Cos[w[t]]))^2),w[0]==0},w,{t,0,2T}]

Die Lösung ist eine interpolierte Funktion (verpackt in eine Ersetzungsregel und gespeichert in der Variable sol). Solche "numerischen Objekte" sind in Mathematica Funktionen, die in vielen Aspekten wie "normale" Funktionsformeln behandelt werden können (Plotting, Ableiten, Integrieren etc.). Unten die grafische Darstellung der Winkelgeschwindigkeit w'[t] (die Befehlsfolge /.sol bedeutet die Anwendung der in der Variable sol gespeicherten Ersetzungsregel auf den Ausdruck w'[t]).

plotw=Plot[Evaluate[w'[t]/.sol],{t,0,2T}, PlotLabel->"Velocity of angle (1/sec)"];

Dieser Plot wird zusammen mit demjenigen der "überstrichenen Fahrfläche in Abhängigkeit der Zeit" dargestellt. Die zeitliche Steigung dieser Fläche ist konstant wegen der Konstanz des Drehimpulses (2. Keplersches Gesetz).

GraphicsGrid[{{Plot[Evaluate[w[t]/.sol],{t,0,2T}, PlotLabel->"Time-Angle"]},{plotw,Plot[((Pi a b t)/T), {t,0,2T}, PlotLabel->"Time-Area (m^2)"]}}]

Differentialgleichung für Fallbewegung mit Reibung aus relativ grosser Höhe: w'==-g((R^2)/(  (s^2)w))-(k/m)w  mit w(0) = (Subscript[w, 0])

Die Differentialgleichung  w'==-g((R^2)/(  (s^2)w))-((k )/m)w mit (parametrischer) Anfangsbedingung w((Subscript[s, 0])) = (Subscript[w, 0])  kann analytisch gelöst werden (aber nicht elementar). Die Resultate (Funktionen) von DSolve werden als Menge von Ersetzungsregeln (Pfeil-Operator -> ) zurückgegeben; diese werden hier in der Variable sol gespeichert, um weiter unten für einen Plot-Befehl wieder verwendet werden zu können.
Man beachte die Zweideutigkeit der allgemeinen Lösung.

sol=DSolve[{w'[s] ==-g((R^2)/(  (s^2)w[s]))-((k )/m)w[s], w[s0]==w0},w,s]//Quiet;
%//TraditionalForm

Der Reibungskoeffizient k ist in diesem System negativ interpretiert. Die Lösung beschreibt ein Geschwindigkeits-Weg-Diagramm.

Plot[{w[s]/.sol[[1]],w[s]/.sol[[2]]}/.{R->6378000.,k->-0.1,m->100.,g->9.81,s0->6378000.+10000, w0->0},{s,6378000.+10000,6378000.},PlotRange->All]

Das (explizite) Euler-Verfahren

Zeit-Winkelfunktion w'==((c( 1-ϵ Cos w ))^2) mit w(0) = 0

Um eine aussagekräftigere (elliptischere) Grafik als mit den realen Erddaten zu erhalten wird die Winkelgleichung für normierte Parameter verwendet (denn für die Erdbahn sind a und b nahezu gleich.)

f[x_,y_]:=1(((1-0.25 Cos[y]))^2);
x0=0; y0=0;epsx=8;epsy=2Pi;

demoEuler1

Das (explizite) lokale Taylor-Verfahren

Zeit-Winkelfunktion w'==((c( 1-ϵ Cos w ))^2) mit w(0) = 0

Um eine aussagekräftigere (elliptischere) Grafik als mit den realen Erddaten zu erhalten wird die Winkelgleichung für normierte Parameter verwendet (denn für die Erdbahn sind a und b nahezu gleich.)

f[x_,y_]:=1(((1-0.25 Cos[y]))^2);
x0=0; y0=0;epsx=8;epsy=2Pi;

demoLocalTaylorMethod1

y'==x (y^(1/3)) mit y(1) = 1

Um eine aussagekräftigere (elliptischere) Grafik als mit den realen Erddaten zu erhalten wird die Winkelgleichung für normierte Parameter verwendet (denn für die Erdbahn sind a und b nahezu gleich.)

f[x_,y_]:=x y^(1/3);
x0=1; y0=1;epsx=2;epsy=2;

demoLocalTaylorMethod1

y'==-100y+99(e^(-x)) mit y(0) = 0

f[x_,y_]:=-100y + 99 E^(-x);
x0=0; y0=0;epsx=2;epsy=2;

demoLocalTaylorMethod1

NDSolve Illustration for classical explicit Runge-Kutta

Special Definition of differential equation:

Needs["DifferentialEquations`NDSolveProblems`"];
Needs["DifferentialEquations`NDSolveUtilities`"];Needs["DifferentialEquations`InterpolatingFunctionAnatomy`"];

Definition of differential equation:

{a,b}={0,8}; (* Interval *)
system=NDSolveProblem[{{y'[x]==1(((1-1/4 Cos[y[x]]))^2)},{y[0]==0},{y},{x,a,b},{},{},{}}]

Reference solution by built-in high precision adaptive method:

sol = NDSolve[system,Method->"StiffnessSwitching",WorkingPrecision->32]

Classical explicit Runge-Kutta coefficients (a-b-c-numbers):

ClassicalRungeKuttaCoefficients[4,p_]:=N[{{{1/2},{0,1/2},{0,0,1}},{1/6,1/3,1/3,1/6},{1/2,1/2,1}},p];

h=0.5; (* Starting step size *)

Solution with classical explicit Runge - Kutta and fixed step size:

sol1a=NDSolve[system,StartingStepSize->h, Method->{"FixedStep",Method->{"ExplicitRungeKutta","DifferenceOrder"->4,"Coefficients"->ClassicalRungeKuttaCoefficients}}]

Solution with classical explicit Euler method and fixed step size:

sol1b=NDSolve[system,StartingStepSize->h, Method->{"FixedStep",Method->{"ExplicitEuler"}}]

Readout of coordinate data for different methods and generate plots:

coords = First[InterpolatingFunctionCoordinates[First[y/.sol1a]]];
gsol1a=ListPlot[Transpose[{coords, First[y/.sol1a][coords]}],PlotStyle->Red];

coords = First[InterpolatingFunctionCoordinates[First[y/.sol1b]]];
gsol1b=ListPlot[Transpose[{coords, First[y/.sol1b][coords]}],PlotStyle->Gray];

Show the plots inluding reference solution:

Show[Plot[y[x]/.sol,{x,a,b},PlotStyle->{Black,Dashed}],gsol1a,gsol1b]

Read out coordinate data for different methods and generate tabular data:

n=21; (* Number of steps = rows in tabular data ≤ Length of interval*)
TableForm[
Take[Transpose[{coords, First[y/.sol1a][coords],First[y/.sol1b][coords],First[y/.sol][coords]}],Min[(b-a)/h+1,n]],TableHeadings->{None,{"(Subscript[t, k])","Classical Runge-Kutta","ExplicitEuler","Reference"}},
TableSpacing->{0.2,2}]

Anhang: Demo-Codes

Der Demo-Code in demoEuler1 enthält u.a. einen numerisch-generischen Lösungsversuch für die entsprechende Differentialgleichung; dieser Versuch kann in gewissen Fällen zu unvollständigen Lösungsmengen führen.

demoEuler1:=
Module[{list={}},
Manipulate[
list=ap/.{x1_,y1_}:>NestList[{h+#[[1]],#[[2]]+h f[#[[1]],#[[2]]]}&,{x1,y1},Floor[epsx/h]];
Quiet@
Show[
If[pr,
VectorPlot[{1,Re[f[x,y]]},
{x,x0-epsx,x0+epsx},{y,y0-epsy,y0+epsy},
Axes->True, VectorPoints->pp,  VectorScale->{Automatic,0.,Automatic},VectorStyle->Darker[Gray],AxesOrigin->If[ao,{0,0},Automatic]
],
VectorPlot[{1,Re[f[x,y]]},
{x,ap[[1]]-epsx,ap[[1]]+epsx},{y,ap[[2]]-epsy,ap[[2]]+epsy},
Axes->True, VectorPoints->pp,  VectorScale->{Automatic,0.,Automatic},VectorStyle->Darker[Gray],AxesOrigin->If[ao,{0,0},Automatic]
]
],
 If[pr,Plot[Evaluate[y[x]/.(ap/.({x1_,y1_}:>NDSolve[{y'[x]==f[x,y[x]],y[x1]==y1},y,{x,x1-epsx,x1+epsx}]))],{x,x0-epsx,x0+epsx}, PlotRange->Automatic,PlotStyle->Thick],
Plot[Evaluate[y[x]/.(ap/.({x1_,y1_}:>NDSolve[{y'[x]==f[x,y[x]],y[x1]==y1},y,{x,x1-epsx,x1+epsx}]))],{x,ap[[1]]-epsx,ap[[1]]+epsx}, PlotRange->Automatic,PlotStyle->Thick],
],
Graphics[{Red,PointSize[Large],Point[ap]}],
ListLinePlot[Take[list,i+1],PlotStyle->{Thick,Red},PlotMarkers->Automatic]
(* PlotRange->If[pr,{{x0-epsx,x0+epsx},{y0-epsy,y0+epsy}},{{ap[[1]]-epsx,ap[[1]]+epsx},Automatic}] *)
],
{{pp,15,"Plotpunkte (Richtungsfeld)"},1,100,1},
Delimiter,
{{h,epsx/5,"Schrittweite h"},epsx/10^6,epsx},
{{i,0,"Anzahl Schritte i" },0,Floor[epsx/h],1},
Delimiter,
{{ao,False,"Achsenkreuz {0,0}"},{True,False},ControlPlacement->Left},
{{pr,False,"Plotbereich fix?"},{True,False},ControlPlacement->Left},
{{ap,{x0,y0},"Anfangspunkt"},{x0-epsx,y0-epsy},{x0+epsx, y0+epsy},ControlPlacement->Left}
]
]

Der Demo-Code in demoLocalTaylorMethod1 enthält u.a. einen numerisch-generischen Lösungsversuch für die entsprechende Differentialgleichung; dieser Versuch kann in gewissen Fällen zu unvollständigen Lösungsmengen führen.

demoLocalTaylorMethod1:=
Module[{maxp=10,temp,hh,tp=0,list={}},

Clear[x,y];
(* Taylor Polynomial: Delta y *)
tp=0;temp=f[x,y[x]];
Do[tp=tp+temp*hh^i/(i!);temp=D[temp,x]/.{y'[x]->f[x,y[x]]},{i,1,maxp}];
tp=CoefficientList[tp,hh];

(*  
Print["Taylorpolnom: ",N[tp]];
Print["CoefficientList: ",CoefficientList[(N[tp]),hh ]]; 
*)

Manipulate[

list=
ap/.{x1_,y1_}:>
NestList[{1.h+#[[1]],(#[[2]]+Take[tp,p+1].Table[hh^i,{i,0,p}]/.{hh->1.h,x->#[[1]],y[x]->#[[2]]})}&,{x1,y1},Floor[epsx/(1.h)]];

Quiet@
Column[
{
Row[
{Show[
If[pr,
VectorPlot[{1,Re[f[x,y]]},
{x,x0-epsx,x0+epsx},{y,y0-epsy,y0+epsy},
Axes->True, VectorPoints->pp,  VectorScale->{Automatic,0.,Automatic},VectorStyle->Darker[Gray],AxesOrigin->If[ao,{0,0},Automatic]
],
VectorPlot[{1,Re[f[x,y]]},
{x,ap[[1]]-epsx,ap[[1]]+epsx},{y,ap[[2]]-epsy,ap[[2]]+epsy},
Axes->True, VectorPoints->pp,  VectorScale->{Automatic,0.,Automatic},VectorStyle->Darker[Gray],AxesOrigin->If[ao,{0,0},Automatic]
]
],
 If[pr,Plot[Evaluate[y[x]/.(ap/.({x1_,y1_}:>NDSolve[{y'[x]==f[x,y[x]],y[x1]==y1},y,{x,x1-epsx,x1+epsx}]))],{x,x0-epsx,x0+epsx}, PlotRange->Automatic,PlotStyle->Thick],
Plot[Evaluate[y[x]/.(ap/.({x1_,y1_}:>NDSolve[{y'[x]==f[x,y[x]],y[x1]==y1},y,{x,x1-epsx,x1+epsx}]))],{x,ap[[1]]-epsx,ap[[1]]+epsx}, PlotRange->Automatic,PlotStyle->Thick],
],
Graphics[{Red,PointSize[Large],Point[ap]}],
ListPlot[Take[list,Min[i+1,Length[list]]],PlotStyle->{PointSize[Small],Red},PlotMarkers->Automatic]
(* PlotRange->If[pr,{{x0-epsx,x0+epsx},{y0-epsy,y0+epsy}},{Automatic,Automatic}] *)
],
Style[TableForm[Join[{{"x","y"}},Take[N[list,12],Min[i+1,Length[list]]]]],Red]
},"    "
],
StringForm["Δy = `1`",Take[tp,p+1].Table[hh^i,{i,0,p}]/.{hh->"h"}]
},Spacings->1
],
{{pp,15,"Plotpunkte (Richtungsfeld)"},1,100,1},
Delimiter,
{{p,1,"Ordnung"},1,maxp,1},
{{h,epsx/5.,"Schrittweite h"},epsx/10^6,epsx},
{{i,0,"Anzahl Schritte i" },0,Floor[epsx/h],1},
Delimiter,
{{ao,False,"Achsenkreuz {0,0}"},{True,False},ControlPlacement->Left},
{{pr,False,"Plotbereich fix?"},{True,False},ControlPlacement->Left},
{{ap,{x0,y0},"Anfangspunkt"},{x0-epsx,y0-epsy},{x0+epsx, y0+epsy},ControlPlacement->Left}
]
]