# Trouble with dependent variables in NDSolve while modelling transparent boundary conditions on a quantum graph with Mathematica

I am trying to create a model of plane wave's propagation on a quantum graph(metric graph with a differential operator, Shrodinger operator in my case, along the edges and continuity condition at the vertices) with given multi-bond transparency boundary conditions at a point 0:  Given initial and boundary conditions look like this: I'm trying to solve it with Mathematica:

initialConditions = {
u11[0, x] == Sin[Pi*(x + 10)/10]^2,
u12[0, x] == 0,
u11[t, -10] == 0,
u11[t, 0] == 0,
u12[t, 0] == 0,
u[t, x] == u11[t, x] + u12[t, x]
}
boundaryConditions = {
-I*(3/2)*((D[ u12[t, x], x] /. x -> 0) - (D[u11[t, x], x] /.
x -> 0)) + ((D[ u12[t, x], t, x] /.
x -> 0) - (D[u11[t, x], t, x] /. x -> 0)) +
3*I*((D[ u12[t, x], t] /. x -> 0) - (D[u11[t, x], t] /.
x -> 0)) == 0,
I*D[u11[t, x], t] + (1/2)*D[u11[t, x], {x, 2}] == 0,
I*D[u12[t, x], t] + (1/2)*D[u12[t, x], {x, 2}] == 0
}
uix = NDSolveValue[{initialConditions, boundaryConditions},
u, {t, 0, 100}, {x, -10, 0}]


But I keep getting this error: If someone has already experienced such troubles in similar problems, I would highly appreciate any help.

Thanks for your attention.

Update

1. All functions should be zero at x = ±10.
2. 𝜓11 is defined at [-10, 0], 𝜓12 at [0, 10].
3. I've split one equation for 𝜓 to two for 𝜓11 and 𝜓12, because otherwise, Mathematica raised an error "system is underdetermined"
4. I've put {x, -10, 0} because otherwise, I'm getting an error, which states, that 0 is not a boundary, so BC can not be defined. So I decided to split the problem into two calculations(maybe in the wrong way).
• Hello, welcome to Mathematica.SE. Then, a few questions: 1. What's the b.c. at $x=\pm 10$? 2. Are $\psi_{12}$ and $\psi_{11}$ both defined on $[-10,10]$? 3. In your code you write I*D[u11[t, x], t] + (1/2)*D[u11[t, x], {x, 2}] == 0, I*D[u12[t, x], t] + (1/2)*D[u12[t, x], {x, 2}] == 0, which is inconsistent with the one in the picture, which one is correct? 4. In the code you write {x, -10, 0}, which is again inconsistent with the picture, which one is correct? – xzczd Apr 13 '20 at 10:49
• Hello, @xzczd, thanks for the answer. 1. All functions should be zero at x = ±10. 2. 𝜓11 is defined at [-10, 0], 𝜓12 at [0, 10]. 3. I've split one equation for 𝜓 to two for 𝜓11 and 𝜓12, because otherwise, Mathematica raised an error "system is underdetermined"4. I've put {x, -10, 0} because otherwise, I'm getting an error, which states, that 0 is not a boundary, so BC can not be defined. So I decided to split the problem into two calculations(maybe in the wrong way). – zanhesl Apr 13 '20 at 13:56
• So you're solving a problem similar to this, in the sense that there's a boundary condition inside the domain? – xzczd Apr 13 '20 at 14:12
• @xzczd yes, I think in that way it is kind of similar. But the main problem with the "dependent variables in boundary conditions" remains. Maybe you can give me a piece of advice on how to reformulate my boundary condition at 0 to solve it with Mathematica? – zanhesl Apr 13 '20 at 15:03

## 1 Answer

As mentioned in the comment above, this problem is somewhat similar to this, and harder, because the b.c. at $$x=0$$ involves first order derivative of $$t$$, which will lead to bdord warning in this case. (Here is an example about dealing with bdord warning. )

So, again, let's discretize the system in $$x$$ direction all by ourselves. I'll use pdetoode for the task.

{lb = -10, mb = 0, rb = 10, tmax = 100};

With[{u = u[t, x]}, eq = I D[u, t] + 1/2 D[u, {x, 2}] == 0;
ic = {u == Sin[Pi (x + 10)/10]^2, u == 0} /. t -> 0;
{bcl, bcm, bcr} = {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}];

points = 25; {gridl, gridr} = Array[# &, points, #] & /@ {{lb, mb}, {mb, rb}};
difforder = 2;
(* Definition of pdetoode isn't included in this post,
please find it in the link above. *)
{ptoofuncl, ptoofuncr} = pdetoode[u[t, x], t, #, difforder] & /@ {gridl, gridr};

del = #[[2 ;; -2]] &;
{odel, oder} = del@#@eq & /@ {ptoofuncl, ptoofuncr};
{odeicl, odeicr} = MapThread[#@#2 &, {{ptoofuncl, ptoofuncr}, ic}];
{odebcl, odebcr} = MapThread[#@#2 &, {{ptoofuncl, ptoofuncr}, {bcl, bcr}}];
odebcm = Equal @@ (#@bcm & /@ {ptoofuncl, ptoofuncr});

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

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

{soll, solr} = MapThread[rebuild, {sollst, {gridl, gridr}}];

sol = {t, x} \[Function] Piecewise[{{soll[t, x], x < mb}}, solr[t, x]];

DensityPlot[sol[t, x] // #, {t, 0, tmax}, {x, lb, rb}, PlotPoints -> 50,
Exclusions -> None, ColorFunction -> "AvocadoColors"] & /@ {Re, Im} // GraphicsRow • Thank you so much! I'm going to try to apply this approach to a three-bonded star graph by myself. – zanhesl Apr 14 '20 at 9:27