Following the one dimensional boundary value problem here, I would like to understand the easiest way to solve a BVP for a coupled system. In the 1D case, BVP can be converted to an initial value problem, and solved using shooting method. It was easily automated using bisection method.

On the other hand, when there are more than one boundary value to guess, bisection is simply not possible. Therefore, the question is, what would be the most efficient method for such problems.

To show the problem concretely, I have the following ODE system.

\begin{align} ( \frac{\tilde \mu_1^2}{B} -1 ) \tilde \Phi_1^2 + \frac{1}{A} {\tilde \Phi_1^{\prime 2}} + \frac{1}{2}\tilde \lambda_1 \tilde \Phi_1^4 + \frac{1}{2}\tilde \lambda_{12} \tilde \Phi_1^2 \tilde \Phi_2^2 +\frac{\tilde \mu_2^2}{B} \tilde \Phi_2^2 + \frac{1}{A} {\tilde \Phi_2^{\prime 2}} - \frac{A'}{\tilde r A^2} + \frac{1}{\tilde r^2 A} - \frac{1}{\tilde r^2} &= 0, \cr ( \frac{\tilde \mu_1 ^2 }{B} + 1 )\tilde \Phi_1^2 + \frac{1}{A} \tilde \Phi_1^{\prime 2} - \frac{1}{2}\tilde \lambda_1 \tilde \Phi_1^4 - \frac{1}{2}\tilde \lambda_{12} \tilde \Phi_1^2 \tilde \Phi_2^2 + \frac{\tilde \mu_2 ^2 }{B} \tilde \Phi_2^2 + \frac{1}{A} \tilde \Phi_2^{\prime 2} -\frac{B'}{\tilde r A B} - \frac{1}{\tilde r^2 A} + \frac{1}{\tilde r^2} &= 0, \cr \frac{1}{A}{\tilde \Phi_1^{\prime\prime }} +\left (\frac{\tilde \mu_1^2}{B} + 1 \right )\tilde \Phi_1 + \tilde \Phi'_1 \left (\frac{B'}{2 AB} - \frac{A'}{2A^2} + \frac{2}{A \tilde r}\right ) - \tilde \lambda_1 \tilde \Phi_1^3 - \frac{\tilde \lambda_{12}}{2} \tilde \Phi_2^2 \tilde \Phi_1 &= 0, \cr \frac{1}{A}{\tilde \Phi_2^{\prime\prime }} +\left (\frac{\tilde \mu_2^2}{B} \right )\tilde \Phi_2 + \tilde \Phi'_2 \left (\frac{B'}{2 AB} - \frac{A'}{2A^2} + \frac{2}{A \tilde r}\right ) - \frac{\tilde \lambda_{12}}{2} \tilde \Phi_1^2 \tilde \Phi_2 &= 0, \end{align} where $\tilde \Phi_1(\tilde r), \tilde \Phi_2(\tilde r), A(\tilde r), B(\tilde r)$ are functions to be solved, with boundary condition \begin{align} A(\infty) & = 1, \cr B(\infty) & = 1, \cr \tilde \Phi_1'(0) & = 0, \cr \tilde \Phi_1(\infty) & = 0, \cr \tilde \Phi_2'(0) & = 0, \cr \tilde \Phi_2(\infty) & = 0, \end{align} So there are $four$ functions, $six$ undetermined coefficients related to all the derivatives. On the other hand, there are $four$ equations constraints, and $six$ boundary values, so it should be solvable. There should be a class of $(\tilde \mu, \tilde \Phi)$, among which I want to select the 'ground state', $i.e.$ those $\tilde \Phi$'s that do not have zeros at finite $\tilde r$.

However, FEM does not handle it as it is not a linear coupled ODE system. Using NDSolve directly takes too long and never spits out any results, although there is no error thrown out either. Any idea?


GRHunterPhi4twoS2[Omega_, Omega2_, \[CapitalLambda]1_, \[CapitalLambda]12_, xEnd1_] := 
Module[{}, Arule = {A[x] -> (1 - 2*(M[x]/x))^(-1), Derivative[1][A][x] -> 
   D[(1 - 2*(M[x]/x))^(-1), x]}; 
eq1 = -Derivative[1][M][x] + x^2*((1/2)*(Omega^2/B[x] - 1)*\[Sigma][x]^2 + (1/2)*(Omega2^2/B[x])*
       \[Sigma]2[x]^2 + (\[CapitalLambda]1/4)*\[Sigma][x]^4 + (\[CapitalLambda]12/4)*\[Sigma][x]^2*\[Sigma]2[x]^2 + Derivative[1][\[Sigma]][x]^2/2/A[x] + 
      Derivative[1][\[Sigma]2][x]^2/2/A[x]) == 0 /. Arule; 
eq2 = Derivative[1][B][x]/A[x]/B[x]/x - (1/x^2)*(1 - 1/A[x]) - (Omega^2/B[x] + 1)*\[Sigma][x]^2 - 
    (Omega2^2/B[x])*\[Sigma]2[x]^2 + (\[CapitalLambda]1/2)*\[Sigma][x]^4 + (\[CapitalLambda]12/2)*\[Sigma][x]^2*\[Sigma]2[x]^2 - 
    Derivative[1][\[Sigma]][x]^2/A[x] - Derivative[1][\[Sigma]2][x]^2/A[x] == 0 /. Arule; 
eq3 = Derivative[2][\[Sigma]][x] + (2/x + Derivative[1][B][x]/2/B[x] - Derivative[1][A][x]/2/A[x])*
     Derivative[1][\[Sigma]][x] + A[x]*((Omega^2/B[x] + 1)*\[Sigma][x] - \[CapitalLambda]1*\[Sigma][x]^3 - 
      (\[CapitalLambda]12/2)*\[Sigma]2[x]^2*\[Sigma][x]) == 0 /. Arule; 
eq4 = Derivative[2][\[Sigma]2][x] + (2/x + Derivative[1][B][x]/2/B[x] - Derivative[1][A][x]/2/A[x])*
     Derivative[1][\[Sigma]2][x] + A[x]*((Omega2^2/B[x])*\[Sigma]2[x]) - (\[CapitalLambda]12/2)*\[Sigma]2[x]*\[Sigma][x]^2 == 0 /. 
  Arule; bc = {M[xStart1] == epsilon, Derivative[1][\[Sigma]][xStart1] == epsilon, 
  Derivative[1][\[Sigma]2][xStart1] == epsilon, B[xEnd1] == 1, \[Sigma][xEnd1] == epsilon, 
  \[Sigma]2[xEnd1] == epsilon}; sollst1 = NDSolveValue[Flatten[{eq1, eq2, eq3, eq4, bc}], 
  {M, B, \[Sigma], \[Sigma]2}, {x, xStart1, xEnd1}, Method -> "StiffnessSwitching", 
  WorkingPrecision -> WorkingprecisionVar, AccuracyGoal -> accuracyVar, 
  PrecisionGoal -> precisionGoalVar, MaxSteps -> Infinity]]
precisionVar = 40; 
psiEndTolerance = 10^(-5); 
epsilon = 10^(-6); 
xStart1 = 10^(-5); 
xEnd1 = 30; 
WorkingprecisionVar = 23; 
accuracyVar = 9; 
precisionGoalVar = 9; 
tmpSol = GRHunterPhi4twoS2[2, 2, 100, 50, xEnd1]
Plot[tmpSol[[4]][tmpx], {tmpx, xStart1, xEnd1}, PlotRange -> All]

P.S: my thought is somehow to convert the above equation set into a matrix representation. Then, maybe one can make use of the built-in Eigensystem[]. In this sense, with the boundary condition fulfilled, it can also be treated as an eigenvalue problem for $(\tilde \mu_1, \tilde \Phi_1)$ and $(\tilde \mu_2, \tilde \Phi_2)$. However, the non-linear terms $\tilde \Phi^n$ make this hard as well.

  • $\begingroup$ My attempt: by setting $A(\tilde r) = 1/(1-2M(\tilde r) /\tilde r)$, I can replace $A(\infty)=1$ to $M(0) = 0$, which makes the system slightly easier to solve. Also, $B(\infty)=1$ can be relaxed to be $B(\infty) = const.$. Now the hard part is $\tilde \Phi_1(\infty)$ and $\tilde \Phi_2(\infty)$. When there is only one $\tilde \Phi$, I could tackle the boundary value problem with shooting method. I doubt the shooting method can ever be used for this case. $\endgroup$
    – Boson Bear
    Apr 18, 2018 at 23:13
  • $\begingroup$ I think this post is very useful, although I have not found a way to implement the $\tilde \Phi^4$ terms yet... $\endgroup$
    – Boson Bear
    Apr 19, 2018 at 0:48
  • $\begingroup$ Unless you provide the code you have so far, I think, it is unlikely that you will receive much help; No one feels like typing all those equations. Why do not you share your NDSolve that take too long? $\endgroup$
    – user21
    Apr 19, 2018 at 5:01
  • $\begingroup$ Thank you for the reminder, @user21. My bad. I have edited the post with the code I have so far. $\endgroup$
    – Boson Bear
    Apr 19, 2018 at 15:38
  • $\begingroup$ Update: another very useful and relevant post is here, given by bbgodfrey. Learning it now and will update if it solves the problem. Getting warm right now. $\endgroup$
    – Boson Bear
    Apr 19, 2018 at 20:21


Your Answer

By clicking “Post Your Answer”, you agree to our terms of service and acknowledge you have read our privacy policy.