1
$\begingroup$

I would like to numerically solve the following system of coupled differential equations:

eq1=h/(2*ma) * PsiA''[x]+ VExt*PsiA[x]+ga*Abs[PsiA[x]]^2 + gab*Abs[PsiB[x]]^2*PsiA[x] == mua*PsiA[x]

eq2=h/(2*mb) * PsiB''[x]+ VExt*PsiB[x]+gb*Abs[PsiB[x]]^2 + gab*Abs[PsiA[x]]^2*PsiB[x] == mub*PsiB[x]

The following are known model parameters:

h ma mb ga gb gab

The term "VExt" is a known function of x, namely:

VExt=-P(Cos[3/2*x/L*2*Pi])^2

The boundary conditions are the periodic ones, i.e.

PsiA[x]==PsiA[x+L]

PsiB[x]==PsiB[x+L]

I would like to know how to find eigenvalues

mua mub

and eigenfunctions

PsiA[x]   PsiB[x]

Thanks in advance.

For example purposes, you can take:

h = 1; ma = 1; mb = 2; ga = 1; gab = 1; gb = 1; L = 4;
$\endgroup$
  • $\begingroup$ Indicate the particular case of functions and parameters. I will show the solution algorithm. $\endgroup$ – Alex Trounev Oct 23 '18 at 15:13
  • $\begingroup$ I've changed my post specifying the form of VExt. Concerning model parameters, what you've chosen is very good. Can you please modify your answer accordingly? $\endgroup$ – AndreaPaco Oct 23 '18 at 15:39
  • $\begingroup$ Can you also briefly explain and comment your code. In particular, I don't understand the role of bc1 and bc2. They seem to be further constraints which I have no reasons to introduce... $\endgroup$ – AndreaPaco Oct 23 '18 at 15:41
  • $\begingroup$ These are the usual conditions for such tasks. To find the eigenvalues, three boundary conditions must be specified, not two. $\endgroup$ – Alex Trounev Oct 23 '18 at 16:11
  • $\begingroup$ Periodic solutions may not exist for every potential. $\endgroup$ – Alex Trounev Oct 23 '18 at 16:45
2
$\begingroup$

I will show the solution algorithm for some set of parameters.

h = 1; ma = 1; mb = 2; ga = 1; gab = 1; gb = 1; L = 4; VExt = \
1;
eq1 = h/(2*ma)*PsiA''[x] + VExt*PsiA[x] + ga*Abs[PsiA[x]]^2 + 
    gab*Abs[PsiB[x]]^2*PsiA[x] == mua*PsiA[x];

eq2 = h/(2*mb)*PsiB''[x] + VExt*PsiB[x] + gb*Abs[PsiB[x]]^2 + 
    gab*Abs[PsiA[x]]^2*PsiB[x] == mub*PsiB[x];

 bc1 = {PsiA[0] == 1, PsiA'[0] == 0};

bc2 = {PsiB[0] == 1, PsiB'[0] == 0};

 PSA = 
 ParametricNDSolveValue[{eq1, eq2, bc1, bc2}, 
  PsiA, {x, 0, L}, {mua, mub}];
PSB = ParametricNDSolveValue[{eq1, eq2, bc1, bc2}, 
  PsiB, {x, 0, L}, {mua, mub}];

FindRoot[{PSA[mua, mub][L] == 1, 
  PSB[mua, mub][L] == 1}, {{mua, -1}, {mub, -1}}]

Out[]= {mua -> -2.84339, mub -> -1.15353}

{Plot[PSA[-2.84339, -1.15353][x], {x, 0, L}, 
  AxesLabel -> {"x", "PsiA"}], 
 Plot[PSB[-2.84339, -1.15353][x], {x, 0, L}, 
  AxesLabel -> {"x", "PsiB"}]}

fig1

In the case of a periodic potential, it is possible to weaken the conditions at the boundary and use a parametric function with 4 parameters.

h = 1; ma = 1; mb = 2; ga = 1; gab = 1; gb = 1; L = 1; 
VExt[x_] := -(Cos[3*x/L*Pi])^2;
eq1 = h/(2*ma)*PsiA''[x] + VExt[x]*PsiA[x] + ga*Abs[PsiA[x]]^2 + 
    gab*Abs[PsiB[x]]^2*PsiA[x] == mua*PsiA[x];

eq2 = h/(2*mb)*PsiB''[x] + VExt[x]*PsiB[x] + gb*Abs[PsiB[x]]^2 + 
    gab*Abs[PsiA[x]]^2*PsiB[x] == mub*PsiB[x];
bc1 = {PsiA[0] == A0, PsiA'[0] == 0};

bc2 = {PsiB[0] == B0, PsiB'[0] == 0};
PSA = ParametricNDSolveValue[{eq1, eq2, bc1, bc2}, 
  PsiA, {x, 0, L}, {mua, mub, A0, B0}, 
  Method -> {"StiffnessSwitching", "NonstiffTest" -> False}, 
  MaxSteps -> Infinity]
PSB = ParametricNDSolveValue[{eq1, eq2, bc1, bc2}, 
  PsiB, {x, 0, L}, {mua, mub, A0, B0}, 
  Method -> {"StiffnessSwitching", "NonstiffTest" -> False}, 
  MaxSteps -> Infinity]

 FindRoot[{PSA[mua, mub, 1.9, 1.7][L] == 1.9, 
  PSB[mua, mub, 1.9, 1.7][L] == 1.7}, {{mua, -18}, {mub, -8}}, 
 Method -> "Secant"]

Out[]= {mua -> -17.8395, mub -> -7.76387}

 {Plot[PSA[-17.8395, -7.76387, 1.9, 1.7][x], {x, 0, L}, 
  AxesLabel -> {"x", "PsiA"}], 
 Plot[PSB[-17.8395, -7.76387, 1.9, 1.7][x], {x, 0, L}, 
  AxesLabel -> {"x", "PsiB"}]}

fig2

$\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.