Hallo Hakan,
Zunächst hat Mma ein Package `Miscellaneous`PhysicalConstants, man
braucht sie also nicht selbst mit Werten zu belegen;
Es geht um die Bewegung eines massiven relativistischen Punktteilchens
in einem gegebenen elektromagnetischen Feld. Dazu muss die Zeitableitung
des relativistischen Impulses (m v/Sqrt[1 - v^2/c^2]) gleich der
Lorentzkraft sein. Eine Veränderung des Feldes durch die bewegte Ladung
wird nicht betrachtet (vgl. Landau/Lifschitz; LB d. theoret. Physik; Bd.
2, Klass. Feldtheorie, Kap. VIII Das Feld bewegter Ladungen). Die
Gleichungen braucht man nicht selbst nach den gesuchten Funktionen
umstellen, das kann Mma machen. Man muss fast nur die Newtonsche
Grundgleichung (dp/dt = F) mit der Lorenzkraft F eintippen und den Rest
von Mma erledigen lassen:
Erst definiert man einige Vektoren; Geschwindigkeit, elektr. und magnet.
Feldstärke:
In[3]:= Remove[vV, vE, vH];
vV[t_] := {v1[t], v2[t], v3[t]};
vE[t_] := {e1[t], e1[t], e3[t]};
vH[t_] := {h1[t], h2[t], h3[t]};
Dann definiert man die Felder (sie sind gegeben bei dieser Aufgabe):
In[18]:= Clear[fieldDefs];
fieldDefs := {e1[t] -> 0, e2[t] -> 0, e3[t] -> 0.003 (* Volt/Meter *),
h1[t] -> 0, h2[t] -> 0, h3[t] -> 0.01 (* Ampere/Meter *)}
und die Anfangsbedingungen:
In[20]:= Clear[initialC];
initialC[tmin_] := {v1[tmin] == 0, v2[tmin] == 0.2, v3[tmin] == 0};
Nun generiert man die Gleichungen und fügt die Anfangbedingungen bei,
weil NDSolve[] das so haben will:
In[22]:= Clear[eqns];
eqns[m_, eC_, c_, mY_, t0_] := Join[Apply[Equal,
Flatten[Simplify[Solve[D[m vV[t]/Sqrt[1 - (vV[t] . vV[t])/c^2], t] - eC
vE[t] - eC mY Cross[vV[t], vH[t]] == {0, 0, 0}, {v1'[t], v2'[t],
v3'[t]}]]], 1] /. fieldDefs, initialC[t0]]
Dazu 2 Bemerkungen: zum einen muss man die Gleichungen neu generieren,
wenn man die fieldDefs und/oder die initalC geändert hat. Das fand ich
den einfachsten Weg, den Input für NDSolve[] zusammenzubringen - Sie
mögen andere Wege finden. Zum anderen wird mit dem Kommando das
Differentialgleichungssystem so umgestellt, dass die Ableitungen der
Geschwindigkeitskomponenten allein auf der linken Seite zu stehen kommen
- die v1'[t] etc. werden selbst als Variable angesehen, nicht etwa die
Differentialgleichungen gelöst. Man kann sich überzeugen, dass dieser
Schritt sinnvoll gelaufen ist:
In[24]:= (* allg. Gleichungscheck *)
Simplify[(D[m vV[t]/Sqrt[1 - (vV[t] . vV[t])/c^2], t] - eC vE[t] - eC mY
Cross[vV[t], vH[t]] == {0, 0, 0}) /. Flatten[Solve[D[m vV[t]/Sqrt[1 -
(vV[t] . vV[t])/c^2], t] - eC vE[t] - eC mY Cross[vV[t], vH[t]] == {0,
0, 0}, {v1'[t], v2'[t], v3'[t]}]]]
Out[25]= True
Bei Check wird in die physikalische Form der Newton(-Lorentz)-Gleichung
die Umstellung nach v1'[t], v2'[t], v3'[t] eingesetzt und geprüft, ob
das System noch identisch erfüllt ist: Okay. Gut, dann kann man jetzt
die Lösung fertigen:
In[26]:= (* find a solution, display it *)
With[{c = 1 (* SpeedOfLight Second/Meter *), eC = 1 (*
ElectronCharge/Coulomb *), m = 1/1000 (* ist die Masse zu klein,
verschwinden die Gleichungen im Chop -> richtige Faktoren ausklammern *),
mY = 1 (* VacuumPermeability/((Volt Second)/(Ampere Meter)) *),
tmin = 0 (* Second *), tmax = 100 (* Second *)},
Print[" eqns: ", eqns[m, eC, c, mY, 0]];
solution = NDSolve[eqns[m, c, eC, mY, tmin], {v1, v2, v3}, {t, tmin, tmax}];
ParametricPlot3D[Evaluate[{v1[t], v2[t], v3[t]} /. solution], {t, tmin,
tmax}, PlotRange -> All]
]
das gibt einen printout der aktuellen Gleichungen (die hängen wie gesagt
von den gegebenen Feldern und Anfangsbedingungen ab) und erzeugt dann
auch ein Bildchen, das im Anhang ist - dort sind die
Geschwindikeitskomponenten abgetragen und die Lichtgeschwindigkeit ist 1
- grösser soll dann auch keine v-Komponente werden.
Zu beachten ist hierbei, dass z.B. mit der richtigen Elektronenmasse
(~10^-30 Kilogramm) die Differentialgleichungen gleich mal verschwinden,
nur die Anfangsbedingungen übrigbleiben und NDSolve[] daraufhin
wortreich meckert. Um es mit den physikalischen Naturkonstanten zu
rechnen, muss man geeignete Faktoren ausklammern auf beiden Seiten. Das
wurde hier jetzt nicht gemacht und bleibt Ihnen überlassen. Macht man
das falsch oder passt nicht auf, dann sieht man entweder nichts wg. der
Kleinheit der Effekte oder NDSolve[] konvergiert nicht, weil
Grössenordnungen nicht zusammenpassen und es numerisch nicht mehr
funktionieren kann.
Gruss nach Berlin
Udo.
Hakan Onel wrote:
Hallo,
ich wollte mir ansehen, wie sich ein relativistisches
geladenes Teilchen in einem Raum mit elektrischen und
magnetischen Feldern verhaelt. Dazu habe ich das
Mathematica (5.1) Notebook, welches auch an diese
eMail angehangen ist, geschrieben.
Darin loese ich die dreidimensionalen Bewegungsgleichung-
en des besagten Teilchens mit Hilfe der NDSolve-Funktion
von Mathematica (5.1) komponentenweise.
Jedoch bricht Mathematica seine Berechnung mit den
Standardparametern leider bei 10000 Iterationschritten ab.
Ein Aufruf von NDSolve mit "MaxSteps -> 'ESC' inf 'ESC'"
fuehrt aber nun dazu, dass der Computer ueber Wochen
beschaeftigt ist, ohne richtig voranzukommen.
Erhoeht man den Wert der Funktion EE_z[t_] (des elektrischen
Feldes in z-Richtung) auf unrealistische Werte (also von
0.003 auf 300000) so stellt Mathematica binnen kuerzester
Zeit das richtige Ergebnis dar.
Nun meine Frage:
Kann ich diesen NDSolve Algor. irgendwie fuer die
realistischen (sich im Notebook befindenden) Werte
beschleunigen??
Oder hat jemand eine andere Idee, wie ich mein Problem
loesen kann?
Vielen Dank und Gruesse aus Berlin
Hakan
