Modelling transparent boundary conditions on a three-bonded quantum graph

I’ve created two issue last year, but unfortunately was able to return to this problem only now.

This question is a continue to this issue.

Essentially I am now trying to apply a solution from the question above to three-bonded star graph (I am trying to apply method, suggested by @xzczd in another issue).

System of equations:

I was able to create such a code (link to pdetoode):

{lb = -20, mb = 0, rb = 20, tmax = 24.3};
func1[x_] = 2/(9*Pi)*Exp[-((x + 10)^2/9) + I*(x + 10)];
With[{u = u1[t, x]}, eq1 = I D[u, t] + 1/2 D[u, {x, 2}] == 0;
  ic1 = {u == func1[x], u == 0} /. t -> 0;
  {bcl1, bcm1, 
    bcr1} = {u == 0 /. 
     x -> lb, -3 I/2 D[u, x] + D[u, t, x] + 3 I D[u, t] /. x -> mb, 
    u == 0 /. x -> rb}];

With[{u = u2[t, x]}, eq2 = I D[u, t] + 1/2 D[u, {x, 2}] == 0;
  ic2 = {u == 0, u == 0} /. t -> 0;
  {bcl2, bcm2, 
    bcr2} = {u == 0 /. 
     x -> lb, -3 I/2 D[u, x] + D[u, t, x] + 3 I D[u, t] /. x -> mb, 
    u == 0 /. x -> rb}];

With[{u = u3[t, x]}, eq3 = I D[u, t] + 1/2 D[u, {x, 2}] == 0;
  ic3 = {u == 0, u == 0} /. t -> 0;
  {bcl3, bcm3, 
    bcr3} = {u == 0 /. 
     x -> lb, -3 I/2 D[u, x] + D[u, t, x] + 3 I D[u, t] /. x -> mb, 
    u == 0 /. x -> rb}];
(*Creating two grids, each corresponds to an edge of the graph
*)
points = 100; {gridl, gridr} = 
 Array[# &, points, #] & /@ {{lb, mb}, {mb, rb}};

difforder = 2;
ptoofunc1 = pdetoode[u1[t, x], t, gridl, difforder];
ptoofunc2 = pdetoode[u2[t, x], t, gridr, difforder];
ptoofunc3 = pdetoode[u3[t, x], t, gridr, difforder];


del = #[[2 ;; -2]] &;


ode1 = del@ptoofunc1@eq1;
ode2 = del@ptoofunc2@eq2;
ode3 = del@ptoofunc3@eq3;


odeic1 = ptoofunc1@ic1;
odeic2 = ptoofunc2@ic2;
odeic3 = ptoofunc3@ic3;


odebc1 = ptoofunc1@bcl1;
odebc2 = ptoofunc2@bcr2;
odebc3 = ptoofunc3@bcr3;


odebcm1 = ptoofunc1@bcm1 == ptoofunc2@bcm2;
odebcm2 = ptoofunc1@bcm1 == ptoofunc3@bcm3;
odebcm3 = ptoofunc2@bcm2 == ptoofunc3@bcm3;

odebc = {odebcm1, odebcm2, odebcm3, 
   With[{sf = 1}, 
    Map[sf # + D[#, t] &, {odebc1, odebc2, odebc3}, {2}]]};
sollst = NDSolveValue[{ode1, ode2, ode3, odeic1, Rest@odeic2, 
     Rest@odeic3, odebc}, {u1 /@ gridl, u2 /@ gridr, u3 /@ gridr}, {t,
      0, tmax}, MaxSteps -> Infinity] // AbsoluteTiming;

But I keep getting error:

I would be grateful for any suggestions considering my problem. thanks for your attention

UPDATE

According to @xzczd suggestion, I’ve tried to correct a few conditions:

• initial conditions are now just ic1 = u == func1[x] /. t -> 0;, ic2 = u == 0 /. t -> 0;, ic3 = u == 0 /. t -> 0;
• $psi_{11}(−0)=psi_{12}(+0)=psi_{13}(+0)$ condition is now implemented as odebcmzero = ptoofunc1@bczero1 == ptoofunc2@bczero2 == ptoofunc3@bczero3;
• only odebcm1 now present in NDSolve

The code looks like this:

{lb = -20, mb = 0, rb = 20, tmax = 24.3};
func1[x_] = 2/(9*Pi)*Exp[-((x + 10)^2/9) + I*(x + 10)];
With[{u = u1[t, x]}, eq1 = I D[u, t] + 1/2 D[u, {x, 2}] == 0;
  ic1 = u == func1[x] /. t -> 0;
  {bcl1, bcm1, bcr1, 
    bczero1} = {u == 0 /. 
     x -> lb, -3 I/2 D[u, x] + D[u, t, x] + 3 I D[u, t] /. x -> mb, 
    u == 0 /. x -> rb, u /. x -> mb}];

With[{u = u2[t, x]}, eq2 = I D[u, t] + 1/2 D[u, {x, 2}] == 0;
  ic2 = u == 0 /. t -> 0;
  {bcl2, bcm2, bcr2, 
    bczero2} = {u == 0 /. 
     x -> lb, -3 I/2 D[u, x] + D[u, t, x] + 3 I D[u, t] /. x -> mb, 
    u == 0 /. x -> rb, u /. x -> mb}];

With[{u = u3[t, x]}, eq3 = I D[u, t] + 1/2 D[u, {x, 2}] == 0;
  ic3 = u == 0 /. t -> 0;
  {bcl3, bcm3, bcr3, 
    bczero3} = {u == 0 /. 
     x -> lb, -3 I/2 D[u, x] + D[u, t, x] + 3 I D[u, t] /. x -> mb, 
    u == 0 /. x -> rb, u /. x -> mb}];
(*Creating two grids, each corresponds to an edge of the graph
*)
points = 50; {gridl, gridr} = 
 Array[# &, points, #] & /@ {{lb, mb}, {mb, rb}};

difforder = 2;
ptoofunc1 = pdetoode[u1[t, x], t, gridl, difforder];
ptoofunc2 = pdetoode[u2[t, x], t, gridr, difforder];
ptoofunc3 = pdetoode[u3[t, x], t, gridr, difforder];


del = #[[2 ;; -2]] &;


ode1 = del@ptoofunc1@eq1;
ode2 = del@ptoofunc2@eq2;
ode3 = del@ptoofunc3@eq3;


odeic1 = ptoofunc1@ic1;
odeic2 = ptoofunc2@ic2;
odeic3 = ptoofunc3@ic3;


odebc1 = ptoofunc1@bcl1;
odebc2 = ptoofunc2@bcr2;
odebc3 = ptoofunc3@bcr3;


odebcm1 = ptoofunc1@bcm1 == ptoofunc2@bcm2;
(*odebcm2 = ptoofunc1@bcm1==ptoofunc3@bcm3;
odebcm3 = ptoofunc2@bcm2==ptoofunc3@bcm3;*)

odebcmzero = 
  ptoofunc1@bczero1 == ptoofunc2@bczero2 == ptoofunc3@bczero3;

odebc = {odebcm1, odebcmzero, 
   With[{sf = 1}, 
    Map[sf # + D[#, t] &, {odebc1, odebc2, odebc3}, {2}]]};
sollst = NDSolveValue[{ode1, ode2, ode3, odeic1, Rest@odeic2, 
     Rest@odeic3, odebc}, {u1 /@ gridl, u2 /@ gridr, u3 /@ gridr}, {t,
      0, tmax}, MaxSteps -> Infinity] // AbsoluteTiming;

I’m still getting error:

but at least it outputs an interpolating functions, which is kind of progress.

FINAL SOLUTION
Thanks to @xzczd’s invaluable assistance, I was able to fix existing problems:

• removed Rest@ for initial conditions
• aded With[{sf = 1}, Map[sf # + D[#, t] &, odebcmzero, {2}] trick

So, fully-working code with the demonstration looks like this (you also should add pdetoode function at the beginning of the file):

{lb = -20, mb = 0, rb = 20, tmax = 24.3};
func1[x_] = 2/(9*Pi)*Exp[-((x + 10)^2/9) + I*(x + 10)];
With[{u = u1[t, x]}, eq1 = I D[u, t] + 1/2 D[u, {x, 2}] == 0;
  ic1 = u == func1[x] /. t -> 0;
  {bcl1, bcm1, bcr1, 
    bczero1} = {u == 0 /. 
     x -> lb, -3 I/2 D[u, x] + D[u, t, x] + 3 I D[u, t] /. x -> mb, 
    u == 0 /. x -> rb, u /. x -> mb}];

With[{u = u2[t, x]}, eq2 = I D[u, t] + 1/2 D[u, {x, 2}] == 0;
  ic2 = u == 0 /. t -> 0;
  {bcl2, bcm2, bcr2, 
    bczero2} = {u == 0 /. 
     x -> lb, -3 I/2 D[u, x] + D[u, t, x] + 3 I D[u, t] /. x -> mb, 
    u == 0 /. x -> rb, u /. x -> mb}];

With[{u = u3[t, x]}, eq3 = I D[u, t] + 1/2 D[u, {x, 2}] == 0;
  ic3 = u == 0 /. t -> 0;
  {bcl3, bcm3, bcr3, 
    bczero3} = {u == 0 /. 
     x -> lb, -3 I/2 D[u, x] + D[u, t, x] + 3 I D[u, t] /. x -> mb, 
    u == 0 /. x -> rb, u /. x -> mb}];
(*Creating two grids, each corresponds to an edge of the graph
*)
points = 100; {gridl, gridr} = 
 Array[# &, points, #] & /@ {{lb, mb}, {mb, rb}};

difforder = 2;
ptoofunc1 = pdetoode[u1[t, x], t, gridl, difforder];
ptoofunc2 = pdetoode[u2[t, x], t, gridr, difforder];
ptoofunc3 = pdetoode[u3[t, x], t, gridr, difforder];


del = #[[2 ;; -2]] &;


ode1 = del@ptoofunc1@eq1;
ode2 = del@ptoofunc2@eq2;
ode3 = del@ptoofunc3@eq3;


odeic1 = ptoofunc1@ic1;
odeic2 = ptoofunc2@ic2;
odeic3 = ptoofunc3@ic3;


odebc1 = ptoofunc1@bcl1;
odebc2 = ptoofunc2@bcr2;
odebc3 = ptoofunc3@bcr3;


odebcm1 = ptoofunc1@bcm1 == ptoofunc2@bcm2;
(*odebcm2 = ptoofunc1@bcm1==ptoofunc3@bcm3;
odebcm3 = ptoofunc2@bcm2==ptoofunc3@bcm3;*)

odebcmzero = 
  ptoofunc1@bczero1 == ptoofunc2@bczero2 == ptoofunc3@bczero3;

odebc = {odebcm1, 
   With[{sf = 1}, Map[sf # + D[#, t] &, {odebcmzero}, {2}]], 
   With[{sf = 1}, 
    Map[sf # + D[#, t] &, {odebc1, odebc2, odebc3}, {2}]]};
sollst = NDSolveValue[{ode1, ode2, ode3, odeic1, odeic2, odeic3, 
     odebc}, {u1 /@ gridl, u2 /@ gridr, u3 /@ gridr}, {t, 0, tmax}, 
    MaxSteps -> Infinity] // AbsoluteTiming;
{soll, solr1, solr2} = 
  MapThread[rebuild, {sollst[[2]], {gridl, gridr, gridr}}];
sol1 = {t, x} [Function] 
   Piecewise[{{soll[t, x], x < mb}}, solr1[t, x]];
sol2 = {t, x} [Function] 
   Piecewise[{{soll[t, x], x < mb}}, solr2[t, x]];
Manipulate[
 Plot[Abs[sol1[t, x]]^2, {x, lb, rb}, 
  AxesLabel -> {x, 
    "|[Psi]!(*SuperscriptBox[(|), (2)]), First-second bond 
propagation"}, PlotRange -> All], {{t, 0, "time"}, 0, tmax, 
  Appearance -> "Labeled"}]
Manipulate[
 Plot[Abs[sol2[t, x]]^2, {x, lb, rb}, 
  AxesLabel -> {x, 
    "|[Psi]!(*SuperscriptBox[(|), (2)]), First-third bond 
propagation"}, PlotRange -> All], {{t, 0, "time"}, 0, tmax, 
  Appearance -> "Labeled"}]

{lb = -20, mb = 0, rb = 20, tmax = 24.3};
func1[x_] = 2/(9*Pi)*Exp[-((x + 10)^2/9) + I*(x + 10)];
With[{u = u1[t, x]}, eq1 = I D[u, t] + 1/2 D[u, {x, 2}] == 0;
  ic1 = u == func1[x] /. t -> 0;
  {bcl1, bcm1, bcr1, 
    bczero1} = {u == 0 /. 
     x -> lb, -3 I/2 D[u, x] + D[u, t, x] + 3 I D[u, t] /. x -> mb, 
    u == 0 /. x -> rb, u /. x -> mb}];

With[{u = u2[t, x]}, eq2 = I D[u, t] + 1/2 D[u, {x, 2}] == 0;
  ic2 = u == 0 /. t -> 0;
  {bcl2, bcm2, bcr2, 
    bczero2} = {u == 0 /. 
     x -> lb, -3 I/2 D[u, x] + D[u, t, x] + 3 I D[u, t] /. x -> mb, 
    u == 0 /. x -> rb, u /. x -> mb}];

With[{u = u3[t, x]}, eq3 = I D[u, t] + 1/2 D[u, {x, 2}] == 0;
  ic3 = u == 0 /. t -> 0;
  {bcl3, bcm3, bcr3, 
    bczero3} = {u == 0 /. 
     x -> lb, -3 I/2 D[u, x] + D[u, t, x] + 3 I D[u, t] /. x -> mb, 
    u == 0 /. x -> rb, u /. x -> mb}];
(*Creating two grids, each corresponds to an edge of the graph
*)
points = 100; {gridl, gridr} = 
 Array[# &, points, #] & /@ {{lb, mb}, {mb, rb}};

difforder = 2;
ptoofunc1 = pdetoode[u1[t, x], t, gridl, difforder];
ptoofunc2 = pdetoode[u2[t, x], t, gridr, difforder];
ptoofunc3 = pdetoode[u3[t, x], t, gridr, difforder];


del = #[[2 ;; -2]] &;


ode1 = del@ptoofunc1@eq1;
ode2 = del@ptoofunc2@eq2;
ode3 = del@ptoofunc3@eq3;


odeic1 = ptoofunc1@ic1;
odeic2 = ptoofunc2@ic2;
odeic3 = ptoofunc3@ic3;


odebc1 = ptoofunc1@bcl1;
odebc2 = ptoofunc2@bcr2;
odebc3 = ptoofunc3@bcr3;


odebcm1 = ptoofunc1@bcm1 == ptoofunc2@bcm2;
(*odebcm2 = ptoofunc1@bcm1==ptoofunc3@bcm3;
odebcm3 = ptoofunc2@bcm2==ptoofunc3@bcm3;*)

odebcmzero = 
  ptoofunc1@bczero1 == ptoofunc2@bczero2 == ptoofunc3@bczero3;

odebc = {odebcm1, 
   With[{sf = 1}, Map[sf # + D[#, t] &, {odebcmzero}, {2}]], 
   With[{sf = 1}, 
    Map[sf # + D[#, t] &, {odebc1, odebc2, odebc3}, {2}]]};
sollst = NDSolveValue[{ode1, ode2, ode3, odeic1, odeic2, odeic3, 
     odebc}, {u1 /@ gridl, u2 /@ gridr, u3 /@ gridr}, {t, 0, tmax}, 
    MaxSteps -> Infinity] // AbsoluteTiming;
{soll, solr1, solr2} = 
  MapThread[rebuild, {sollst[[2]], {gridl, gridr, gridr}}];
sol1 = {t, x} [Function] 
   Piecewise[{{soll[t, x], x < mb}}, solr1[t, x]];
sol2 = {t, x} [Function] 
   Piecewise[{{soll[t, x], x < mb}}, solr2[t, x]];
Manipulate[
 Plot[Abs[sol1[t, x]]^2, {x, lb, rb}, 
  AxesLabel -> {x, 
    "|[Psi]!(*SuperscriptBox[(|), (2)]), First-second bond 
propagation"}, PlotRange -> All], {{t, 0, "time"}, 0, tmax, 
  Appearance -> "Labeled"}]
Manipulate[
 Plot[Abs[sol2[t, x]]^2, {x, lb, rb}, 
  AxesLabel -> {x, 
    "|[Psi]!(*SuperscriptBox[(|), (2)]), First-third bond 
propagation"}, PlotRange -> All], {{t, 0, "time"}, 0, tmax, 
  Appearance -> "Labeled"}]