Nonlinear system of ODEs with boundary conditions

Question

I'm trying solve this problem:

g'(r) = a(r)g(r)/r, (1/r)a'(r) = g(r)^2-1

which have the following boundary conditions: a(0)=n and g(infinity)=1, where n>=1.

I have problem with the error for n=1[hence, failed for the other n values as well.]:

"Infinite expression 1/0.^2 encountered"

using DSolve directly.

A work colleague solved for the values n=1,2,5,10 using the software MAPLE. But, I don't know to use MAPLE and I would like to get these results by Mathematica. Below are the plots found by MAPLE.

Answer 1

Another way:

eps = 10^-10; end = 10;
sol[n_] := NDSolve[
{g'[r] == a[r]*g[r]/r, a'[r]/r == g[r]^2 - 1, a[eps] == n, g[end] == 1}
,{a[r], g[r]}, {r, eps, end}, Method -> {"Shooting", 
"StartingInitialConditions" -> {a[end] == eps, g[end] == (1 - eps)}}];

Plot[Evaluate[g[r] /. sol[#] & /@ {1, 2, 5, 10}], {r, eps, end},PlotRange -> All, 
PlotLegends -> {"n=1", "n=2", "n=5", "n=10"}, PlotLabels -> "g(r)"]
Plot[Evaluate[a[r] /. sol[#] & /@ {1, 2, 5, 10}], {r, eps, end}, PlotRange -> All, 
PlotLegends -> {"n=1", "n=2", "n=5", "n=10"}, PlotLabels -> "a(r)"]

For larger values e.g.end=20:

eps = 10^-20; end = 20;
sol[n_] := NDSolve["Code the same",WorkingPrecision -> 25];
Plot["Code the same"] // Quiet(*Error messages form InterpolatingFunction to Quiet *)

Answer 2

Error in boundary condition corrected

To solve this ODE system, begin by determining its behavior near r == 0.

Series[{g'[r] - (a[r] g[r])/r, a'[r]/r + 1 - g[r]^2}, {r, 0, 2}] // Normal;
Thread[Flatten@CoefficientList[r % /. a[0] -> n, r] == 0];
Solve[%, {a'[0], a''[0], a'''[0], a''''[0], g[0], g'[0], g''[0], g'''[0]}]

(* a'[0] -> 0, a''[0] -> -1, a'''[0] -> 0, a''''[0] -> 0,
   g[0] -> 0, g'[0] -> 0, g''[0] -> 0, g'''[0] -> 0  *)

In other words, g[r] vanishes to all orders at the origin for most values of n. However, if n is a positive integer, the same set of commands leaves Derivative[n][g][0] undetermined. So, g[eps] is approximately proportional to eps^n. Next, solve the equations numerically by shooting, using the method previously employed here and elsewhere on this site. (As noted in my comment above, DSolve returns unevaluated when applied to these equations.) A typical calculation for n == 5 is given by

eps = 10^-5; end = 20;
sp = ParametricNDSolveValue[{g'[r] == a[r] g[r]/r, a'[r]/r == g[r]^2 - 1, 
    a[eps] == n0, g[eps] == eps^n0 gp, 
    WhenEvent[g[r] > 101/100, {bool = 1, "StopIntegration"}], 
    WhenEvent[a[r] < 0, {bool = 0, "StopIntegration"}]}, 
    {a[r], g[r]}, {r, eps, end + 5}, {gp, wp0, n0}, 
    WorkingPrecision -> wp0, Method -> "StiffnessSwitching"];
bl = 0; bu = 2; imax = 100; wp = 75; n = 5;
Do[bool = -1; bmiddle = (bl + bu)/2; st = sp[bmiddle, wp, n]; 
    rm = (First[st] /. r -> "Domain")[[1, 2]]; 
    If[bool == 0, bl = bmiddle, bu = bmiddle]; ip = i; 
    If[bool == -1, Return[]], {i, imax}]; 
s[n] = st;

For larger values of n, wp must be increased. For instance, with n = 10, wp = 90 is required. Results for n = {1, 2, 5, 10} are plotted below.

Plot[{s[1], s[2], s[5], s[10]}, {r, eps, 10}, PlotRange -> All, AxesLabel -> {r, "a, g"}, 
    ImageSize -> Large, LabelStyle -> {Black, Bold, Medium}]

The solutions takes no more than several seconds each.