Modelling transparent boundary conditions on a quantum graph with Mathematica NDSolve and finite-element method [duplicate]

Continue to this issue

I am now trying to apply a solution from the question above to three-bonded star graph. My idea is to take the third bond at [10,20] because I am using only one function u[t,x]. Link to pdetoode function:

System of equations: Method:

But on doing calculations in Mathematica, I keep getting errors and miscalculations(as for example, dimensions of the result of NDSolve are obviously off, it should correspond to number of points):

(*
Its pdetoode function, link in the post above, you can skip itenter image description here
*)
Clear[fdd, pdetoode, tooderule, pdetoae, rebuild]
fdd[{}, grid_, value_, order_, periodic_] := value;
fdd[a__] := NDSolveFiniteDifferenceDerivative@a;

pdetoode[funcvalue_List, rest__] :=
pdetoode[(Alternatives @@ Head /@ funcvalue) @@ funcvalue[[1]],
rest];
pdetoode[{func__}[var__], rest__] :=
pdetoode[Alternatives[func][var], rest];
pdetoode[front__, grid_?VectorQ, o_Integer, periodic_: False] :=
pdetoode[front, {grid}, o, periodic];

pdetoode[func_[var__], time_, {grid : {__} ..}, o_Integer,
periodic : True | False | {(True | False) ..} : False] :=
With[{pos = Position[{var}, time][[1, 1]]},
With[{bound = #[[{1, -1}]] & /@ {grid},
pat = Repeated[_, {pos - 1}],
spacevar = Alternatives @@ Delete[{var}, pos]},
With[{coordtoindex =
Function[coord,
Piecewise[{{1, PossibleZeroQ[# - #2[[1]]]}, {-1,
PossibleZeroQ[# - #2[[-1]]]}}, All] &, {coord, bound}]]},
tooderule@
Flatten@{((u : func) |
Derivative[dx1 : pat, dt_, dx2___][(u : func)])[x1 : pat,
t_, x2___] :> (Sow@coordtoindex@{x1, x2};

fdd[{dx1, dx2}, {grid},
Outer[Derivative[dt][u@##]@t &, grid],
"DifferenceOrder" -> o,
PeriodicInterpolation -> periodic]),
inde : spacevar :>
With[{i = Position[spacevar, inde][[1, 1]]},
Outer[Slot@i &, grid]]}]]];

tooderule[rule_][pde_List] := tooderule[rule] /@ pde;
tooderule[rule_]@Equal[a_, b_] :=
Equal[tooderule[rule][a - b], 0] //.
eqn : HoldPattern@Equal[_, _] :> Thread@eqn;
tooderule[rule_][expr_] := #[[Sequence @@ #2[[1, 1]]]] & @@
Reap[expr /. rule]

pdetoae[funcvalue_List, rest__] :=
pdetoae[(Alternatives @@ Head /@ funcvalue) @@ funcvalue[[1]], rest];
pdetoae[{func__}[var__], rest__] :=
pdetoae[Alternatives[func][var], rest];

pdetoae[func_[var__], rest__] :=
Module[{t},
Function[pde, #[
pde /. {Derivative[d__][u : func][inde__] :>
Derivative[d, 0][u][inde, t], (u : func)[inde__] :>
u[inde, t]}] /. (u : func)[i__][t] :> u[i]] &@
pdetoode[func[var, t], t, rest]]

rebuild[funcarray_, grid_?VectorQ, timeposition_: 1] :=
rebuild[funcarray, {grid}, timeposition]

rebuild[funcarray_, grid_, timeposition_?Negative] :=
rebuild[funcarray, grid, Range[Length@grid + 1][[timeposition]]]

rebuild[funcarray_, grid_, timeposition_: 1] /;
Dimensions@funcarray === Length /@ grid :=
With[{depth = Length@grid},
ListInterpolation[
Transpose[
Map[DeveloperToPackedArray@#["ValuesOnGrid"] &, #, {depth}],
Append[Delete[Range[depth + 1], timeposition], timeposition]],
Insert[grid, Flatten[#][[1]]["Coordinates"][[1]],
timeposition]] &@funcarray]

(*
constants
*)
{lb = -10, mb = 0, rb = 10, rb2 = 20, tmax = 67.3};
func2[x_] = Sin[Pi (x + 10)/10]^2

With[{u = u[t, x]}, eq = I D[u, t] + 1/2 D[u, {x, 2}] == 0;
ic = {u == func2[x], u == 0, u == 0} /. t -> 0;
{bcl, bcm, bcm2, bcr,
bcr2} = {u == 0 /.
x -> lb, -3 I/2 D[u, x] + D[u, t, x] + 3 I D[u, t] /.
x -> mb, -3 I/2 D[u, x] + D[u, t, x] + 3 I D[u, t] /. x -> rb,
u == 0 /. x -> rb, u == 0 /. x -> rb2}];
(*Creating grids, each corresponds to an edge of the graph
*)
points = 25; {gridl, gridr, gridr2} =
Array[# &, points, #] & /@ {{lb, mb}, {mb, rb}, {rb, rb2}};
difforder = 2;
(*Creating ode for each edge
*)
{ptoofuncl, ptoofuncr, ptoofuncr2} =
pdetoode[u[t, x], t, #, difforder] & /@ {gridl, gridr, gridr2};

del = #[[2 ;; -2]] &;
(*Calculating ode's on both grids at each individual point
*)
{odel, oder, oder2} =
del@#@eq & /@ {ptoofuncl, ptoofuncr, ptoofuncr2};
(*Calculating initial conditions on grids
*)
{odeicl, odeicr, odeicr2} =
MapThread[#@#2 &, {{ptoofuncl, ptoofuncr, ptoofuncr2}, ic}];
(*Calculating boundary conditions on grids
*)
{odebcl, odebcr, odebcr2} =
MapThread[#@#2 &, {{ptoofuncl, ptoofuncr, ptoofuncr2}, {bcl, bcr,
bcr2}}];
(*Calculating boundary conditions at middle point
*)
odebcm = {ptoofuncl[bcm] == ptoofuncr[bcm],
ptoofuncl[bcm] == ptoofuncr2[bcm2]};

odebc = {odebcm,
With[{sf = 1},
Map[sf # + D[#, t] &, {odebcl, odebcr, odebcr2}, {2}]]};

sollst = NDSolveValue[{odel, odeicl, oder, Rest@odeicr, oder2,
Rest@odeicr2, odebc}, {u /@ gridl, u /@ gridr, u /@ gridr2}, {t,
0, tmax}, MaxSteps -> Infinity]; // AbsoluteTiming

{soll, solr, solr2} =

sol = {t, x} \[Function] Piecewise[{{soll[t, x], x < mb}}, solr[t, x]];
Manipulate[
Plot[Abs[sol[t, x]]^2, {x, lb, rb},
AxesLabel -> {x,
"|\[Psi]\!$$\*SuperscriptBox[\(|$$, $$2$$]\)"}], {{t, 0, "time"},
0, tmax, Appearance -> "Labeled"}]

DensityPlot[sol[t, x] // #, {t, 0, tmax}, {x, lb, rb},
PlotPoints -> 50, Exclusions -> None] & /@ {Re, Im} // GraphicsRow


So is my approach correct, or I should calculate it differently?

• “My idea is to take the third bond at [10,20] because I am using only one function u[t,x].” I don't think this is the correct way to go, because $\psi _{12}(+10)=0$ is not consistent with b.c. at $x=0$. Defining 3 equations is necessary. May 17, 2020 at 2:46
• @xzczd so I had better define three individual functions like u11, u12, u13 on each bond and solve the system using them? May 17, 2020 at 8:08
• Yes. BTW there's a typo in the screenshot of formula, you write 2 $\psi_{12}$ in the definition of $\psi$. May 17, 2020 at 8:34
• How is this related to the finite element method? Right now you seem to use a finite difference scheme. Are you planing on using FEM later? May 18, 2020 at 12:49