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$
11
  • $\begingroup$ Indicate the particular case of functions and parameters. I will show the solution algorithm. $\endgroup$ Oct 23, 2018 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, 2018 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, 2018 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$ Oct 23, 2018 at 16:11
  • $\begingroup$ Periodic solutions may not exist for every potential. $\endgroup$ Oct 23, 2018 at 16:45

1 Answer 1

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 and acknowledge that you have read and understand our privacy policy and code of conduct.

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