E .  Gilli
01 30 864

\(\(\(
\)\(Needs["LinearAlgebra`Tridiagonal`"]

    BeginPackage["pertubation`", {"LinearAlgebra`MatrixManipulation`", "LinearAlgebra`Tridiagonal`"}]\)\)\)

Clear["pertubation`*"];
RunCalc;
(Subscript[ω, 0]) = 1;
(Subscript[ω, 1]) = 0.1;
(Subscript[ω, 2]) = 0.2;
Δt = 0.01;
tmax = 40;
dist0[n_] := n + 1/2;
f2[nt_] := 
        Cos[(Subscript[ω, 2]) nt*Δt];
M = 10;
initialstate[mm_]  :=  
        Table[(Subscript[δ, 1, i]), {i, 1, 
            mm}];  
PlotX;
PlotP;
PlotES[nes_];
PlotClass;
Begin["`Private`"]
Clear["pertubation`Private`*"]
\(RunCalc\  := \((
	ℏ = 1; 
	
        m = 1; 
	δ\_\(a_, b_\) := KroneckerDelta[a, b]; 
	
        Hmat[nt_] := \ 
          Table[ℏ\ \(ω\_0\) \((dist0[
                    n])\)\ δ\_\(m, n\) + ℏ\ ω\_1\ f2[
                  nt]\ δ\_\(Abs[m - n], 1\)\ Sqrt[Max[{m, n}]], {n, 0, 
              M - 1}, {m, 0, M - 1}]; 
	
        Mat1\  = \ 
          IdentityMatrix[M] + I/ℏ*
              Hmat[nt - 1] Δt/2; 
	
        Mat2\  = 
          IdentityMatrix[M] - I/ℏ*
              Hmat[nt - 1] Δt/2; 
	
        a1 = Table[Mat1[\([n, n + 1]\)], {n, 1, M - 1}]; 
	
        a2 = Table[Mat1[\([n, n]\)], {n, 1, M}]; 
	
        a3 = Table[Mat1[\([n + 1, n]\)], {n, 1, M - 1}]; 
	
        Clear[cu]; 
\ 
        cu[0] = \ initialstate[M];  (*Anfangszustand*) 
	
        cu[nnt_Integer] := \(cu[nnt] = \((nt = nnt; \ 
              res = \ TridiagonalSolve[a1, a2, a3, Mat2 . cu[nt - 1]]; 
              nt =. ; \ res)\)\); 
	Clear[c]; 
	
        c[nt_Integer] := \(c[nt] = 
            cu[nt]/Sqrt[cu[nt] . Conjugate[cu[nt]]]\); 
	
        Table[c[100*i], {i, 1, 
            Ceiling[tmax/\((100*Δt)\)]}]; 
	
        P[nt_] := Abs[c[nt]^2]; 
	
        PT[t_] := P[Round[t/Δt]]; 
	
        a\_\(m_, n_\) := Sqrt[Max[m, n]]\ δ\_\(m + 1, n\); 
	
        adagger\_\(m_, n_\) := 
          Sqrt[Max[m, n]] δ\_\(n + 1, m\); 
	
        Ma = Table[
            a\_\(m, n\), {m, 0, M - 1}, {n, 0, M - 1}]; 
	
        Madagger = 
          Table[adagger\_\(m, n\), {m, 0, M - 1}, {n, 0, 
              M - 1}]; 
	
        XMat = Sqrt[ℏ/\((2\ m\ ω\_0)\)] \((Ma + 
                Madagger)\); 
	ClearAll[XM, XP]; 
	
        XM[nt_] := \(XM[nt] = Conjugate[c[nt]] . XMat . c[nt]\); 
	
        Table[XM[100*i], {i, 1, tmax/\((100*Δt)\)}]; 
	
        XP[t_] := Chop[XM[Round[t/Δt]]]; 
	
        Clear[x, e0, q]; 
	q = 1; 
	
        e0 = ω\_1/q\ *Sqrt[ℏ\ 2\ m\ ω\_0]; 
	{x[t_]} = 
          x[t] /. DSolve[{\(x''\)[t] + \(ω\_0\^2\) 
                      x[t] == \(-q\)\ e0/m\ Cos[ω\_2\ t], 
                x[0] == 0, \(x'\)[0] == 0}, x[t], 
              t]; 
	
        Imp\  = \ I*
            Chop[Hmat[t] . XMat - XMat . Hmat[t]]; 
	
        ImpVal[t_] := \(ImpVal[t] = 
            Conjugate[c[Round[t\/Δt]]] . Imp . 
              c[Round[t\/Δt]]\); "Ready!")\)\)
\(PlotES[
        nes_] := \((GA = 
          GraphicsArray[
            Table[{Plot[\(PT[t]\)[\([k]\)], {t, 0, tmax}, 
                  DisplayFunction -> Identity, PlotLabel -> k - 1], 
                Plot[\(PT[t]\)[\([k + 1]\)], {t, 0, tmax}, 
                  DisplayFunction -> Identity, 
                  PlotLabel -> k]}, {k, 1, nes, 2}]]; 

        Show[GA])\)\)
PlotX := 
      Plot[XP[t], {t, 0, tmax}, AxesLabel -> {"t", "<X(t)>"}, 
        PlotLabel -> "Zeitlicher Verlauf des Erwartungswertes:"]
PlotClass := 
      Plot[{x[t] /. ω2 -> 0.4, XP[t]}, {t, 0, tmax}, 
        PlotStyle -> {RGBColor[1, 0, 0], RGBColor[0, 0, 1]}]
PlotP := 
      Plot[{ImpVal[t], XP[t]}, {t, 0, tmax}, 
        PlotStyle -> {RGBColor[1, 0, 0], RGBColor[0, 0, 1]}, 
        PlotLabel -> "Zeitlicher Verlauf des Impulses", 
        AxesLabel -> {"t", "<P(t)>"}]
End[]
EndPackage[]

Anwendungsbeispiel:

Aufrufen des Packages mittels: <<pertubation`, bzw Needs["pertubation`"]

Dannach Definition der Variablen (Optional) (Defaultwerte):

(Subscript[ω, 0]) = 
        1;    
(Subscript[ω, 1]) = 0.1;    
(Subscript[ω, 2]) = 
      0.2;    (*Frequenz der Stoerkraft*) 
Δt = 0.01;
tmax = 40;
dist0[n_] := n + 1/2;
M = 10;

Durchfuehrung der Rechnung:

RunCalc

Plotfunktionen:

PlotX

PlotP

PlotClass

PlotES[10]