# 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[5] == 0.99, A[0.01] == 0,
F[5] == 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. – MMM 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[5] == 0.99, A[0.01] == 0, F[5] == 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[5]}] // 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. – MMM 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 at 15:22