Tag: boundary

differential equations – Solve boundary value problem with NDSolve. How to print out approximations to a solution?

I solve particular boundary-value-problem for ODE with NDSolve “Shooting” method. The case is that solution is attained very slow, that seems to mean that boundary-value-problem which is supposed to be well-defined enough, in fact, is ill at some place. So i try figure out. First step is to see concrete values of a produced approximations to a solution. What language constructs should i use for that?

Simple example. Suppose we consider particle motion in vertical plane. We throw our particle from initial point with coordinates {x0,y0} := {0,0} and initial trajectory angle 45 degrees. And try to achieve point with coordinates {x1,y1} := {1,0} by varying particle initial velocity. We don’t know two things here: initial velocity and a duration of a motion. Here is how this toy problem can be presented and solved in mathematica:

gravity = 10;
bvpsol = NDSolve(
    {
     {
      (* ODE system (5th order) *)
      {x''(u) / c(u)^2 == 0,
       y''(u) / c(u)^2 == -gravity,
       c'(u) == 0},
      (* boundary conditions (5 items) *)
      {x(0) == y(0) == 0,
       x(1) == 1,
       y(1) == 0,
       x'(0) == y'(0)}
      }
     }, {x(u), y(u), c(u)}, u, 
    Method -> {"Shooting", 
      "StartingInitialConditions" -> {x(0) == y(0) == 0, 
        x'(0) == y'(0) == 1, c(0) == 1}}) // Flatten;

{dxdu, dydu} = D({x(u), y(u)} /. bvpsol, u);
{vx0, vy0} = ({dxdu, dydu} / c(u) /. bvpsol) /. {u -> 0};
duration = c(u) /. bvpsol /. {u -> RandomReal()};

ivpsol = NDSolve({
    (* ODE system *)
    {x''(t) == 0, y''(t) == -gravity},
    (* initial values *)
    {x(0) == y(0) == 0, x'(0) == vx0, y'(0) == vy0}
    }, {x(t), y(t)}, {t, 0, duration});

ParametricPlot({x(t), y(t)} /. ivpsol, {t, 0, duration}, 
 GridLines -> Automatic, AspectRatio -> 1/2)

comprehensive parabolic trajectory

Question: Now what options or language construct should i use to see approximations which are produced NDSolve while solving boundary-value-problem?

list manipulation – How to define more general boundary conditions for ConvolutionLayer?

In Mathematica 12, I only see the option to define "PaddingSize" for a ConvolutionLayer in order to get Dirichlet boundary conditions (i.e. padding with 0’s outside of the input image). In order to obtain some other type of boundary conditions (say, for example, Neumann), my best idea so far is to wrap the layer in a NetChain and use a NetEncoder({"Function", ...}) there to extend the input in whatever way I desire. However, this is cumbersome, probably inefficient and makes it hard to use such a layer inside more complicated models. Any ideas how to proceed better?

Sample code:

The following convolution essentially computes the derivative of its 1D input:

conv = ConvolutionLayer("Biases" -> None, "Weights" -> {{{-1, 0, 1}}}, "PaddingSize" -> 1)

However, we immediately see that the Dirichlet boundary conditions provided by "PaddingSize" -> 1 introduce some noise at the beginning and ending of the output, since

dims = {1, 10};
input = ConstantArray(1, dims);
conv@input

produces {{1., 0., 0., 0., 0., 0., 0., 0., 0., -1.}} instead of the more natural {{0., 0., 0., 0., 0., 0., 0., 0., 0., 0.}} (corresponding to Neumann boundary conditions). To solve this, we define

netenc = NetEncoder({"Function", {Join({First@First@#}, Flatten@#, {Last@First@#})} &, dims + {0, 2}});

which first extends the input with the Neumann boundary conditions, and then

net = NetChain({ConvolutionLayer("Biases" -> None, "Weights" -> {{{-1, 0, 1}}})}, "Input" -> netenc);

to compute the derivative in the interior of the extended region (i.e. with no "PaddingSize"). Now net@input returns {{0., 0., 0., 0., 0., 0., 0., 0., 0., 0.}} as desired… but is there a better way to do the same?

differential equations – Solving PDE with Neumann boundary conditions

I want to solve a pde on a triangular domain, such that you can’t leave the domain (which I think would just be Neumann boundary conditions).

My pde takes the following form:

$frac{partial u}{partial x}(A(x,y)u(x,y)) + frac{partial u}{partial y}(B(x,y)u(x,y)) + frac{partial^2 u}{partial x^2}(C(x,y)u(x,y)) + frac{partial^2 u}{partial x partial y}(E(x,y)u(x,y)) + frac{partial^2 u}{partial y^2}(F(x,y)u(x,y)) = 0.$

However, I’m struggling to solve this numerically in mathematica. When I use the code below, mathematica doesn’t do anything, just returning exactly what I wrote. My code is as follows, where A,B,C,E,F are functions of x and y I have already defined:

domain = Triangle({{0, 0}, {1, 0}, {0, 1}});
eqn = D(A(x,y)*u(x, y), x) + 
    D(B(x, y)*u(x, y), y) + 
    D(C(x,y)*u(x, y), x, x) + 
    D((E(x,y))*u(x, y), x, y) + 
     D((F(x, y)*u(x, y), y, y) == 0;
s = NDSolve(eqn, NeumannValue(0, ...), u, {x, y} (Element) domain);

I also require that the integral of $u(x,y)$ over the domain be equal to 1 – is there a way of including this condition in the pde solver?

brownian motion – Limit of an integral / Boundary behaviour of a Gaussian convolution / single layer potential

Let $k(t,x)$ be the transition density of Brownian motion $$ k(t,x) := frac{1}{sqrt{2 pi t}} exp left{ frac{-x^2}{2t} right} , quad t geq 0, x in {mathbb R.}$$

Question
Let $0 < x < a$. Show that
$$ lim_{x nearrow a}int_0^t frac{a-x}{s} k(s,x-a)k(t-s,a)ds = k(t,a).$$
Can someone offer some intuition as to why this is true?

mp.mathematical physics – How to proceed in this Boundary value problem where Eigen values are calculated numerically?

While solving a boundary value problem (background provided in the Context section) I reach the following variable separated two equations ($F(x)$ and $G(y)$)

begin{eqnarray}
lambda_h F”’ – 2 lambda_h beta_h F” + left( (lambda_h beta_h – 1) beta_h – mu right) F’ + beta_h^2 F &=& 0, tag 1\
V lambda_c G”’ – 2 V lambda_c beta_c G” + left( (lambda_c beta_c – 1) V beta_c + mu right) G’ + V beta_c^2 G &=& 0,tag 2
end{eqnarray}
with some separation constant $mu in mathbb{R}$.

with boundary conditions as:

For G: $G'(0)=0, G(0)=0$ and $frac{G”(1)}{G'(1)}=beta_c$

For F: $F'(0)=0$ and $frac{F”(1)}{F'(1)}=beta_h$

The non-homogeneous condition in $F$ is: $beta_h e^{-beta_c y}G'(y)F(0)=1$

Since the boundary conditions for $(2)$ are all homogeneous, I have calculated the eigenvalues $mu_i$ for $(2)$ numerically. For some realistic parameters like βc = 0.921, λc = 1.775*10^-4, V=1, these eigen values are

0.834041, 0.845661, 0.864286, 0.888675, 0.916951, 0.94696, 0.977271, 1.0079, 1.03972, 1.07361, 1.11015,...

Now since these eigenvalues are numeric in nature, I cannot figure out how to move forward with this problem.

I know from standard PDE problems that these eigen values should be utilized to build the $G$ solution and then should be employed in the $F$ non-homogeneous condition to determine the constants. But how am I supposed to even use orthogonality here ? I would really appreciate if someone could give a step-wise way forward for such problems where the EVs are numeric.

Context

I had the following system of PDEs
$$frac{partial theta_h}{partial x}+beta_h (theta_h-theta_w) = 0 tag A$$

$$frac{partial theta_c}{partial y} + beta_c (theta_c-theta_w) = 0 tag B$$

$$lambda_h frac{partial^2 theta_w}{partial x^2} + lambda_c Vfrac{partial^2 theta_w}{partial y^2}-frac{partial theta_h}{partial x} – Vfrac{partial theta_c}{partial y} = 0 tag C$$

with the boundary conditions ($beta_h, beta_c, V, lambda_h, lambda_c$ are constants)

$$theta_w(0,y)=1, theta_w(x,0)=0$$
$$frac{partial theta_w(1,y)}{partial x}=frac{partial theta_w(x,1)}{partial y}=0$$

$$theta_h(0,y)=1, theta_c(x,0)=0$$
Using the transformation $theta_{h1}(x,y):=theta_h (x,y)-1$ (so that $theta_w(0,y)=0$ which is needed to get an additional homogeneous condition on $F$) where

begin{eqnarray}
theta_{h1}(x,y) &=& beta_h e^{-beta_h x} int e^{beta_h x} (theta_w(x,y)-1) , mathrm{d}x,\
theta_c(x,y) &=& beta_c e^{-beta_c y} int e^{beta_c y} theta_w(x,y) , mathrm{d}y.
end{eqnarray}

which are substituted in $(C)$ to get:

$$lambda_h frac{partial^2 theta_w}{partial x^2} + lambda_c V frac{partial^2 theta_w}{partial y^2} +( -beta_h – V beta_c )theta_w +beta_h^2 e^{-beta_h x} int e^{beta_h x} theta_w(x,y) mathrm{d}x + beta_c^2 e^{-beta_c y}int e^{beta_c y} theta_w(x,y)mathrm{d}y = 0 tag D$$

Using the ansatz $theta_w(x,y) = e^{-beta_h x} f(x) e^{-beta_c y} g(y)$ on $(D)$ we obtain two linear third-order ODEs ($1,2$) with constant coefficients for $F(x) := int f(x) , mathrm{d}x$ and $G(y) := int g(y) , mathrm{d}y$ given in the original question

differential equations – Boundary condition NeumannValue is not valid

I have this code, entirely copy-pasted from Heat Transfer Tutorial. It just says the Neumann boundary condition is not valid. Any idea why? I looked around the site, but none of the suggestions I have seen so far helped. It gets up to Out(30), with the output of the Manipulate(…) in the tutorial showing up fine. The problem is then to set up a basic heat transfer model, with a simple rectangle domain. Except this “basic heat transer example” does not work for me. 

Needs("NDSolve`FEM`")
ClearAll(HeatTransferModel)
HeatTransferModel(T_, X_List, k_, (Rho)_, Cp_, Velocity_, Source_) :=
Module({V, Q, a = k}, 
  V = If(Velocity === "NoFlow", 
  0, (Rho)*Cp*Velocity.Inactive(Grad)(T, X));
  Q = If(Source === "NoSource", 0, Source);
  If(FreeQ(a, _?VectorQ), a = a*IdentityMatrix(Length(X)));
    If(VectorQ(a), a = DiagonalMatrix(a));
   (*Note the-sign in the operator*)
  a = PiecewiseExpand(Piecewise({{-a, True}}));
  Inactive(Div)(a.Inactive(Grad)(T, X), X) + V - Q)

 HeatTransferModel(T(x), {x}, k, (Rho), Cp, "NoFlow", "NoSource")

 TimeHeatTransferModel(T_, TimeVar_, X_List, k_, (Rho)_, Cp_, 
   Velocity_, Source_) := (Rho)*Cp*D(T, {TimeVar, 1}) + 
   HeatTransferModel(T, X, k, (Rho), Cp, Velocity, Source)

 TimeHeatTransferModel(T(t, x), t, {x}, k, (Rho), Cp, "NoFlow", "NoSource")

left = 0;
right = 1/5;
tend = 600;
(CapitalOmega) = Line({{left}, {right}});
Subscript((Rho), air) = 
QuantityMagnitude(ThermodynamicData("Air", "Density"));
Subscript(Cp, air) = QuantityMagnitude(ThermodynamicData("Air", "IsobaricHeatCapacity"));
Subscript(k, air) =  QuantityMagnitude(ThermodynamicData("Air", "ThermalConductivity"));
SmoothedStepFunction(fmin_, fmax_, ts_, m_) := 
Function(t, (fmin*Exp(ts*m) + fmax*Exp(m*t))/(Exp(ts*m) + Exp(m*t)))
Manipulate(
Show(Plot(SmoothedStepFunction(fmin, fmax, ts, m)(t), {t, 0, tend}, 
Sequence(PlotTheme -> "Detailed", PlotLegends -> None, 
PlotLabel -> Style("Smoothed step function: R(t)", 18), 
FrameLabel -> {"t", ""}, 
Epilog -> {Red, PointSize(Medium), 
  Point({{0, fmin}, {tend, fmax}})}, PlotRange -> {-0.5, 5.5}, 
ImageSize -> Medium)), 
Graphics({{Directive(Dashed, Orange, Thick), 
  Line({{ts, -0.5}, {ts, 5.5}})}, {Directive(Red), 
  Arrowheads(0.03), Arrow({{ts + 50, 0}, {ts, -0.5}})}, 
  Text(Style("!(*SubscriptBox((f), (min)))", 15, Red), {10,fmin + 0.5}), 
  Text(Style("!(*SubscriptBox((f), (max)))", 15, Red), {585,fmax - 0.5}), 
  Text(Style("!(*SubscriptBox((t), (s)))", 15,Red), {ts + 70, 0.2})})), 
  Sequence({{fmin, 0}, -0.5, 0.5}, {{fmax, 5}, 4.5, 5.5}, {{ts, 300}, 
     200, 400}, {{m, 0.025}, 0.015, 0.035}, SaveDefinitions -> True))

(* Here is where it breaks - The above works fine. Also, ignore the 
 ugly formatting in Show(...). I just copy-pasted it, was a bit lazy*)
 (CapitalOmega)2D = Rectangle({0, 0}, {0.02, 0.01});
  parameters = {k -> 3, h -> 50, (Rho) -> 3690,Cp -> 880, 
          (CurlyEpsilon) -> 0.7, Subscript(T, hot)-> 600, Subscript(T, amb) -> 323};
  Subscript((CapitalGamma), temp) = 
  DirichletCondition(T(x, y) == Subscript(T, hot),x == 0 || x == 0.02);
  Subscript((CapitalGamma), convective) = NeumannValue(h*(Subscript(T, amb) - T(x, y)), y == 0.01);
  (Sigma) = First(UnitConvert(Quantity("StefanBoltzmannConstant")));
   Subscript((CapitalGamma), radiation) =   
    NeumannValue((CurlyEpsilon)*(Sigma)*(Subscript(T, amb)^4 - T(x, y)^4), y == 0.01);
   pde = {HeatTransferModel(T(x, y), {x, y}, k, (Rho), Cp, "NoFlow", 
  "NoSource") == 
  Subscript((CapitalGamma), convective) + 
  Subscript((CapitalGamma), radiation), 
   Subscript((CapitalGamma), temp)} /. parameters;
  Tfun = NDSolveValue(pde, T, {x, y} (Element) (CapitalOmega)2D)
  Legended(ContourPlot(Tfun(x, y), {x, y} (Element) (CapitalOmega)2D, 
       Sequence(PlotRange -> MinMax(Tfun("ValuesOnGrid")), 
       ColorFunction -> 
       ColorData({"TemperatureMap", MinMax(Tfun("ValuesOnGrid"))}), 
       ContourStyle -> Opacity(0.1), ColorFunctionScaling -> False, 
       Contours -> 30, PlotPoints -> 41, FrameLabel -> {"x", "y"}, 
       PlotLabel -> Style("Temperature Field: T(x,y)", 18), 
       AspectRatio -> Automatic, ImageSize -> 500)), 
       BarLegend({"TemperatureMap", MinMax(Tfun("ValuesOnGrid"))}, 50, 
       LegendLabel -> Style("(K)", Opacity(0.6))))

I do not know enough Mma to even try to attempt to fix this. Thanks

real analysis – Question on “the extension” of an estimate on the boundary

Let $U$ be a compact manifold in $mathbb R^3$ and consider $f:U to mathbb R_{+}$ such that

$f le C$ for all $xin {g=0}cap {x:dist(x,{g>0})ge 0 }$ for some $gin mathcal C(U)$ and $Cin mathbb R$.

The following question arose in my head and I can’t decide whether is true or not:

If $f$ was also continuous on $U$, then could we claim that $f le C$
as well on the boundary, i.e for all $xin partial {g=0}cap
{x:dist(x,{g>0})ge 0 }$? If not, then is there any regularity
assumption one could impose on $f$ in order for the estimate to hold
also on the boundary?

Any help will be much appreciated.

Thanks in advance!

plotting – Obtaining the boundary line of a RegionPlot

I have the following code

Clear(Rg3, SetRg3, Ig3, SetIg3, κ1, Setκ1, Γ, SetΓ, κ2, Setκ2, g1, Setg1, g2, Setg2, r1)
Setκ1 = 1; SetΓ = 0.01; Setκ2 = 20;
r1 = RegionPlot(SetRg3 = 0; SetIg3 = 0; 
NMG3 = {{-Γ/
  2, -I*g1, -I*Rg3 + Ig3}, {-I*g1, -κ1/
  2, -I*g2}, {-I*Rg3 - Ig3, - I*g2, -κ2/2}} /. {Rg3 -> 
  SetRg3, Ig3 -> SetIg3, κ1 -> 
  Setκ1, Γ -> 
  SetΓ, κ2 -> Setκ2, g1 -> Setg1, 
 g2 -> Setg2};
EigensysNMG3 = Eigensystem(NMG3, Cubics -> True); {Chop(Min({((Abs(
       Normalize(
         EigensysNMG3((2))((1))) /. {κ1 -> 
          Setκ1, Γ -> 
          SetΓ, κ2 -> Setκ2, 
         g3 -> Setg3})).Transpose({{1, 0, 0}})), ((Abs(
       Normalize(
         EigensysNMG3((2))((2))) /. {κ1 -> 
          Setκ1, Γ -> 
          SetΓ, κ2 -> Setκ2, 
         g3 -> Setg3})).Transpose({{1, 0, 0}})), ((Abs(

       Normalize(
         EigensysNMG3((2))((3))) /. {κ1 -> 
          Setκ1, Γ -> 
          SetΓ, κ2 -> Setκ2, 
         g3 -> Setg3})).Transpose({{1, 0, 0}}))}) >= 
 RankedMax({((Abs(
        Normalize(
          EigensysNMG3((2))((1))) /. {κ1 -> 
           Setκ1, Γ -> 
           SetΓ, κ2 -> Setκ2, 
          g3 -> Setg3})).Transpose({{1, 0, 0}}))((
    1)), ((Abs(
        Normalize(
          EigensysNMG3((2))((2))) /. {κ1 -> 
           Setκ1, Γ -> 
           SetΓ, κ2 -> Setκ2, 
          g3 -> Setg3})).Transpose({{1, 0, 0}}))((
    1)), ((Abs(
        Normalize(
          EigensysNMG3((2))((3))) /. {κ1 -> 
           Setκ1, Γ -> 
           SetΓ, κ2 -> Setκ2, 
          g3 -> Setg3})).Transpose({{1, 0, 0}}))((1))}, 
  2))}, {Setg1, 0.01, 10}, {Setg2, 0.01, 12}, PlotRange -> Full, PlotLegends -> Automatic, PlotPoints -> Automatic, PlotStyle -> Directive(Red, Opacity(0.35)), PlotRangePadding -> None, BoundaryStyle -> {Black, Thick})

which generates the following region plot
enter image description here

My goal is to only plot the line bounding the red and white region. How should I go about going that? I tried using ContourPlot and threw in the following in addition to the code above:

ContourPlot({Chop(Min({((Abs(
      Normalize(
        EigensysNMG3((2))((1))) /. {κ1 -> 
         Setκ1, Γ -> 
         SetΓ, κ2 -> Setκ2, 
        g3 -> Setg3})).Transpose({{1, 0, 0}})), ((Abs(
      Normalize(
        EigensysNMG3((2))((2))) /. {κ1 -> 
         Setκ1, Γ -> 
         SetΓ, κ2 -> Setκ2, 
        g3 -> Setg3})).Transpose({{1, 0, 0}})), ((Abs(
      Normalize(
        EigensysNMG3((2))((3))) /. {κ1 -> 
         Setκ1, Γ -> 
         SetΓ, κ2 -> Setκ2, 
        g3 -> Setg3})).Transpose({{1, 0, 0}}))}) == 
RankedMax({((Abs(
       Normalize(
         EigensysNMG3((2))((1))) /. {κ1 -> 
          Setκ1, Γ -> 
          SetΓ, κ2 -> Setκ2, 
         g3 -> Setg3})).Transpose({{1, 0, 0}}))((1)), ((Abs(
       Normalize(
         EigensysNMG3((2))((2))) /. {κ1 -> 
          Setκ1, Γ -> 
          SetΓ, κ2 -> Setκ2, 
         g3 -> Setg3})).Transpose({{1, 0, 0}}))((1)), ((Abs(

       Normalize(
         EigensysNMG3((2))((3))) /. {κ1 -> 
          Setκ1, Γ -> 
          SetΓ, κ2 -> Setκ2, 
         g3 -> Setg3})).Transpose({{1, 0, 0}}))((1))}, 
 2))}, {Setg1, 0.01, 10}, {Setg2, 0.01, 12}, PlotPoints -> Automatic)

but I was returned with a white blank plot. The inequality for the RegionPlot was a >= and I tried doing the ContourPlot for =. Appreciate any constructive help that I can take. Thanks for reading.

differential equations – Modelling transparent boundary conditions on a quantum graph with Mathematica NDSolve and finite-element method

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:
equations

Method:
enter image description here

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__) := NDSolve`FiniteDifferenceDerivative@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, 
        MapThread(
         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(Developer`ToPackedArray@#("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} = 
  MapThread(rebuild, {sollst, {gridl, gridr, gridr2}});


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

errors

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