Lösung von Differentialgleichungen

Initialisierung

SetOptions[Plot,BaseStyle->{FontFamily->"Arial",FontSize->14},PlotStyle->{Thickness[0.01]}];
SetOptions[DensityPlot,BaseStyle->{FontFamily->"Arial",FontSize->14}];

Harmonischer Oszillator

freier harmonischer Oszillator

solution0=DSolve[{y''[t]+y[t]==0,y[0]==1,y'[0]==0},y[t],t]

Plot[y[t]/.solution0,{t,0,50},PlotRange->All,Frame->True,Axes->True,FrameLabel->{"Zeit","Amplitude"}]

gedämpfter harmonischer Oszillator

DSolve[y''[t]+y'[t]/10+y[t]==0,y[t],t]
DSolve[y''[t]+0.1y'[t]+y[t]==0,y[t],t]

solution1=DSolve[{y''[t]+0.1y'[t]+y[t]==0,y[0]==1,y'[0]==0},y[t],t]

Plot[y[t]/.solution1,{t,0,50},PlotRange->All,Frame->True,Axes->True,FrameLabel->{"Zeit","Amplitude"}]

gedämpfter getriebener harmonischer Oszillator

solution2=DSolve[{y''[t]+0.1y'[t]+y[t]==a Cos[0.3 t],y[0]==1,y'[0]==0},y[t],t]

Plot[y[t]/.solution2/.a->0.1,{t,0,50},PlotRange->All,Frame->True,FrameLabel->{"Zeit","Amplitude"}]

Natürliche Oszillationsfrequenz geht nach Einschwingzeit in Frequenz der treibenden Kraft über.

DensityPlot[y[t]/.solution2,{t,0,50},{a,0,1},PlotPoints->40,FrameLabel->{"Zeit","Amplitude der treibenden Kraft"},ColorFunction->"AvocadoColors"]

parametrisch gedämpfter getriebener harmonischer Oszillator

solution3=DSolve[{y''[t]+0.1Cos[0.5 t]y'[t]+y[t]==0.1 Cos[0.3 t],y[0]==1,y'[0]==0},y[t],t]

Analytisch nicht mehr (von Mathematica) lösbar. Numerische Lösung mittels NDSolve:

solution3=NDSolve[{y''[t]+0.1Cos[0.5 t]y'[t]+y[t]==0.1Cos[0.3 t],y[0]==1,y'[0]==0},y[t],{t,0,50}]

Plot[y[t]/.solution3,{t,0,50},PlotRange->All,Frame->True,FrameLabel->{"Zeit","Amplitude"}]

Numerische Lösung der zeitunabhängigen Schrödingergleichung

Das Potential V[x] und die Wellenfunktion ψ[x] können im Ort diskretisiert werden.

ψ[x] wird dann durch ein Array  {ψ[(Subscript[x, 1])],ψ[(Subscript[x, 2])],...} dargestellt und gleiches gilt für das Potential V[x].

Die Schrittweite dx der Diskretisierung ist durch  dx = ((Subscript[x, i+1])-(Subscript[x, i])) gegeben.

Die zweite Ableitung der Wellenfunktion ψ''[x] kann in erster Ordnung durch folgenden Näherung geschrieben werden:

ψ''[x] -> ψ''[(Subscript[x, i])] = (1/(dx^2))(ψ[(Subscript[x, i+1])]-2ψ[(Subscript[x, i])]+ψ[(Subscript[x, i-1])])

Die Schrödingergleichung mit dem Hamiltonoperator H=-((ℏ^2)/(2m))((Subscript[∇, i])^2)+V[(Subscript[x, i])] als Matrixgleichung lautet

H ψ = E ψ.

Die Lösung dieser Eigenwertgleichung liefert die Energien (Eigenwerte) und die zugehörigen Wellenfunktionen (Eigenvektoren).

Clear[ℏ,m,dx,V];
size=5;
ExampleMatrix=(DiagonalMatrix[Table[-((ℏ^2)/(2m))((-2)/(dx^2))+V[x],{size}]]+DiagonalMatrix[Table[-((ℏ^2)/(2m))(1/(dx^2)),{size-1}],1]+DiagonalMatrix[Table[-((ℏ^2)/(2m))(1/(dx^2)),{size-1}],-1])//MatrixForm

Beispiel 1 - endlicher Potentialtopf

ℏ=1; 
m=1;
a=1; (* Breite des Potentialtopfs *)
V=100; (* Höhe des Potentialtopfs *)
ndivs=1000; (* Diskretisierung des Ortes*)
boundary=2a; (* wieviel Platz wird ausserhalb des Potentialwalls miteingeschlossen *)
dx=2boundary/ndivs; (* Schrittweite *)

(* 'piecewise'- Definition von Funktionen *)

Vfunc1[x_]=(Piecewise[{{0, x<-(a/2)||(a/2)<x}, {-V, -a/2≤x≤a/2}}])

Varray1=Table[Vfunc1[x],{x,-boundary,boundary,dx}];
ListPlot[Varray1,Filling->1.1Min[Varray1],Frame->True,PlotRange->All]

Definition des Differentialoperators

D2op[A_,dim_]:=DiagonalMatrix[Table[-2 A,{dim}]]+DiagonalMatrix[Table[1 A,{dim-1}],1]+DiagonalMatrix[Table[1 A,{dim-1}],-1]

Aufstellen der Matrix :

MProblem1=D2op[-((ℏ^2)/(2m (dx^2))),ndivs+1]+DiagonalMatrix[Varray1]//N;

Lösen der Eigenwertgleichung :

answers1=Eigensystem[MProblem1];
(* sortierte Eigenwerte und Eigenfunktionen *)
{evals,evects}=Sort[Transpose[answers1],#1[[1]]<#2[[1]]&]//Transpose;
(* Anzahl der gebundenen Zustände *)
nboundstates=Select[evals,#<0&]//Length

Wähle nur gebundene Zustände aus :

boundstateenergies=Take[evals,nboundstates];

ListPlot[ Table[evects[[i]]+(i-1)/2,{i,1,nboundstates}],Joined->True,Filling->Table[n->(n-1)/2,{n,1,nboundstates}],PlotRange->All,ImageSize->600,Frame->True]

Graphische Darstellung :

(* Potentialtopf *)
potplot=ListPlot[Varray1,Filling->1.1Min[Varray1],Joined->True,Frame->True,PlotRange->All,GridLines->{None,boundstateenergies},GridLinesStyle->Directive[Thickness[0.002],Red,Dashed]];

(* Eigenzustände auf Höhe der Eigenenergien *)
ScaleFactor=30;
stateplot=ListPlot[Table[ScaleFactor evects[[i]]+evals[[i]],{i,1,nboundstates}],Joined->True,Filling->Table[n->evals[[n]],{n,1,nboundstates}],PlotStyle->Red,FillingStyle->Directive[Red,Opacity[0.3]],PlotRange->All,ImageSize->600,Frame->True,Axes->False];
Show[stateplot,potplot]

Aufgabe 2 - SGL des Harmonischen Oszillators

Berechnen und zeigen Sie die Wellenfunktionen eines harmonischen Oszillatorpotentials V(x) = - 100 +  (x^2).

ℏ=1;
m=1;
a=1; (* Frequenz  *)
V=100; (* Offset*)
ndivs=1000; (* Diskretisierung des Ortes*)
boundary=5a; (* wieviel Platz wird ausserhalb  miteingeschlossen *)
dx=2boundary/ndivs; (* Schrittweite *)

D2op[A_,dim_]:=DiagonalMatrix[Table[-2 A,{dim}]]+DiagonalMatrix[Table[1 A,{dim-1}],1]+DiagonalMatrix[Table[1 A,{dim-1}],-1]

VfuncHO[x_,a_,V_]:=-V+a (x^2)

VarrayHO=Table[VfuncHO[x,a,V],{x,-boundary,boundary,dx}];
ListPlot[VarrayHO,Filling->1.02Min[VarrayHO],Frame->True,PlotRange->All]

MProblemHO=D2op[-((ℏ^2)/(2m (dx^2))),ndivs+1]+DiagonalMatrix[VarrayHO]//N;

answersHO=Eigensystem[MProblemHO];
{evalsHO,evectsHO}=Sort[Transpose[answersHO],#1[[1]]<#2[[1]]&]//Transpose;

(* nur die untersten 10 Energieniveaus *)
stateenergies=Take[evalsHO,10];

ListPlot[ Table[evectsHO[[i]]+(i-1)/2,{i,1,10}],Joined->True,Filling->Table[n->(n-1)/2,{n,1,10}],PlotRange->All,ImageSize->600,Frame->True]

potplot=ListPlot[VarrayHO,Filling->1.1Min[VarrayHO],Joined->True,Frame->True,PlotRange->All,GridLines->{None,stateenergies},GridLinesStyle->Directive[Thickness[0.002],Red],Axes->False];
ScaleFactor=5;
stateplot=ListPlot[Table[ScaleFactor evectsHO[[i]]+evalsHO[[i]],{i,1,10}],Joined->True,Filling->Table[n->evalsHO[[n]],{n,1,10}],PlotStyle->Red,FillingStyle->Directive[Red,Opacity[0.3]],PlotRange->{-101,-75},ImageSize->600,Frame->True,Axes->False];
Show[stateplot,potplot,PlotRange->{Automatic,{-V,-0.85V}}]

Äquidistante Energienivaus (bis auf numerische Fehler):

ListPlot[tab=Table[stateenergies[[i+1]]-stateenergies[[i]],{i,1,9}],PlotRange->{Automatic,Automatic},PlotStyle->PointSize[0.02],GridLines->{None,{Mean[tab]}},Frame->True,Axes->False,FrameLabel->{"i","(Subscript[E, i + 1])-(Subscript[E, i])"}]