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} =
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"}]
• Several problems I can spot at the moment: 1. Definitions of ic1, ic2, ic3 are incorrect, because we've introduced u1, u2 and u3, the code should be modified accordingly. 2. Similarly, we need to code the condition $\psi_{11}(-0)=\psi_{12}(+0)=\psi_{13}(+0)$ explicitly. 3. The most serious problem: something is wrong with the conditions at the node. It's not too hard to notice four conditions i.e. $\psi_{11}(-0)=\psi_{12}(+0)=\psi_{13}(+0)$ and one of odebcm1, odebcm2, odebcm3 already determines a solution. Apr 23 '21 at 15:58
• @xzczd, thank you for your immensely useful tips, I've tried to make corrections you suggested (considering initial conditions - in my mind it should be just one equation for each function, because at the beginning of the whole process wave is present only on the first bond for u1, and not present on two other bonds for u2 and u3). It seems to push the problem a bit further, because it outputs interpolating function, but the issue is still present. Apr 24 '21 at 9:59
• 1. Rest before odeic1 and odeic2 should be removed now. 2. We need the With[{sf = 1}, Map[sf # + D[#, t] &, …, {2}] trick for odebcmzero. Apr 24 '21 at 10:28
• @xzczd yes, that absolutely nails it! Greatest thanks for your assistance! Apr 24 '21 at 11:14

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} =
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"}]