4
$\begingroup$

I am trying to solve Klein-Gordon equation : $$(\frac{\partial^2}{\partial t^2}-\frac{\partial^2}{\partial x^2}+1)\psi(x,t)=0$$ in a new coordinate system where $x\to y=\frac{x}{L(t)}$.

Here is my code :

ydum1 = x/L1[t];
expr1 = 1/c^2 D[ψ[ydum1, t], {t, 2}] - 
     D[ψ[ydum1, t], {x, 2}] + (m^2 c^2)/
      h^2 ψ[ydum1, t] /. ψ[ydum1, t] -> ψ[y, t] /. 
   x -> y L1[t] // Expand
m = 1;
c = 1;
h = 1;
ω1 = 1;
L1[t_] := 2 + Sin[ω1 t];
ic = {ψ[y, 0] == 
   Sqrt[2 (m c)/(c Sqrt[m^2 c^2 + π^2 h^2/L1[0]^2] L1[0])]
     Sin[ π y], 
  D[ψ[y, t], 
     t] == (-y L1'[t]/
       L1[t] D[Sqrt[
         2 (m c)/(c Sqrt[m^2 c^2 + π^2 h^2/L1[t]^2] L1[t])]
          Sin[ π y] Exp[-I c Sqrt[m^2 c^2 + π^2 h^2/L1[t]^2]
            t], y] + 
      D[Sqrt[2 (m c)/(c Sqrt[m^2 c^2 + π^2 h^2/L1[t]^2] L1[t])]
         Sin[ π y] Exp[-I c Sqrt[m^2 c^2 + π^2 h^2/L1[t]^2]
           t], t]) /. t -> 0};
sol1 = NDSolveValue[{expr1 == 0, 
   DirichletCondition[ψ[y, t] == 0, y >= 1], ψ[0, t] == 0, 
   ic}, ψ, {y, 0, 1}, {t, 0, 10}]

Turning back to $(x,t)$ coordinate :

xsol1 = {x, t} \[Function] sol1[x/L1[t], t];

Now I am going to define charge density : $J^0(x,t) = -\frac{i}{2}(\psi^*\frac{\partial \psi}{\partial t}-\frac{\partial \psi^*}{\partial t}\psi)$

xsolC = Conjugate@*xsol1;
dxsol = Derivative[0, 1][xsol1];
dxsolC = Conjugate@*dxsol;
J0 = -I/2 (xsolC[#1, #2] dxsol[#1, #2] - 
      dxsolC[#1, #2] xsol1[#1, #2]) &;

The charge $\int_0^{L1(t)} J^0(x,t) dx =\int_0^{L1(t)}-\frac{i}{2}(\psi^*\frac{\partial \psi}{\partial t}-\frac{\partial \psi^*}{\partial t}\psi) dx$ is a conserved quantity. However, when I plot integral (I am going to do the integral manually because I don't know how to integrate J0 with code) :

n = 300;
result = 
  Table[Last@Accumulate[Array[1/n Abs@J0[#, t] &, n, {0, L1[t]}]], {t,
     0, 10, 10/n}];
ListPlot[result]

enter image description here

This value is far from conserved. I am guessing that the problem lies in the fact that $J^0(x,t) = -\frac{i}{2}(\psi^*\frac{\partial \psi}{\partial t}-\frac{\partial \psi^*}{\partial t}\psi)$ involves multiplication of 2 quantities that come from the numerical solution. So the small error in the numerical solution is magnified.

e.g. $(\text{sol} + \text{error})(\text{sol} + \text{error}) = 2\text{sol}(\text{error}) + ...$ The error is magnified by the solution.

Is there anyway I can fix this?

Reference Taken from Modern Quantum Mechanics by J.J. Sakurai 2nd ed., page 490 enter image description here

I did make a mistake by having an extra overall minus sign, but this shouldn't affect conservation of the quantity.

I further provide proof that the charge is indeed conserved : enter image description here

-----EDIT-----

I have tested this on a simpler problem where analytical solution is easily achievable in order to trace the error. You can look up the analytical solution from equation 3.41 in the following paper http://i-rep.emu.edu.tr:8080/xmlui/bitstream/handle/11129/1302/SulaimanRafea.pdf?sequence=1. I am going to pick the solution with negative exponential sign.

ϕ[x_, t_] := Sqrt[2/en] Sin[π x] Exp[-I en t];
ϕc[x_, t_] := Sqrt[2/en] Sin[π x] Exp[I en t];
j = -I (ϕc[x, t] D[ϕ[x, t], t] - 
    D[ϕc[x, t], t] ϕ[x, t])

Output:

-4 Sin[π x]^2

j is clearly conserved since it's time independent.

Numerical treatment : The equation along with its b.c. and i.c. is written in the following code :

kge = 1/c^2 D[ψ[x, t], {t, 2}] - 
   D[ψ[x, t], {x, 2}] + ψ[x, t];
ic = {ψ[x, 0] == Sqrt[2/Sqrt[1 + π^2]] Sin[ π x], 
  D[ψ[x, t], t] == 
    D[Sqrt[2/Sqrt[1 + π^2]]
       Sin[ π x] Exp[-I Sqrt[1 + π^2] t], t] /. t -> 0}
ss = NDSolveValue[{kge == 0, 
   DirichletCondition[ψ[x, t] == 0, x >= 1], ψ[0, t] == 0, 
   ic}, ψ, {x, 0, 1}, {t, 0, 10}]

Define charge density, J :

ssC = Conjugate@*ss;
dss = Derivative[0, 1][ss];
dssC = Conjugate@*dss;
J = -I/2 (ssC[#1, #2] dss[#1, #2] - dssC[#1, #2] ss[#1, #2]) &;

I am going to plot the charge (integral of J over all space)

n = 200;
λ = 
  Table[Last@Accumulate[Array[1/n Abs@J[#, t] &, n, {0, 1}]], {t, 0, 
    1, 1/n}];
ListPlot[λ]

enter image description here

And it's not conserved.

I am going to compare analytical and numerical solutions by plotting

Manipulate[
 Plot[{Re@ss[x, t], Re@ϕ[x, t]}, {x, 0, 1}, PlotRange -> {-1, 1},
   PlotLegends -> {"Numerical", "Analytical"}], {t, 0, 10}]

enter image description here

Although the numerical solution matches pretty well with the analytical one, the numerical solution kind of "jerks" as it move in time, which affects its time derivative.

Plot[Re@dss[0.5, t], {t, 0, 10}]

enter image description here

Since the charge density takes into account the time derivative $J^0(x,t) = -\frac{i}{2}(\psi^*\frac{\partial \psi}{\partial t}-\frac{\partial \psi^*}{\partial t}\psi)$, this is probably what has caused the error such that the charge is not conserved. I would greatly appreciate it if anyone can help me fix this.

-----EDIT 2-----

Back to the original problem with TensorProductGrid method :

ic2 = {\[Psi][y, 0] == 
    Sqrt[2 (m c)/(c Sqrt[m^2 c^2 + \[Pi]^2 h^2/L1[0]^2] L1[0])]
      Sin[ \[Pi] y], 
   D[\[Psi][y, t], t] == -I Sqrt[1 + \[Pi]^2/4] Sqrt[1/Sqrt[
      1 + \[Pi]^2/4]] Sin[\[Pi] y] /. t -> 0};
sol1 = NDSolveValue[{expr1 == 
    0, \[Psi][x, t] == 0 /. {{x -> 0}, {x -> 1}}, ic2}, \[Psi], {y, 0,
    1}, {t, 0, 10}]

I receive warning :

NDSolveValue::eerr: Warning: scaled local spatial error estimate of 41.164657086731026 at t = 10. in the direction of independent variable y is much greater than the prescribed error tolerance. Grid spacing with 25 points may be too large to achieve the desired accuracy or precision. A singularity may have formed or a smaller grid spacing can be specified using the MaxStepSize or MinPoints method options.

Plot3D[{Re@sol1[y, t], Im@sol1[y, t]}, {y, 0, 1}, {t, 0, 10}, 
 AxesLabel -> {y, t, \[Psi]}, PlotRange -> {-1, 1}, 
 PlotLabel -> "Wavefunction \[Psi](y,t)", 
 PlotLegends -> {"Real", "Im"}]

enter image description here

And unfortunately the solution goes wild at later time. The fact that NDSolveValue evaluates beyond y=1 is suspicious to me, which was why I added DirichletCondition[ψ[y, t] == 0, y >= 1] before.

Manipulate[
 Plot[{Re@sol1[y, t], Im@sol1[y, t]}, {y, 0, 2}, 
  PlotLabel -> "Wavefunction \[Psi](y,t)", 
  PlotRange -> {{0, 2}, {-1.5, 1.5}}, 
  PlotLegends -> {"Real", "Im"}], {t, 0, 10}]

enter image description here

Any kind of help is much appreciated.

$\endgroup$
13
  • $\begingroup$ What happens if you decrease the max step size? $\endgroup$
    – user21
    Feb 22, 2022 at 12:35
  • 1
    $\begingroup$ Are you sure the formula for charge density is correct? Can you add a reference or something? $\endgroup$
    – xzczd
    Feb 22, 2022 at 12:49
  • $\begingroup$ Echoing the previous comment of @xzczd - are you sure you should have partial derivatives w.r.t. t there and instead they should be x. $\endgroup$
    – 1729taxi
    Feb 22, 2022 at 13:44
  • $\begingroup$ @user21, by max step size, do you mean the n in my code? I have tried varying it but the shape of the curve that I get doesn't change. @xzczd, @1729taxi, Yes, the formula is correct (except for an overall minus sign), I have provided a reference in the question above. $\endgroup$ Feb 22, 2022 at 15:05
  • $\begingroup$ MaxStepSize is an option for NDSolve. $\endgroup$
    – user21
    Feb 22, 2022 at 20:29

1 Answer 1

2
$\begingroup$

OK, since both FiniteElement and TensorProductGrid are not good at handling the problem, let me show a solution using pdetoode:

(* Definition of expr1 and ic2 are the same as yours. *)
bc = ψ[x, t] == 0 /. {{x -> 0}, {x -> 1}};

domain = {0, 1}; points = 25; difforder = 2;
grid = Array[# &, points, domain];

(* Definition of pdetoode isn't included in this post,
   please find it in the link above. *)
ptoofunc = pdetoode[ψ[y, t], t, grid, difforder];
del = #[[2 ;; -2]] &;
ode = del@ptoofunc[expr1 == 0];
odeic = ptoofunc@ic2;
odebc = With[{sf = 0}, ptoofunc@diffbc[{t, 2}, sf]@bc];
sollst = NDSolveValue[{ode, odeic, odebc}, ψ /@ grid, 
                      {t, 0, 10}]; // AbsoluteTiming

solmid = rebuild[sollst, grid, 2];
sol = {x, t} |-> solmid[x/L1[t], t];

reg = DiscretizeRegion@ImplicitRegion[0 < t < 5 && 0 < x/L1[t] < 1, {x, t}]

Plot3D[{Re@sol[x, t], Im@sol[x, t]}, {x, t} ∈ reg, AxesLabel -> {x, t, ψ}, 
 PlotRange -> All, PlotLabel -> "Wavefunction ψ(x,t)", 
 PlotLegends -> Placed[{"Real", "Im"}, {After, Center}], PlotPoints -> 100]   

Mathematica graphics

solc = Conjugate@*sol;
dsol = Derivative[0, 1][sol];
dsolc = Conjugate@*dsol;
J0 = -I/2 (solc[#1, #2] dsol[#1, #2] - dsolc[#1, #2] sol[#1, #2]) &;
help[x_?NumericQ, t_] := Abs@J0[x, t]
ListPlot[Quiet@
  Table[NIntegrate[help[x, t], {x, 0, L1[t]}, 
    Method -> {Automatic, SymbolicProcessing -> 0}], {t, 0, 1, 1/20}]]

Mathematica graphics

The conservation will be better if you make points larger.

$\endgroup$
2
  • $\begingroup$ Now everything runs perfectly, can't thank you enough!! I have some questions, what is the odebc for? I saw in your other post you have set sf=100, what is this quantity? $\endgroup$ Feb 24, 2022 at 16:41
  • 1
    $\begingroup$ @ForacleFunacle I've discussed this in detail in this post: mathematica.stackexchange.com/a/127411/1871 In short, sf is amount to "ScaleFactor" sub-option of NDSolve. Since the i.c.s and b.c.s are consistent in your case, we can set it to 0. Notice odebc together with odeic[[All, {1, -1}]] is equivalent to bc. $\endgroup$
    – xzczd
    Feb 24, 2022 at 16:50

Your Answer

By clicking “Post Your Answer”, you agree to our terms of service and acknowledge you have read our privacy policy.

Not the answer you're looking for? Browse other questions tagged or ask your own question.