# Numerical solution of coupled ODEs with boundary conditions

I have to solve the following set of ODEs and just can't get good results using Mathematica $$r\frac{d}{dr}\left(\frac{1}{r}\frac{d}{dr}A(r)\right)-\xi^2F(r)^2\left(A(r)-1\right)=0$$ $$\frac{1}{r}\frac{d}{dr}\left(r\frac{d}{dr}F(r)\right)-\frac{n^2}{r^2}F(r)\left(A(r)-1\right)^2-\frac{1}{2}F(r)\left(F(r)^2-1\right)=0$$ There are two parameters $\xi\in\mathbb{R}^+$ and $n\in\mathbb{Z}$. There also are boundary conditions at $0$ and at $\infty$. $$F(0)=A(0)=0$$ $$F(\infty)=A(\infty)=1$$ These gave me the first puzzle: how to implement boundary conditions at infinity? I solved this by taking a finite number (in my case 5) instead of infinity. Do you think it is a good approach or is there a more elegant way to implement boundary conditions both at zero and infinity?
Next I want to solve these equations numerically. Here is the Mathematica code for $\xi=1$ and $n=1$

paramND = {
r D[1/r D[A[r], r], r] - \[Xi]^2 F[r]^2 (A[r] - 1) == 0,
1/r D[r D[F[r], r], r] - n^2/r^2 F [r] (A[r] - 1)^2 -
1/2 F[r] (F[r]^2 - 1) == 0,
A == 0.99, A[0.01] == 0,
F == 0.99, F[0.01] == 0} /. {\[Xi] -> 1, n -> 1};

numSol = NDSolve[paramND, {A[r], F[r]}, {r, 0.01, 5}] // Flatten;

paramPlot = {A[r], F[r]} //. numSol;

Plot[Evaluate@paramPlot, {r, 0.01, 5}, PlotRange -> {0, 1},
PlotStyle -> {{Black, Thick}, {Red, Thick}},
PlotLegends -> {"A(r)", "F(r)"}]


For $\xi=1$ and $n=1$ this works fine and I get a result that I expect: The problem is that the numerical solution only works if both parameters $\xi$ and $n$ are close to $1$. Changing these parameters to values different from $1$ produces somewhat awkwrd results:

$\xi=3$, $n=1$: $\xi=1$, $n=5$: As you can see, the last two plots don't even respect the boundary conditions! Why is this so? All solutions should look more or less like the first plot, which clearly isn't the case.
I suppose that I lack some fundamental knowledge about how to use NDSolve to produce results which I expect. Can you suggest a good strategy in order to obtain good solutions to these ODEs? I've also tried using the "Shooting" algorithm, which, however, didn't give any useful result either. Maybe you can point out to me what the exact reason is, why NDSolve fails here so I know what to cure. Up to now it was rather a try and error procedure.
So, to sum up, my two questions are:

1. How to implement boundary conditions at infinity?
2. How to choose a good strategy for NDSolve to obtain results one expects?
• I have tried to solve ur sys for the specific values you mentioned with maple and the output is accurate. – zhk Oct 13 '13 at 3:09

(2) Compactification

Update

Ok, my original post was rather compact. Let us start from the begining. You functions $A$ and $F$ map $[0,∞)$ into the reals. But you can think of the preimage as the result of another bijective function that maps say $[0,1)$ into $[0,∞)$. Since the domain of that second function can be easily modified to include $1$, the compound function will have its domain (at least formally) compactified. The question is: how is it done in Mathematica.

Let us use $r=ArcTanh(ρ)$ as the compactifying function.

Plot[ArcTanh[x], {x, 0, 1}, PlotRange -> All] Define $B(ρ)=A(r(ρ))$ and $G(ρ)=F(r(ρ))$. That is, the functions $B$ and $G$ are the functions $A$ and $F$ "in $ρ$-space".

The following code shows how to go from $r$ space into $ρ$ space and back:

{F[r], F'[r], F''[r]} /. F -> Function[{x}, G[Tanh[x]]] /. r -> ArcTanh[\[Rho]]
% /. G -> Function[{x}, F[ArcTanh[x]]] /. ρ -> Tanh[r] // FullSimplify[#, Assumptions -> r \[Element] Reals] &


{G[ρ], (1 - ρ^2) G'[ρ], -2 ρ (1 - ρ^2) G'[ρ] + (1 - ρ^2)^2 (G'')[ρ]}

{F[r], F'[r], (F'')[r]}

The above code basically computes $\frac{d F}{ d r}$ as $\frac{d G}{ d ρ} \frac{d ρ}{ d r}$ and a similar formula for $F''(r)$ in terms of $ρ$. You might want to check those results by hand to convince yourself of the replacements done in the differential equations:

eq = {r D[1/r D[A[r], r], r] - ξ^2 F[r]^2 (A[r] - 1) == 0,
1/r D[r D[F[r], r], r] - n^2/r^2 F[r] (A[r] - 1)^2 -
1/2 F[r] (F[r]^2 - 1) == 0} /. {
F -> Function[{x}, G[Tanh[x]]],
A -> Function[{x}, B[Tanh[x]]]
} /. r -> ArcTanh[ρ] // FullSimplify[#, Assumptions -> ρ \[Element] Reals] &


Similarly, the boundary conditions become

bc = {A == 0.99, A[0.01] == 0, F == 0.99, F[0.01] == 0} /. {
F -> Function[{x}, G[Tanh[x]]],
A -> Function[{x}, B[Tanh[x]]]
} // N


{B[0.999909] == 0.99, B[0.00999967] == 0., G[0.999909] == 0.99, G[0.00999967] == 0.}

The next modifications of your code are:

hNumSol = NDSolve[Flatten[{
eq /. {ξ -> 1, n -> 1},
bc
}], {B, G}, {ρ, Tanh[.001], Tanh}] // Flatten;
hParamPlot = {B[Tanh[r]], G[Tanh[r]]} //. hNumSol;


If you then plot

Plot[Evaluate@hParamPlot, {r, .01, 5}, PlotRange -> All,
PlotStyle -> {{Blue, Thick}, {Green, Thick}},
PlotLegends -> {"A(r)", "F(r)"}]


you will obtain exactly the same graph as in your question.

So far, we have reproduced your calculation but we have used $ρ$-space. In order to impose boundary conditions at infinity, impose them at $ρ=1$:

max = 1;
bc = {B[max] == 0.99, B[0.001] == 0., G[max] == 0.99, G[0.001] == 0};
hNumSol = NDSolve[Flatten[{eq /. {ξ -> 1, n -> 1}, bc}], {B, G},
{ρ, 0.001, max}] // Flatten;


The new plot misbehaves even though ξ=1 and n=1: Therefore, your problem is not due the inability to impose the boundary conditions at infinity. A solution worth exploring is to change the functions for which you are solving.

• I do not have time to work the solution, but now that you have the equations in ρ-space, find the $p$ and $q$ such that the $b$ defined by $B(ρ) = ρ + ρ^p (1-ρ)^q b(q)$ is a function that stays finite (and non-zero) at $ρ=0$ and $ρ=1$. Repeat for $g$ in $G(ρ) = ρ + ρ^s (1-ρ)^t g(ρ)$. Once you got the boundary behavior out of the way, get the differential equations in terms of $b$ and $g$. Try to solve those. – Hector Oct 13 '13 at 0:38
• Well, as @Hector has shown that your bvp has no issue for choosing finite value at the far field condition (infinity). The problem here is with your parameters values. Suppose, your bvp is in dimensionless form then you are free to choose whatever value you want for the parameters involved. But in your case, you not that lucky to do so. The best way will be to find a way, to select a specific domain for each parameter for which the bvp has a stable solution. – zhk Oct 13 '13 at 0:48
• Thanks @Hector for teaching me how to implement compactification in Mathematica, it's very useful, and answers entirely my first question. As you both Hector and MMM mentioned, the inability to find a correction numerical solution is caused by the numerical instability of the ODEs depending on the chosen parameters. Hector, you're suggesting to rewrite the ODEs in terms of $b(\rho)$ and $g(\rho)$, could you maybe comment on how you arrived at this specific ansatz? – Stan Oct 13 '13 at 20:53
• @MMM I didn' t quite understand what you meant by dimensionlessness in this contest, in my case all variables are dimensionless. – Stan Oct 13 '13 at 20:54
• Why are you using $0.99$ instead of $1$ in the boundary conditions? – becko Mar 21 '16 at 15:22

This is the most difficult of the nearly two dozen nonlinear ODE separatrix computations that I have encountered on Mathematica.SE. Nonetheless, it can be can be solved by a systematically refined search for initial conditions that maximize the range in r over which the ODE system can be integrated before clearly departing from the separatrix.

Series expansion at r == 0

To begin, we need the behavior of the dependent variables at very small r. The leading term for each is given by the OP in a comment above, and the next terms can be obtained by expanding the ODEs as power series at r == 0, equating the coefficients of each power of r to zero, and solving the resulting polynomial equations. (The next order term in the series expansions is desired for precision consistency with the high WorkingPrecision needed to integrate the ODEs.) These steps can be packaged as a simple function.

params[n_] := Module[{ser, sereq},
ser = Series[Unevaluated[{r D[1/r D[A[r], r], r] - ξ^2 F[r]^2 (A[r] - 1),
1/r D[r D[F[r], r], r] - n^2/r^2 F [r] (A[r] - 1)^2 - 1/2 F[r] (F[r]^2 - 1)}] /.
{A[r] -> Sum[ca[i] r^i, {i, 2, n + 4}],
F[r] -> Sum[cf[i] r^i, {i, n, n + 4}]}, {r, 0, n + 2}] // Normal;
sereq = Thread[DeleteCases[Flatten@(Most /@ (CoefficientList[#, r] & /@ ser)), 0] == 0];
{Sum[ca[i] r^i, {i, 2, n + 3}], Sum[cf[i] r^i, {i, n, n + 3}]} /.
Factor@Flatten@Solve[sereq, Flatten@{Table[ca[i], {i, 3, n + 3}],
Table[cf[i], {i, n + 1, n + 3}]}] /. {ca -> a, cf[n] -> f}]


where a and f are free parameters to be determined to obtain the separatrix. params[] yields for the first four values of n,

params /@ Range
(* {{a r^2 - 1/8 f^2 r^4 ξ^2, f r - 1/16 (1 + 4 a) f r^3},
{a r^2, f r^2 - 1/24 (1 + 16 a) f r^4},
{a r^2, f r^3 - 1/32 (1 + 36 a) f r^5},
{a r^2, f r^4 - 1/40 (1 + 64 a) f r^6}} *)


By observation the general term is

{a r^2 - If[n == 1, f^2 r^4 ξ^2/8, 0], f r^n (1 - (1 + 4 a n^2) r^2/(8 (n + 1)))}


Solution for n == 1, ξ == 1

The question provides a solution for n == 1, ξ == 1. Although the code in the question does not yield the plot shown for Mathematica 11.1.1, probably because NDSolve has changed over the past four years, the plot does provide a useful point of comparison. (A similar situation occurred for the answer by xzczd to question 110626, as described in its associated comments.)

Define the solution to the ODE system as

r0 = 1/10000; inf = 15; ξ = 1; n = 1;
numSol = ParametricNDSolveValue[
{r D[1/r D[A[r], r], r] - ξ^2 F[r]^2 (A[r] - 1) == 0,
1/r D[r D[F[r], r], r] - n^2/r^2 F [r] (A[r] - 1)^2 - 1/2 F[r] (F[r]^2 - 1) == 0,
A[r0] == a r0^2 - If[n == 1, f^2 r0^n ξ^2/8, 0],
F[r0] == f r0^n (1 - (1 + 4 a n^2) r0^2/(8 (n + 1))),
A'[r0] == 2 a r0 - If[n == 1, r0^3 f^2 ξ^2/2, 0],
F'[r0] == f  r0^(n - 1) (n - (2 + n) (1 + 4 a n^2) r0^2/(8 (1 + n))),
WhenEvent[A'[r] < 0 || F'[r] < 0 || A[r] > 1 || F[r] > 1, "StopIntegration"]},
{A, F}, {r, r0, inf}, {a, f}, WorkingPrecision -> 30, MaxSteps -> 20000];


Because the desired solutions monotonically approach A == 1 and F == 1 at large r, WhenEvent[] is used to terminate a calculation not doing so. The systematic search and refinement of the parameters a and f is accomplished by

dist[a_?NumericQ, f_?NumericQ] := (First@numSol[a, f])["Domain"][[1, 2]]


which determines how far in r an integration proceeds before stopping, and

Clear[paramFind];
paramFind[a1_, a2_, f1_, f2_, inc_, lp_, nm_: 1, ft_: 0] :=
Module[{mx = {{a1 + a2, f1 + f2, 0}/2}, ahw = (a2 - a1)/2, fhw = (f2 - f1)/2},
Do[mx = First[mx];
tab = Flatten[ParallelTable[{a, f, dist[a, f + ft]}, {a, mx[] - ahw,
mx[] + ahw, 2 ahw/inc}, {f, mx[] - fhw, mx[] + fhw, 2 fhw/inc}], 1];
ahw = 2 ahw/inc; fhw = 2 fhw/inc;
mx = TakeLargestBy[tab, Last, nm], {i, lp}]; mx]


It generates a Table of integration distances for the specified ranges of a and f, determines the maximum distance, and if requested repeats this process lp times, increasing the accuracy of the desired parameters by of order 1/inc with each iteration. It returns the final tab (as a side effect), and the largest nm distances and their parameters. For instance,

mmm = paramFind[1/10, 5/10, 3/10, 9/10, 50, 1, 20]


searches a large range of parameters, returning a, f, and the maximum distance integrated (as well as the next nineteen largest values in tab).

(* {63/250, 3/5, 3.61273127329305003526188484581} *)


The corresponding integration curves are given by

s = numSol @@ Most[First@mmm];
Plot[{(First@s)[r], (Last@s)[r]}, {r, r0, (First@s)["Domain"][[1, 2]]},
PlotRange -> {0, 1}, PlotStyle -> {{Black, Thick}, {Red, Thick}},
PlotLegends -> Placed[{"A(r)", "F(r)"}, {.9, .5}], AxesLabel -> {r, "A, F"}] which, not surprisingly, is not very close to the desired result. tab itself can be visualized with

ListDensityPlot[tab, ColorFunction -> "Rainbow", PlotRange -> All] The topology of dist in {a, f} space is a narrow ridge. It is much easier to find the peak on this ridge by aligning one axis of the search Table with the ridge. This is accomplished by

pltlim = MinMax[First /@ mmm];
ffit = Rationalize[LinearModelFit[Most /@ mmm, a, a, Weights -> Last /@ mmm] // Normal, 0]
Show[
ListPlot[Table[Style[Most[mmm[[i]]], ColorData[
"Rainbow", (Last[mmm[[i]]] - Last[mmm[[-1]]])/(Last[mmm[]] -
Last[mmm[[-1]]])]], {i, Length[mmm]}], DataRange -> pltlim],
Plot[ffit, Flatten@{a, pltlim}]]

(* 49771269/224512015 + (76623081 a)/51020704 *) Returning ffit to paramFind[] then give a much more useful result for tab.

mmm = paramFind[57/250, 67/250, -1/100, 1/100, 50, 1, 1, ffit];
mmm = (# + {0, ffit /. a -> First@#, 0}) & /@ mmm
ListDensityPlot[tab, ColorFunction -> "Rainbow", PlotRange -> All] which easily accommodates systematic refinement.

mmm = paramFind[57/250, 67/250, -1/100, 1/100, 10, 15, 1, ffit];
mmm = (# + {0, ffit /. a -> First@#, 0}) & /@ mmm

(* {{381469726561/1525878906250, 301261684232301644236719949/499387950864892578125000000,
15.0000000000000000000000000000}} *)


Plotting the integration results for this parameter set, using the code given above, yields The desired asymptotic values are reached at about r == 8. This entire calculation required less than five minutes of computer time. For completeness, here is another depiction of the ridge, this time using Plot3D. Although prettier, it required a good initial guess and a large fraction of an hour of computer time to obtain a less accurate result than that above.

Plot3D[dist[a, f], {a, 2499/10000, 2501/10000}, {f, 60316/100000, g0336/100000},
PlotPoints -> 51, MaxRecursion -> 3, PlotRange -> All, Mesh -> None]
TakeLargestBy[(% // InputForm)[[1, 1, 1]], Last, 1] (* {{0.25, 0.603262, 8.80276}} *)


Solution for n == 1, ξ == 3

Increasing ξ makes the computation more challenging, because the ridge becomes much narrower. Using the same code as above yields an initial plot of tab (after a bit of searching), from which ffit is obtained after an iteration as

(* 11817596/452896041 + (92677142 a)/127177533 *)


with the final result,

(* {{47443/40000, 34190793340749568069048951/38398800799897902000000000,
4.54199136850410415171442304534}} *) Perhaps, the integration could have been carried farther with WorkingPrecision ->45, but the curve seems good enough.

Incidentally, the corresponding optimum parameters for ξ == 3/2 (computed as a test case) are,

(* {{172416/390625, 5298661942732265830186933/7730236542592271015625000,
7.46251057225238295282438028696}} *)


Solution for n == 5, ξ == 1

Increasing n instead of ξ relative to the base case does not so much narrow the ridge as increase the needed values of WorkingPrecision -> 45, MaxSteps -> 30000 in numSol at n == 5 even to find the ridge. So, the computations become somewhat slower. In addition, "LocationMethod" -> "LinearInterpolation" is added to WhenEvent[] to eliminate warning messages. The final result, obtain using otherwise the same code as above, is

(* {{1/20, 26918753177298531983608529/5975307063483195750000000000,
15.0000000000000000000000000000000000000000000}} *) Parameters for n between 2 and 4 are,

(* {{1/8, 26137069653429078458702681/110681880730882350000000000,
15.0000000000000000000000000000}} *)
(* {{8333333333/100000000000, 13767675775706146512163493/
188286666204471770000000000, 15.0000000000000000000000000000}} *)
(* {{3124999999/50000000000, 374053654379840009082877269/
19447580662313216000000000000, 15.0000000000000000000000000000}} *)


It is curious that optimal values of a are simple rational numbers for ξ == 1, namely,

{1/4, 1/8, 1/12, 1/16, 1/20}


for n == Range. If this observation generally is true, solving ξ == 1 problems becomes much easier.

• I've come across a number of your answers like this while wandering around the site, and I'm always very impressed with them. Do you go into the detail of the separatrix theory somewhere? – KraZug Jul 10 '17 at 14:39
• I've run into similar equations in the study of topological defects in physics. (In fact, I suspect the OP was looking at either vortices in some kind of fluid or cosmic strings.) I solved my equations by going to MATLAB, which can solve boundary value problems via relaxation rather than via shooting. Mathematica can be cajoled into using finite element methods for ODEs, for which it uses relaxation methods instead; but it's a bit clunky. – Michael Seifert Jul 10 '17 at 15:45
• ... and now I remember that Mathematica's FEM algorithms can't handle non-linear ODEs. So that method wouldn't be useful here. – Michael Seifert Jul 10 '17 at 15:53
• @MichaelSeifert Thanks for your comments. It is too bad that Mathematica does not provide ODE solution by either homotopy or relaxation. You will find one separatrix problem solved by inverting the finite-difference approximation of an ODE here. What do you make of the observation in my last sentence? – bbgodfrey Jul 10 '17 at 19:21
• @KraZug Thanks for your kind comment. I have not written a general approach to solving separatrix problems. I have found: (1) if the ODE is autonomous, use that fact to reduce its order by one; (2) use Shooting, if good initial guesses are available; (3) as a last resort, search for parameters that maximize the distance of stable integration, as here and in 147207. – bbgodfrey Jul 10 '17 at 19:37
1. How to implement boundary conditions at infinity?

This can be troublesome in some cases, but for your equation the "taking a finite number instead of infinity" approach isn't bad, because the solutions tend to constants fast.

1. How to choose a good strategy for NDSolve to obtain results one expects?

### Pre-v12 Solution

Before v12 "Shooting" method is the only available method for solving nonlinear boundary value problem (BVP) of ordinary differential equation (ODE) in NDSolve, and your equation set turns out to be another example that "Shooting" method can't handle very well, as already shown in other answers.

So, what to do? Actually an unimpressive hint can be found in the tutorial Numerical Solution of Boundary Value Problems (BVP):

The advantage of the shooting method is that it takes advantage of the speed and adaptivity of methods for initial value problems. The disadvantage of the method is that it is not as robust as finite difference or collocation methods: some initial value problems with growing modes are inherently unstable even though the BVP itself may be quite well posed and stable.

Therefore finite difference method (FDM) may be suitable for your problem, let's try it out. I'll use pdetoae for the generation of difference equation:

zero = 0;
inf = 10;
eq = Map[r^2 # &, {r D[1/r D[A[r], r], r] - ξ^2 F[r]^2 (A[r] - 1) == 0,
1/r D[r D[F[r], r], r] - n^2/r^2 F[r] (A[r] - 1)^2 - 1/2 F[r] (F[r]^2 - 1) ==
0}, {2}] // Expand;
bc = {A[inf] == 1, A[zero] == 0, F[inf] == 1, F[zero] == 0};
domain = {zero, inf};
points = 25;
grid = Array[# &, points, domain];
difforder = 4;
(*Definition of pdetoae isn't included in this post,
ptoafunc = pdetoae[{A, F}[r], grid, difforder];
del = #[[2 ;; -2]] &;
ae[ξ_, n_] = del /@ ptoafunc@eq;

initialguess = 1;
valξ = 3; valn = 1;
solrule = FindRoot[{ae[valξ, valn], bc},
Flatten[#, 1] &@Table[{var[r], initialguess}, {var, {A, F}}, {r, grid}]];

{sol@A, sol@F} = ListInterpolation[#, domain] & /@ Partition[solrule[[All, -1]], points];

Plot[{sol[A][r], sol[F][r]}, {r, zero, inf}, PlotRange -> All,
PlotLabel -> Row@{ξ == valξ, ", ", n == valn}]  OK, FDM is indeed suitable for your problem.

### V12 Solution

Since v12, FiniteElement method of NDSolve supports nonlinear problem. This BVP is no exception.

eq, bc, etc. are the same as those in Pre-v12 Solution section:

{solfem@A, solfem@F} =
NDSolveValue[{eq, bc} /. {ξ -> valξ, n -> valn}, {A, F}, {r, zero, inf},
Method -> FiniteElement, InitialSeeding -> {A@r == 1, F@r == 1}];

Plot[{solfem[A][r], solfem[F][r]}, {r, zero, inf}, PlotRange -> All,
PlotLabel -> Row@{ξ == valξ, ", ", n == valn}]

• I was trying to find Nielsen-Olesen string numerical solution and this method is probably the simplest for a mathematica beginner like myself (just copy-paste the code with a few adjustments). Thank you! – Kosm Nov 25 '19 at 10:57
• @kosm If you're in v12, FiniteElement method may be a better choice. See my update. – xzczd Nov 25 '19 at 11:18