Chaotische Bewegung im harmonischen Oszillator

© Jürgen Schnack, 2019/10/15

Klassische Molekulardynamik - Argoncluster - Dynamik

Energien in meV,
Längen in A,
Geschwindigkeiten in A/ps,
Massen in u

Definitionen

Teilchenzahl, Potentialparameter, Plotoptionen

Teilchenzahl=4;
M=39.948`;
ϵ=10.3`;
σ=3.405`;
Rand=20.`;
PlotRand=Rand;
TMax=50;

Hamilton-Funktion und verallgemeinerte Gradienten in P und X

P und X werden als Array definiert, damit weiß Mathematica, was das für Objekte sind

Array[P[t],3*Teilchenzahl];
Array[X[t],3*Teilchenzahl];

Hamilton-Funktion in P und X

H[P_,X_]=(0.1036446/(2 M)) (∑_(i=1)^(3*Teilchenzahl))(P[i][t]*P[i][t])      +4*ϵ*
(∑_(i=1)^Teilchenzahl)((∑_(j=i+1)^Teilchenzahl)(((σ^12)/((((((X[3*i-2][t]-X[3*j-2][t]))^2)+  (((X[3*i-1][t]-X[3*j-1][t]))^2)+(((X[3*i][t]-X[3*j][t]))^2)))^6))-
((σ^6)/((((((X[3*i-2][t]-X[3*j-2][t]))^2)+  (((X[3*i-1][t]-X[3*j-1][t]))^2)+(((X[3*i][t]-X[3*j][t]))^2)))^3))))+
1.0*(∑_(i=1)^(3*Teilchenzahl))(X[i][t]*X[i][t])

Definiere Ableitungsfelder

PGrad=Table[D[H[P,X],P[i][t]],{i,1,3*Teilchenzahl}];

XGrad=Table[D[H[P,X],X[i][t]],{i,1,3*Teilchenzahl}];

Dynamik A

Differentialgleichungen

EOMP=Table[(P[i])'[t]==-XGrad[[i]],{i,1,3*Teilchenzahl}];

EOMX=Table[(X[i])'[t]==PGrad[[i]],{i,1,3*Teilchenzahl}];

Anfangsbedingungen:

Ort0={{0.0,0.0,10.0},{0.0,0.0,-15.0},{7.0,2.0,0.0},{0.0,20.0,0.0}};

(Subscript[X, 0])=Flatten[Ort0]

(Subscript[P, c])=0.0;
Impulse0=Table[(Subscript[P, c])*{RandomReal[]-0.5,RandomReal[]-0.5,RandomReal[]-0.5},{k,1,Teilchenzahl}];
(Subscript[P, 0])=Flatten[Impulse0]

ABP=Table[P[i][0]==(Subscript[P, 0])[[i]],{i,1,3*Teilchenzahl}];
ABX=Table[X[i][0]==(Subscript[X, 0])[[i]],{i,1,3*Teilchenzahl}];

Alle Gleichungen

Equations=Join[EOMP,EOMX,ABP,ABX];

Alle Impulse und Koordinaten

AlleP=Table[P[i],{i,1,3*Teilchenzahl}];
AlleX=Table[X[i],{i,1,3*Teilchenzahl}];

AlleKoordinaten=Join[AlleP,AlleX]

DGL lösen

SolutionA=NDSolve[Equations,AlleKoordinaten,{t,0,TMax},MaxSteps->Infinity];

Dynamik B

Differentialgleichungen

EOMP=Table[(P[i])'[t]==-XGrad[[i]],{i,1,3*Teilchenzahl}];

EOMX=Table[(X[i])'[t]==PGrad[[i]],{i,1,3*Teilchenzahl}];

Anfangsbedingungen:

Ort0={{0.0,0.0,10.0+1.0},{0.0,0.0,-15.0},{7.0,2.0,0.0},{0.0,20.0,0.0}};

(Subscript[X, 0])=Flatten[Ort0]

(Subscript[P, c])=0.0;
Impulse0=Table[(Subscript[P, c])*{RandomReal[]-0.5,RandomReal[]-0.5,RandomReal[]-0.5},{k,1,Teilchenzahl}];
(Subscript[P, 0])=Flatten[Impulse0]

ABP=Table[P[i][0]==(Subscript[P, 0])[[i]],{i,1,3*Teilchenzahl}];
ABX=Table[X[i][0]==(Subscript[X, 0])[[i]],{i,1,3*Teilchenzahl}];

Alle Gleichungen

Equations=Join[EOMP,EOMX,ABP,ABX];

Alle Impulse und Koordinaten

AlleP=Table[P[i],{i,1,3*Teilchenzahl}];
AlleX=Table[X[i],{i,1,3*Teilchenzahl}];

AlleKoordinaten=Join[AlleP,AlleX]

DGL lösen

SolutionB=NDSolve[Equations,AlleKoordinaten,{t,0,TMax},MaxSteps->Infinity];

Plots

Bewegung darstellen

PlotOrte=Join[Table[ParametricPlot3D[Table[X[3 (i-1)+j][t]/.SolutionA[[1]],{j,1,3}],{t,0.1`,TMax},DisplayFunction->Identity,PlotRange->{{-Rand,Rand},{-Rand,Rand},{-Rand,Rand}},
BaseStyle->{FontFamily->"Helvetica",FontSize->18},PlotStyle->Red],{i,1,Teilchenzahl}],Table[ParametricPlot3D[Table[X[3 (i-1)+j][t]/.SolutionB[[1]],{j,1,3}],{t,0.1`,TMax},DisplayFunction->Identity,PlotRange->{{-Rand,Rand},{-Rand,Rand},{-Rand,Rand}},
BaseStyle->{FontFamily->"Helvetica",FontSize->18},PlotStyle->Blue],{i,1,Teilchenzahl}]];
FirstFrame=Show[PlotOrte,DisplayFunction->$DisplayFunction]

ColFrames=Table[FirstFrame,{i,1,TMax}];
For[i=1,i<TMax,i++,
ti=i*1.0;
ColFrames[[i]]=Show[Graphics3D[Join[
Table[{PointSize[0.05],Hue[0.0],Point[{X[3*(j-1)+1][ti]/.SolutionA[[1]],
X[3*(j-1)+2][ti]/.SolutionA[[1]],X[3*(j-1)+3][ti]/.SolutionA[[1]]}]},{j,1,Teilchenzahl}],
Table[{PointSize[0.05],Hue[0.60],Point[{X[3*(j-1)+1][ti]/.SolutionB[[1]],
X[3*(j-1)+2][ti]/.SolutionB[[1]],X[3*(j-1)+3][ti]/.SolutionB[[1]]}]},{j,1,Teilchenzahl}],
{PointSize[0.00001],Hue[0.1],Point[{-PlotRand,-PlotRand,-PlotRand}]},
{PointSize[0.00001],Hue[0.1],Point[{PlotRand,PlotRand,PlotRand}]}],Axes->True]]];

ListAnimate[ColFrames]

Distanz (Problem Ljapunov-Exponent)

DistX=Table[{ti,Sum[(N[X[3 (k-1)+j][ti]/.SolutionA[[1]]]-N[X[3 (k-1)+j][ti]/.SolutionB[[1]]])^2,{j,1,3},{k,1,Teilchenzahl}]},{ti,1,TMax}];
ListPlot[DistX]

DistXP=Table[{ti,Sum[(N[X[3 (k-1)+j][ti]/.SolutionA[[1]]]-N[X[3 (k-1)+j][ti]/.SolutionB[[1]]])^2,{j,1,3},{k,1,Teilchenzahl}]+
Sum[(N[P[3 (k-1)+j][ti]/.SolutionA[[1]]]-N[P[3 (k-1)+j][ti]/.SolutionB[[1]]])^2,{j,1,3},{k,1,Teilchenzahl}]},{ti,1,TMax}];
ListPlot[DistXP]