1
$\begingroup$

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: enter image description here

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: enter image description here

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: enter image description here

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"}]
$\endgroup$
4
  • 1
    $\begingroup$ 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. $\endgroup$
    – xzczd
    Apr 23 '21 at 15:58
  • $\begingroup$ @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. $\endgroup$
    – zanhesl
    Apr 24 '21 at 9:59
  • 1
    $\begingroup$ 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. $\endgroup$
    – xzczd
    Apr 24 '21 at 10:28
  • $\begingroup$ @xzczd yes, that absolutely nails it! Greatest thanks for your assistance! $\endgroup$
    – zanhesl
    Apr 24 '21 at 11:14
2
$\begingroup$

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"}]
$\endgroup$

Your Answer

By clicking “Post Your Answer”, you agree to our terms of service, privacy policy and cookie policy

Not the answer you're looking for? Browse other questions tagged or ask your own question.