Find the parameter in NDSolve giving the desired solution

Question

EDIT: I've edited many times because the system gives errors when the body is modified all at once. I've simplified the problem. The simplification consists of having reduced the number of parameters from 2 to 1.

I find myselef in choosing by hand the right parameter that let ParametricNDSolve return the desidered solution. I'm sure there is an easy alternative, but I can't find it anywhere. I thank anyone who can help.

I have a two ODEs system depending on 1 parameter $k$. Thanks to ParametricNDSolve I get the desidered family of solutions $h[k](x)$ and $v[k](x)$, as expected. Furthermore, I know (from books) that the solution $v[k]$ behaves at large $x$ like $v[k](x)\sim \frac{A_{k}}{x}+\frac{B_{k}}{x^2}$, so I can fit the tail of the function $v[k]$ using NonlinearModelFit for any given value of $k$. I'm eventually interested in finding the value of the parameter $k$ with the particular property that the fit gives $A_{k}=0$ and $B_{k}\neq 0$. How can this be done?

For now I found by hand that the value $k=7.94$ kinda works, but an automated method would be much better. Here is the code:

rh = 1
rhb = rh*(1 + 0.01)
q = 2
m = -2
mu = 31

system = {
  f'[x] == g[x], 
  g'[x] == -2 g[x]/x + 
    2*(q^2)*f[x]*v[x]^2/(x^2 - ((rh)^3)/x),
  v'[x] == a[x],
  a'[x] == -((2 x + ((rh)^3)/x^2)/(x^2 - ((rh)^3)/x))*a[x] - 2*a[x]/x - (q^2)*(f[x]/(x^2 - ((rh)^3)/x))^2*v[x] + (m/(x^2 -(rh)^3/x))*v[x],
  (*initial conditions below*)
  f[rhb] == 0 , 
  f[100000000] == mu, 
  v[rhb] == k, 
  a[rhb] == ((m)/(3*rh))*(k)}; 
sol = ParametricNDSolve[system, {f, v}, {x, rhb, 1000000}, {k}];

datascalar[k_] = 
  Table[{i, Evaluate[v[k][x] /. sol] /. x -> i}, {i, 200, 100000}];
nlmscalar[k_] = NonlinearModelFit[datascalar[k], A/x + B/x^2, {A, B}, x]

As can be read from the last line of the code, I can make a fit for every value of $k$! For example if I choose $k=7.94$ (found by hand) i get

nlmscalar[7.94]["BestFitParameters"]
{B -> 203.147, A -> -0.00218232}

and $A$ is "close" to zero, but it's not zero. I just need an algorithm that finds the value of $k$ for which the fit gives A -> 0.. This value is supposed to be near $7.94$. Said otherwise, I want the code to vary $k$ until the fit gives {B -> "whatever it is", A -> 0.}.

Answer 1

We can use NMinimize to fit parameters as follows

rh = 1;
rhb = rh*(1 + 0.01);
q = 2;
m = -2;
mu = 31;

system = {f'[x] == g[x], 
   g'[x] == -2 g[x]/x + 2*(q^2)*f[x]*v[x]^2/(x^2 - ((rh)^3)/x), 
   v'[x] == a[x], 
   a'[x] == -((2 x + ((rh)^3)/x^2)/(x^2 - ((rh)^3)/x))*a[x] - 
     2*a[x]/x - (q^2)*(f[x]/(x^2 - ((rh)^3)/x))^2*
      v[x] + (m/(x^2 - (rh)^3/x))*v[x],(*initial conditions below*)
   f[rhb] == 0, g[rhb] == e, v[rhb] == k, a[rhb] == ((m)/(3*rh))*(k)};
sol = ParametricNDSolve[system, {f, v}, {x, rhb, 10000}, {k, e}];
solm = 
 NMinimize[{Norm[
    Table[A/x + B/x^2 - Evaluate[v[k, e][x] /. sol] /. x -> i, {i, 
      200, 10000, 100}]], {k > 0, e > 0, A < 0, B > 0}}, {k, e, A, B}]

(*Out[]= {1.72976*10^-8, {k -> 0.310498, e -> 2.41152, A -> -0.12763, 
  B -> 0.796307}}*)

Note, this solution is differ from the hand made k=2,e=8 and also more precise. Visualization

p = {k, e} /. 
  solm[[2]]; Plot[{v[p[[1]], p[[2]]][x] /. sol, (A/x + B/x^2) /. 
   solm[[2]]}, {x, 200, 10000}, PlotStyle -> {Blue, Dashed}, 
 Frame -> True, PlotLegends -> {"NDSolve", "NMinimize"}]

Update 1. Let consider system of equations (6)-(7) from the paper Building an AdS/CFT superconductor. This system we transfer in to new one by introducing two functions $\Psi =u(y), \Phi =v(y)$ and mapping interval $(1,\infty)$ to $(0,1)$ with $y=1/r$, we have

sys={(2 u[y])/(1/y^2 - m y) + (v[y]^2 u[y])/(1/y^2 - m y)^2 + 
   2 y^3 Derivative[1][u][y] - 
   y^2 (2 y + (2/y + m y^2)/(1/y^2 - m y)) Derivative[1][u][y] + 
   y^4 (u^\[Prime]\[Prime])[y] == 
  0, -((2 v[y] u[y]^2)/(1/y^2 - m y)) + y^4 (v^\[Prime]\[Prime])[y] ==
   0} 

Note, in this model the Hawking temperature of the black hole is $T=\frac{3 m^{1/3}}{4 \pi L^{4/3}}$, where $L=1$ is the AdS radius, and we put $m=1$ as well. Boundary conditions for this system we put as follows

bc1 = {v[y] == 0, u[y] == 3 y0 u'[y]/2} /. 
  y -> y0; bc0 = { v[y] == mu - y rho, u[y] == A y + B y^2} /. y -> 0;

We are looking for some solutions with $A=0$. To solve this problem numerically we introduce new boundary conditions with two parameters $p_1=v'(1),p_2=u(1)$, then we normalize $v/p_1,u/p_2$, finally we have parametric function in the form

eps = $MachineEpsilon; y0 = 1 - 10^-6; m = 1; 
sol[p1_, p2_] := 
 NDSolve[{(2 u[y])/(1/y^2 - m y) + (
     p1^2 v[y]^2 u[y])/(1/y^2 - m y)^2 + 2 y^3 Derivative[1][u][y] - 
     y^2 (2 y + (2/y + m y^2)/(1/y^2 - m y)) Derivative[1][u][y] + 
     y^4 (u^\[Prime]\[Prime])[y] == 
    0, -((2  p2^2 v[y] u[y]^2)/(1/y^2 - m y)) + 
     y^4 (v^\[Prime]\[Prime])[y] == 0, v[y0] == 0, v'[y0] == 1, 
   u[y0] == 1, u'[y0] == 2/3/y0}, {u, v}, {y, eps, y0}]

We can plot function u'[0.] to visualize region where $A=0$

Plot3D[u'[eps] /. sol[p1, p2][[1]], {p1, 0.1, 5}, {p2, 0.1, 5}, 
 ColorFunction -> "Rainbow", Mesh -> None, PlotRange -> All, 
 Boxed -> False, PlotTheme -> "Marketing", PlotPoints -> 50, 
 AxesLabel -> Automatic]

Using ContourPlot we can retrieve data with A=0(see this post)

plot=ContourPlot[(Evaluate[u'[eps] /. sol[p1, p2][[1]]]) == 0, {p1, 0.1, 
   5}, {p2, .1, 5}] // Quiet

points = plot // 
  Cases[#, GraphicsComplex[points_, ___] :> points, Infinity] &

Update 2. In a case of original system we have

sys = {f'[x] == g[x], 
  g'[x] == -2 g[x]/x + 2*(q^2)*k^2 *f[x]*v[x]^2/(x^2 - ((rh)^3)/x), 
  v'[x] == a[x], 
  a'[x] == -((2 x + ((rh)^3)/x^2)/(x^2 - ((rh)^3)/x))*a[x] - 
    2*a[x]/x - (q^2)*g0^2 (f[x]/(x^2 - ((rh)^3)/x))^2*
     v[x] + (m/(x^2 - (rh)^3/x))*v[x]}; bc2 = {f[rhb] == 0, 
  g[rhb] == 1, v[rhb] == 1, a[rhb] == ((m)/(3*rh))*(1)};

sol2[p1_, p2_] := 
 NDSolve[{sys, bc2} /. {q -> 2, m -> -2, rh -> 1, rhb -> 11/10, 
    k -> p1, g0 -> p2}, {f, g, v, a}, {x, 11/10, 10}]

s2 = sol2[1, 1][[1]]

We map this system on the unit interval with using y=1/x as follows

eq = {-y^2 Derivative[1][f][y] == 
    g[y], -y^2 Derivative[1][g][y] == -2 y g[y] + (
     2 q^2  k^2 f[y] v[y]^2)/(
     1/y^2 - rh^3 y), -y^2 Derivative[1][v][y] == 
    a[y], -y^2 Derivative[1][a][
      y] == -2 y a[y] - ((2/y + rh^3 y^2) a[y])/(1/y^2 - rh^3 y) + (
     m v[y])/(1/y^2 - rh^3 y) - (
     q^2 g0^2 f[y]^2 v[y])/(1/y^2 - rh^3 y)^2};
bc1 = {f[rhb] == 0, g[rhb] == 1, v[rhb] == 1, 
   a[rhb] == ((m)/(3*rh))*(1)};
par = {q -> 2, m -> -2, rh -> 1, rhb -> 10/11};

eps = $MachineEpsilon; 
sol[p1_, p2_] := 
 NDSolve[{eq, bc1} /. {q -> 2, m -> -2, rh -> 1, rhb -> 10/11, 
    k -> p1, g0 -> p2}, {f, g, v, a}, {y, eps, 10/11}]

s1 = sol[1, 1][[1]]

We can compare two numerical solutions in one point, for example, at p1=1,p2=1,x=1/y=3 and have got same results

{f[x], g[x], v[x], a[x]} /. s2 /. x -> 3

(*Out[]= {1.51741, 0.624487, 0.384472, -0.166605}*)

{f[x], g[x], v[x], a[x]} /. s1 /. x -> 1/3

(*Out[]= {1.51741, 0.624487, 0.384472, -0.166605}*)

Then we compute function v'[0.] to visualize region where A=0

ContourPlot[(Evaluate[v'[eps] /. sol[p1, p2][[1]]]) == 0, {p1, 0.01, 
   1.5}, {p2, .01, 1.5}] // Quiet

This picture shows chaotic behavior of the numerical solution.