Need help solving a system of two 1st order nonlinear differential equations

Question

The original system of equations reads:

$\begin{cases} f'(r) + f(r) \left(a(r) - \frac{1}{r}\right) = 0,\\ f^2(r) + a'(r) + \frac{a(r)}{r} - 1 = 0\,, \end{cases}$

with boundary conditions $f(0) = 0\,,\,\, f(\infty) = 1\,.$ I know that it allows non-divergent nontrivial solutions $\left(0 \leq f(r) \leq 1\,,\,\, f(r \to 0) \sim r\right.,$ and approaching exponentially to $1$ at large $r$$\left. \right)$.

Due to the boundary condition at infinity I tried the following change of variables: $r = \frac{1}{x} - 1$, which transforms it to the following system:

$\begin{cases} -x^2 f'(x) + f(x) \left(a(x) - \frac{x}{1 - x}\right) = 0,\\ f^2(x) - x^2 a'(x) + a(x) \frac{x}{1 - x} - 1 = 0\,, \end{cases}$

with boundary conditions $f(0) = 1\,,\,\, f(1) = 0\,,$ which Matematica (NDSolve) can't solve:

NDSolve[{-x^2 f'[x] + f[x] (a[x] - x/(1 - x)) == 0, f[x]^2 - x^2 a'[x] + a[x] x/(1 - x) == 1,
        f[0] == 1, f[1] == 0}, {f,a}, {x, 0, 1}]

One can eliminate $a(x)$ to obtain a 2nd-order differential equation for $f(x)$:

$f^2(x) + x^4\left( \frac{f'(x)}{f(x)}\right)^2 - x^4 \frac{f''(x)}{f(x)} +x^3\left( \frac{1}{1 - x} - 2\right)\frac{f'(x)}{f(x)} - 1 = 0\,,$ for which NDSolve again has no answer:

NDSolve[{f[x]^2 + x^4 (f'[x]/f[x])^2 - x^4 f''[x]/f[x] + x^3 (1/(1 - x) - 2) f'[x]/f[x] - 1 == 0,
      f[0] == 1, f[1] == 0}, f, {r, 0, 1}]

Any suggestions how to solve this problem in Mathematica (and not Matlab)?

Answer 1

There is a problem with boundary conditions. Changing them and fixing a few typos you can get:

eq = {-f[x]^2 + f[x]^4 + 
     x^4 Derivative[1][f][x]^2 + (x^3) /(-1 + x)
       f[x] ((1 - 2 x) f'[x] - (-1 + x) x f''[x]) == 0, f[1/10] == 1, 
   f'[1/10] == -1/10};
eq // Column // TraditionalForm

s = NDSolve[eq, f, {x, 1/10, 9/10}];
Plot[Evaluate[f[x] /. s], {x, 1/10, 9/10}, PlotRange -> All]

Usually people first try to consider some expansions, approximations, etc. at the boundaries to understand the behaviour. So the system solvable - just understand it first and set up proper BCs. Perhaps there is an error somewhere in derivation of final equation that it stops corresponding to your BCs?

Answer 2

Starting from the two equations in the question,

eq1 = D[f[r], r] + f[r] (a[r] - 1/r);
eq2 = f[r]^2 + D[a[r], r] + a[r]/r - 1;

eliminate a[r].

arule = Solve[eq1 == 0, a[r]][[1, 1]];
darule = D[%, r] // Simplify;
eq3 = -f[r] (eq2 /. {arule, darule}) // Simplify
(* f[r] - f[r]^3 + Derivative[1][f][r]/r - Derivative[1][f][r]^2/f[r] + 
   Derivative[2][f][r] *)

Next, formulate the r = 0 boundary condition by expanding f[r] as a power series there and requiring that f[0] == 0.

Series[eq3 /. Thread[Rule[#, # /. f -> Function[{r}, Sum[c[i] r^i, {i, 4}]]] & /@ 
    {f[r], Derivative[1][f][r], Derivative[2][f][r]}], {r, 0, 2}] // Normal
(* c[2] + r (c[1] - c[2]^2/c[1] + 4 c[3]) + 
   r^2 (c[2] + c[2]^3/c[1]^2 - (4 c[2] c[3])/c[1] + 9 c[4]) *)

Evidently, c[1] is arbitrary, while c[2] (and c[4]) is zero.

From question 91854 there are a variety of related ways to proceed. A particularly simple approach, due to Michael E2, is as follows.

r0 = 10^-5; r1 = 50; 
sEvent = ParametricNDSolveValue[{eq3 == 0, f[r0] == f'[r0] r0, 
    f'[r0] == b, WhenEvent[f'[r] < 0, "StopIntegration"], 
    WhenEvent[f[r] > 1, "StopIntegration"]}, f, {r, r0, r1}, {b}, WorkingPrecision -> 45];
dist[b_?NumericQ] := sEvent[b]["Domain"][[1, 2]]
FindMaximum[dist[b], {b, 8529/10000}, WorkingPrecision -> 45];
{dmax, b0} = {First@%, b /. Last[%]}
(* {19.7939826947595284038626883012992096669595719, 
    0.853177865433588923966572213040425596262318851} *)

Plot[sEvent[b0][r], {r, r0, dmax}, PlotRange -> All, AxesLabel -> {r, f}]

For completeness, a[r] can be obtained from

arule[[2]] // Apart
(* 1/r - Derivative[1][f][r]/f[r] *)

Plot[1/r - sEvent[b0]'[r]/sEvent[b0][r], {r, r0, dmax}, AxesLabel -> {r, a}]