# Problem involving a system of nonlinear coupled ODE's with adjustable boundary

The problem is to solve the following system ODEs:

I. $$\ \ \ \ \ \ \ \dfrac{4}{r}[1+a(r)]\left(\dfrac{dH}{dr}\right)^2+\dfrac{dG}{dr}=0$$

II. $$\ \ \ \ \dfrac{1}{r}\left(\dfrac{dA}{dr}\right)+f(r)+k^2G(r)=0$$

III. $$\ \ \ \ \left(\dfrac{dG}{dr}\right)^2+4k^2\left(G(r)+\dfrac{f(r)}{k^2}\right)^2\left(\dfrac{dH}{dr}\right)^2-(1.6)H(r)\left(\dfrac{dH}{dr}\right)^2=0$$

IV. $$\ \ \ \ \dfrac{d^2f}{dr^2}+\dfrac{1}{r}\dfrac{df}{dr}+\dfrac{1}{r}\dfrac{dA}{dr}-4k^2f(r)\left(\dfrac{dH}{dr}\right)^2=0$$

where $$k>0$$. The boundary conditions are

$$H(0)=1, \ \ A(0)=0, \ \ G(L)=0, \ \ f'(0)=0, \ \ f(L)=0 \ \ [\text{and} \ \ f'(L)=0]$$

where $$L$$ is the boundary, which must be adjusted in order to satisfy the boundary conditions above for each value of parameter $$k$$ chosen. Note that we have 6 (six) boundary conditions which must be satisfied, but we need only 5 (five). This sixth condition needs to be guaranteed in some way and this is my difficulty. The condition $$f'(L)=0$$ is necessary due to be a physical problem. Further, this physical problem requires also:

$$H(L)=H'(L)=0, \ \ A'(L)=0$$.

Thus, assuming these conditions mentioned, a brief analysis of the ODEs allows conclude that:

The ODE I ensures $$H'(0)=0$$ and $$G'(L)=0$$ wiht 1 (one) free parameter $$a(L)=const.$$

The ODE II ensures $$A'(0)=0$$ and 2 (two) free parameters $$G(0)=const.$$ and $$f(0)=const.$$

The ODE III ensures $$G'(0)=0$$.

It was possible to find a numerical solution using the shooting method as follow, but already advance that in this case the boundary condition A'[L]=0 is not satisfied as I want. It follow that,

ode1 = 4 (h'[r] ^2) (a[r] + 1)/r == -G'[r];
ode2 = a'[r]/r == -(k^2) (G[r]) - f[r];
ode3 = (G'[r])^2 + 4 (k^2) (h'[r]^2) (G[r] + f[r]/k^2)^2 == (8/5) h[r] (h'[r]^2);
ode4 = f''[r] + f'[r]/r + (a'[r]/r) == 4 (k^2) (f[r]) ((h'[r])^2);


Simplifying ode1 and ode3 can avoid problems with complex numbers. This is done using Solve[{ode1, ode3}, {h'[r], G'[r]}] // FullSimplify. The new set of equations is reduced to form

ode13 = {Derivative[1][h][r] == -((r Sqrt[2k^2 h[r] - 5 (f[r] + k^2 G[r])^2])/(2 Sqrt[5]k(1+a[r]))),Derivative[1][G][r] == (r(-2 k^2 h[r] + 5 ( f[r] + k^2 W[r])^2))/(5k^2 (1 + a[r]))};
ode2 = a'[r]/r == -(k^2) (G[r]) - f[r];
ode4 = f''[r] + f'[r]/r + (a'[r]/r) == 4 (k^2) (f[r]) ((h'[r])^2);


Applying the algorithm for k=1, as example, we have:

eps = 10^-10; L = 2.667; k = 1;
ode13 = {Derivative[1][h][r] == -((r Sqrt[2k^2 h[r] - 5 (f[r] + k^2 G[r])^2])/(2 Sqrt[5]k(1+a[r]))),
Derivative[1][G][r] == (r(-2 k^2 h[r] + 5 ( f[r] + k^2 G[r])^2))/(5k^2 (1 + a[r]))};
ode2 = a'[r]/r == -(k^2) (G[r]) - f[r];
ode4 = f''[r] + f'[r]/r + (a'[r]/r) == 4 (k^2) (f[r]) ((h'[r])^2);
bcs = {h[eps] == 1, a[eps] == 0, G[L] == 0, f'[eps] == 0, f[L] == 0};
sol = NDSolve[{ode13, ode2, ode4, bcs}, {a, h, f, G}, {r, eps, L},
Method -> {"Shooting","StartingInitialConditions" -> {f[0] == -0.340181, G[0] == 0.736905}}];
Plot[Evaluate[Re[{h[r], h'[r]}] /. sol], {r, eps, L}, PlotLegends -> {"h", "h'"}, PlotRange -> All]
Plot[Evaluate[Re[{a[r], a'[r]}] /. sol], {r, eps, L}, PlotLegends -> {"a", "a'"}, PlotRange -> All]
Plot[Evaluate[Re[{f[r], f'[r]}] /. sol], {r, eps, L}, PlotLegends -> {"f", "f'"}, PlotRange -> All]
Plot[Evaluate[Re[{G[r], G'[r]}] /. sol], {r, eps, L}, PlotLegends -> {"G", "G'"}, PlotRange -> All]


In this case, the values of f[0] and G[0] were adjusted together with the value of L in order to satisfy the boundary conditions. With the exception of f'[L]=0, all the other boundary conditions has been satisfied.

Then, main question is:

Is it possible to obtain a numerical solution by assuring the boundary condition f'[L]=0 in some way together with the others already reached in the solution shown above?

Remark: Note that in the above solution, it was only the condition f[L]=0 was used and not f'[L]=0. I tried to substitute one for the other but not works, that is, it does not satisfy the problem.

• It is necessary to give a link to the original topic and all solutions mathematica.stackexchange.com/questions/180269/… – Alex Trounev Oct 24 '18 at 6:45
• For each additional boundary condition, a parameter is required. For example, you can pick k for which f'[L]=0 – Alex Trounev Oct 24 '18 at 6:57
• How do you obtain h[L] == 0 mathematically from the equations? – bbgodfrey Oct 27 '18 at 1:19
• @acsantos With Normal@Series[eq, {r, r0, 0}] /. {f[r0] -> 0, G[r0] -> 0} with r0 an arbitrary value of r > 0 yields a'[r0] -> 0 but not h[r0] -> 0, unless r0 is equal to the particular value of L given in your question. In other words, you appear to be choosing L to cause h[L] to vanish. However, there appears to be a branch point in the solution at that value of L, and integrating past is does not seem possible. In any case, I have spent some time attempting to find a way to obtain f'r0] -> 0 but without success. Probably, it cannot be done. – bbgodfrey Oct 28 '18 at 2:36
• @acsantos By the way, Solve[{ode1, ode3}, {h'[r], G'[r]}] has three solutions. Are you sure that you have chosen the correct one? – bbgodfrey Oct 28 '18 at 2:38

In principle, it should be possible to use bisection search to find a suitable parameter L. Define

ClearAll[F, L];
F[L_?NumericQ] :=
f'[L] /. NDSolve[{ode13, ode2, ode4, bcs}, {a, h, f, G}, {r, eps, L},
Method -> {"Shooting",
"StartingInitialConditions" -> {f[0] == -0.340181,
G[0] == 0.736905}}][[1]]


and run

FindRoot[F[L], {L, L0, L1}, Method -> "Brent"]


where L0 and L1 have to bracket a solution, i.e., F[L0] and F[L1] need to have opposite signs. But in your particular case, this seems to be impossible: NDSolve produces bogus for values greater than L = 2.667 and for values below 2.667, we have

LList = Subdivide[0.0000001, 2.667, 100];
data = F /@ LList;
ListLinePlot[Transpose[{LList, data}], AxesLabel -> {"L", "f'[L]"},
PlotRange -> {0, Max[data] 1.1}]


This sort of implies that the only solution to your problem is L=0.

• Dear @HenrikShumacher, Could I take the derivative of any of the equations, for example Eq. II, in order to introduce the boundary condition f'[L]=0` and, thus, try to solve the problem? Is it acceptable to do this? – sotnasac Oct 25 '18 at 1:01