5
$\begingroup$

I am trying to solve two coupled non-linear differential equations for $F(r)$ and $h(r)$: $$ \begin{aligned} F''-F(F^2-1)/r^2- Fh^2&=0 \\ h''+2h'/r-2F^2h/r^2+\beta^2/2 h(1-h^2)&=0 \end{aligned} $$ I know the solutions' behaviour near the origin (left): $$ \begin{aligned} F&=1+a r^2+... \\ h&=b r+... \end{aligned} $$ and at the infinity (right): $$ \begin{aligned} F&=Ae^{-r}+... \\ h&=1-B e^{-\beta r}/r+... \end{aligned} $$ where $a$, $b$, $A$ and $B$ are free parameters, $\beta$ is fixed.

To solve the system of equations I do shooting forward from $r=r_1\ll 1$:

r1 = 0.01;
r2 = 10;
mid = 4;
beta=1;

profFun1 = 
 ParametricNDSolveValue[{F''[r] - F[r] (F[r]^2 - 1)/r^2 - 
     F[r] h[r]^2 == 0, 
   h''[r] + 2/r h'[r] - 2 F[r]^2 h[r]/r^2 + 
     beta^2/2 h[r] (1 - h[r]) (1 + h[r]) == 0, 
   F[r1] == 1 + a r1^2 + 1/10 (3 a^2 + b^2) r1^4, 
   F'[r1] == 2 a r1 + 4/10 (3 a^2 + b^2) r1^3, 
   h[r1] == b r1 + b/10 (4 a - beta) r1^3, 
   h'[r1] == b + 3 b/10 (4 a - beta) r1^2}, {F, h}, {r, r1, mid}, {a, 
   b}, Method -> "ExplicitRungeKutta"]

and backward from $r=r_2\gg 1$ using the substitution $u=1/r$:

profFun2 = 
 ParametricNDSolveValue[{ 
   F1''[u] + 2/u F1'[u] - F1[u] (F1[u]^2 - 1)/u^2 - 
     F1[u] h1[u]^2/u^4 == 0, 
   h1''[u] - 2 F1[u]^2 h1[u]/u^2 + 
     beta^2/2 /u^4 h1[u] (1 - h1[u]) (1 + h1[u]) == 0, 
   F1[r1] == A Exp[-1/r1], F1'[r1] == A Exp[-1/r1] 1/r1^2, 
   h1[r1] == 1 + B Exp[-beta /r1] r1, 
   h1'[r1] == B (Exp[-beta /r1] + beta Exp[-beta/r1]/r1)}, {F1, 
   h1}, {u, r1, 1/mid}, {A, B}]

So, I get the solutions (for the left and the right):

FFl[a_, b_, r_] := profFun1[a, b][[1]][r]
hhl[a_, b_, r_] := profFun1[a, b][[2]][r]
FFr[A_, B_, r_] := profFun2[A, B][[1]][1/r]
hhr[A_, B_, r_] := profFun2[A, B][[2]][1/r]

Now I need to find such a, b, A, B that in a point r=mid the functions from the left and the right are matched smoothly. For this I also create derivatives of the functions in the point r=mid

FFld[a_, b_] := D[FFl[a, b, r], r] /. r -> mid
hhld[a_, b_] := D[hhl[a, b, r], r] /. r -> mid
FFrd[A_, B_] := D[FFr[A, B, r], r] /. r -> mid
hhrd[A_, B_] := D[hhr[A, B, r], r] /. r -> mid

For matching I have these four equations (for functions and their derivatives)

{FFl[a, b, mid] == FFr[A, B, mid], hhl[a, b, mid] == hhr[A, B, mid], 
FFld[a, b] == FFrd[A, B], hhld[a, b] == hhrd[A, B]}

So, my target is to find these four parameters to match the function from the left and the right. For my further calculations, I need just functions $F(r)$ and $h(r)$ independent on the free parameters.

I have tried something like this

FindRoot[{FFl[a, b, mid] == FFr[A, B, mid], hhl[a, b, mid] == hhr[A, B, mid], 
FFld[a, b] == FFrd[A, B], hhld[a, b] == hhrd[A, B]}, {{a,-0.3}, {b,0.6},{A,1.},{B,0.05}}]

but it does not work. However, I had already done this algorithm with FindRoot for one non-linear differential equation (with two free parameters) and it worked.

I would really appreciate suggestions, how I can get these parameters. I also used Solve instead of FindRoot, it does not work too. Maybe I should "say" something additional to FindRoot for helping it. Maybe there are some better alternatives to FindRoot for such task?

I would be also very thankful for any other ideas for solving such a system of differential equations. I recently got another system with four equations and accordingly with eight free parameters. I found this idea to solve such differential equations by matching the solutions from the left and the right side in some article, but do not know how to realise it in Mathematica. It would be great for me to learn it.

Thank you very much for reading.

$\endgroup$
0

1 Answer 1

7
$\begingroup$

I would be also very thankful for any other ideas for solving such a system of differential equations.

Then why not new-in-v12 nonlinear FEM of NDSolve?:

r1 = 0;
r2 = 9;
beta = 1;

sol=NDSolveValue[{F''[r] - F[r] (F[r]^2 - 1)/r^2 - F[r] h[r]^2 == 0, 
     h''[r] + 2/r h'[r] - 2 F[r]^2 h[r]/r^2 + 
         beta^2/2 h[r] (1 - h[r]) (1 + h[r]) == 0, 
     F[r1] == 1, h[r1] == 0, F[r2] == 0, h[r2] == 1}, {F, h}, {r, r1, r2}, 
 Method -> FiniteElement, InitialSeeding -> {F[r]==1-r/r2,h[r]==r/r2}]

Plot[sol[t] // Through // Evaluate, {t, r1, r2}, PlotRange -> All]

enter image description here

Notice I've made r2 a bit smaller. With better initial guess (it's set by the InitialSeeding option) one should be able to set larger r2.

"But what's wrong with my original attempt?" There're several issues:

  1. Evaluation order isn't controlled properly. Just execute FFl[a, b, mid] == FFr[A, B, mid] and observe the output, you'll see what's wrong. Add ?NumericQ to proper positions or add Evaluated -> False to FindRoot will resolve the problem.

  2. Value of mid isn't proper, it should be set to 1, as suggested by the asymptotic solutions.

  3. You've used r1 rather than r2 in profFun2. This is equivalent to setting r2=100, which is way too large.

After correcting all these, one can obtain:

{a -> -0.340902, b -> 0.73182, A -> 3.17522, B -> -1.90734}

Aside from these critical problems, your implementation is unnecessarily verbose. The following is mine, compare it with yours carefully:

r1 = 1/100;
r2 = 10;
mid = 1;
β = 1;

eq = {F''[r] - F[r] (F[r]^2 - 1)/r^2 - F[r] h[r]^2 == 0, 
      h''[r] + 2/r h'[r] - 2 F[r]^2 h[r]/r^2 + β^2/2 h[r] (1 - h[r]) (1 + h[r]) == 0};

FL[r_] = 1 + a r^2 + 1/10 (3 a^2 + b^2) r1^4;
hL[r_] = b r + b/10 (4 a - β) r^3;
FR[r_] = A Exp[-r];
hR[r_] = 1 + B Exp[-β r] /r;

sys1 = {eq, F[r1] == FL[r1], F'[r1] == FL'[r1], h[r1] == hL[r1], h'[r1] == hL'[r1]};
sys2 = {eq, F[r2] == FR[r2], F'[r2] == FR'[r2], h[r2] == hR[r2], h'[r2] == hR'[r2]};

varmid = {F, h, F', h'}[mid] // Through

profFun1 = ParametricNDSolveValue[sys1, varmid, {r, r1, mid}, {a, b}];
profFun2 = ParametricNDSolveValue[sys2, varmid, {r, mid, r2}, {A, B}];

solrule = FindRoot[
  profFun1[a, b] - profFun2[A, B], {{a, -0.003}, {b, 0.738}, {A, 1}, {B, 0.05}}]

ListLinePlot@
 Flatten@{NDSolveValue[sys1 /. solrule, {F, h}, {r, r1, mid}], 
          NDSolveValue[sys2 /. solrule, {F, h}, {r, mid, r2}]}

enter image description here

$\endgroup$
5
  • $\begingroup$ @xzczf Very tricky solution! I tried to vary r2 but the solution changes dramatically (r2=10,11). What could be the reason, I would expect a more robust behavior. $\endgroup$ Oct 14, 2020 at 9:05
  • $\begingroup$ @UlrichNeumann Rule of thumb: using too large a value to approximate infinity when dealing with b.c. at infinity can easily cause problem. Using a better initial guess (for example, the solution obtained with r2=9) should help. $\endgroup$
    – xzczd
    Oct 14, 2020 at 9:08
  • 2
    $\begingroup$ I got a more robust result setting InitialSeeding -> { F[r] == (1 - r/r2)^2 (1 + 2 r/r2), h[r] == 1 - (1 - r/r2)^2} $\endgroup$ Oct 14, 2020 at 9:27
  • $\begingroup$ Thank you very much for helping and correction! I also did some substitution r=x/(1-x), then you don't have the infinity and your solution works well for small and large beta! Thanks againg! $\endgroup$
    – Julia
    Oct 15, 2020 at 8:19
  • 1
    $\begingroup$ @Julia As to variable change, you may want to try this: mathematica.stackexchange.com/q/80241/1871 $\endgroup$
    – xzczd
    Oct 15, 2020 at 8:26

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.