-1
$\begingroup$

Rotation speed of Milky Way looks like this until about 25 kpc

enter image description here

I know the following points to describe this curve

{{0.1, 0.566016}, {0.11, 0.575781}, {0.121, 0.5875}, {0.133, 
  0.600781}, {0.146, 0.620703}, {0.161, 0.653906}, {0.177, 
  0.703516}, {0.195, 0.767969}, {0.214, 0.834375}, {0.236, 
  0.889844}, {0.259, 0.929297}, {0.285, 0.954688}, {0.314, 
  0.969531}, {0.345, 0.977344}, {0.38, 0.980469}, {0.418, 
  0.979297}, {0.459, 0.975391}, {0.505, 0.96875}, {0.556, 
  0.960547}, {0.612, 0.95}, {0.673, 0.936719}, {0.74, 
  0.921094}, {0.814, 0.905078}, {0.895, 0.889844}, {0.985, 
  0.876953}, {1.083, 0.866016}, {1.192, 0.855859}, {1.311, 
  0.846875}, {1.442, 0.838672}, {1.586, 0.830859}, {1.745, 
  0.823828}, {1.919, 0.818359}, {2.111, 0.814453}, {2.323, 
  0.813281}, {2.555, 0.816016}, {2.81, 0.823047}, {3.091, 
  0.833594}, {3.4, 0.848438}, {3.74, 0.867188}, {4.114, 
  0.885156}, {4.526, 0.896484}, {4.979, 0.904688}, {5.476, 
  0.914453}, {6.024, 0.926563}, {6.626, 0.935547}, {7.289, 
  0.937891}, {8.018, 0.933594}, {8.82, 0.924609}, {9.702, 
  0.916016}, {10.672, 0.914844}, {11.739, 0.926172}, {12.913, 
  0.948438}, {14.204, 0.970703}, {15.625, 0.975391}, {17.187, 
  0.961719}, {18.906, 0.950391}, {20.797, 0.952734}, {22.876, 
  0.959375}, {25.164, 0.951953}}

We could Interpolate and get an Interpolating fonction vc[R] to solve the following problem :

Knowing e=0.995, I am trying now to get the density Rho[m]from the equation

vc[R]^2 == Integrate[(4 Pi Sqrt[1 - e^2] m^2 Rho[m])/ Sqrt[-e^2 m^2 + R^2], {m, 0, R}] This equation comes from Galactic Dynamics Binney Tremaine 2008. (2.132)

But solutions given elswhere Galactic rotation curve are working only with an analytical function vc[R], and I dont know why ??

Any idea ?

Thanks in advance

Improve this question
$\endgroup$
7
  • $\begingroup$ Your new dataset doesn't contain point R=0! $\endgroup$
    – Ulrich Neumann
    Commented Oct 26, 2023 at 16:51
  • $\begingroup$ Have a look at the transformed integral equation. The expression v[R]^2/(4Pi R^2) must be finite for all values of R. $\endgroup$
    – Ulrich Neumann
    Commented Oct 26, 2023 at 17:03
  • $\begingroup$ @ Ulrich Yes,you are right and I noticed these two points but I cant find something working. with this Interpolating function. .... $\endgroup$
    – Rémy Galli
    Commented Oct 26, 2023 at 17:27
  • $\begingroup$ 1. Please add the link to that "solutions given elswhere" to the body of your question. 2. My solution under your previous question fully applys to this new vc, which part are you having difficulty? (Of course you need to supplement the missing data at R==0. ) $\endgroup$
    – xzczd ♦
    Commented Oct 27, 2023 at 1:42
  • 2
    $\begingroup$ @RémyGalli In a case of our Galaxy the density has a singular part $m\delta(R)$, where $m$ is the black hole mass (see Nobel prize 2020). Therefore equation for vc[R] should be update. $\endgroup$
    – Alex Trounev
    Commented Oct 27, 2023 at 3:07

1 Answer 1

Reset to default
3
$\begingroup$

dataset should include point R=0. That's why I added point {0,0}.

With the transformed integralequation v[R]^2=Integrate[...,{s,0,e}] the solution evaluates.

e = 0.995;
 
    v = Interpolation[{{0, 0}, {0.1, 0.566016}, {0.11, 0.575781}, {0.121, 
     0.5875}, {0.133, 0.600781}, {0.146, 0.620703}, {0.161, 
     0.653906}, {0.177, 0.703516}, {0.195, 0.767969}, {0.214, 
     0.834375}, {0.236, 0.889844}, {0.259, 0.929297}, {0.285, 
     0.954688}, {0.314, 0.969531}, {0.345, 0.977344}, {0.38, 
     0.980469}, {0.418, 0.979297}, {0.459, 0.975391}, {0.505, 
     0.96875}, {0.556, 0.960547}, {0.612, 0.95}, {0.673, 
     0.936719}, {0.74, 0.921094}, {0.814, 0.905078}, {0.895, 
     0.889844}, {0.985, 0.876953}, {1.083, 0.866016}, {1.192, 
     0.855859}, {1.311, 0.846875}, {1.442, 0.838672}, {1.586, 
     0.830859}, {1.745, 0.823828}, {1.919, 0.818359}, {2.111, 
     0.814453}, {2.323, 0.813281}, {2.555, 0.816016}, {2.81, 
     0.823047}, {3.091, 0.833594}, {3.4, 0.848438}, {3.74, 
     0.867188}, {4.114, 0.885156}, {4.526, 0.896484}, {4.979, 
     0.904688}, {5.476, 0.914453}, {6.024, 0.926563}, {6.626, 
     0.935547}, {7.289, 0.937891}, {8.018, 0.933594}, {8.82, 
     0.924609}, {9.702, 0.916016}, {10.672, 0.914844}, {11.739, 
     0.926172}, {12.913, 0.948438}, {14.204, 0.970703}, {15.625, 
     0.975391}, {17.187, 0.961719}, {18.906, 0.950391}, {20.797, 
     0.952734}, {22.876, 0.959375}, {25.164, 0.951953}}
   ];

Rmin = 0; Rmax = 25;
ni = 25;

FEM approach:

Needs["NDSolve`FEM`"]

riGrid[n_] := 
Block[{dt, \[CapitalDelta]t0, \[CapitalDelta]t1, zw, sol}, 
dt = Map[\[CapitalDelta]t0 + (\[CapitalDelta]t1 - \[CapitalDelta]t0)(#/Rmax)^2 &, Subdivide[0, Rmax, n]];
zw = Accumulate[dt];
sol = Solve[{zw[[1]] == Rmin,zw[[-1]] == Rmax}, 
{\[CapitalDelta]t0,\[CapitalDelta]t1}][[1]];
zw /. sol // N
]



ri = riGrid[ni];(*problem adapted grid*)
(*{0., 0.00452489, 0.0226244, 0.0633484, 0.135747, 0.248869,0.411765,0.633484, 0.923077, 1.28959, 1.74208, 2.28959,2.94118, 3.70588, 4.59276, 5.61086, 6.76923, 8.07692, 9.54299, 11.1765,12.9864,14.9819, 17.1719, 19.5656, 22.1719, 25.}*)



\[Rho]i =Table[\[Rho][k], {k, 1, Length[ri]}];

netz = ToElementMesh[Map[{#} &, ri]];
\[Phi]i =Map[ElementMeshInterpolation[netz, #]&,IdentityMatrix[Length[ri]]];

Plot[Evaluate[Through[\[Phi]i[R]]], {R, Rmin, Rmax}];
int[R_?NumericQ] :=Block[{ },v[R]^2 - 4 Pi R^2 Sqrt[1 - e^2]/e^3
NIntegrate[
s^2/Sqrt[1 - s^2] Through[\[Phi]i[s R/e]], {s, 0, e}, 
Method -> {"FiniteElement", "LocalAdaptive"}[[1]]] . \[Rho]i]

zwint = Table[int[R], {R, ri}] /. 0. -> 0 // Quiet;

mini0 = NMinimize[zwint . zwint, \[Rho]i]
(*{-1.73472*10^-18, {\[Rho][1] -> 0.0352952, \[Rho][2] -> 0.206768, \[Rho][3] -> 0.170064, \[Rho][4] -> 0.178969, \[Rho][5] -> 0.171873,
\[Rho][6] -> 
0.161654, \[Rho][7] -> 0.138237, \[Rho][8] -> 0.10491, \[Rho][9] ->
0.0706109, \[Rho][10] -> 0.0444486, \[Rho][11] -> 
0.0275268, \[Rho][12] -> 0.0172977, \[Rho][13] -> 
0.0111636, \[Rho][14] -> 0.0074253, \[Rho][15] -> 
0.00507642, \[Rho][16] -> 0.00356313, \[Rho][17] -> 
0.00255434, \[Rho][18] -> 0.00187092, \[Rho][19] -> 
0.00139474, \[Rho][20] -> 0.0010538, \[Rho][21] -> 
0.000810122, \[Rho][22] -> 0.000628789, \[Rho][23] -> 
0.000496635, \[Rho][24] -> 0.000392479, \[Rho][25] -> 
0.000317462, \[Rho][26] -> 0.000254499}}*)


ListPlot[ Transpose[{ri, \[Rho]i} /. mini0[[2]]],GridLines -> {ri,None}]

enter image description here

Interpolation&Gaussian quadrature :

swGauss =GaussianQuadratureWeights[  15, 0, e] ; 



JJ = Sum[{1, 1/(1 + (R/2)^4)  }[[-1]] 
(v[R]^2  -4 Pi R^2   Sqrt[1 - e^2]/e^3 
Total@Map[#[[2]] #[[1]]^2 /Sqrt[1 - #[[1]]^2]Interpolation[Transpose[{Rationalize[ri], \[Rho]i}]][#[[1]]R/e] & ,swGauss] )^2, {R, ri}];

mini00 = NMinimize[{JJ  ,Map[# > 0 &, \[Rho]i] }, \[Rho]i] 
ListPlot[Rest@Transpose[{ri, \[Rho]i} /. mini00[[2]]],GridLines -> {ri, None}]

enter image description here

Improve this answer
$\endgroup$
7
  • $\begingroup$ I get negative values : could give the {ri,rhoi} please ? Thanks $\endgroup$
    – Rémy Galli
    Commented Oct 26, 2023 at 20:34
  • $\begingroup$ Orbital velocity is about 7,650 km/s in central part of Galaxy since there is the black hole Sagittarius A*. :) $\endgroup$
    – Alex Trounev
    Commented Oct 27, 2023 at 4:52
  • $\begingroup$ @RémyGalli Rerun my code with fresh kernel. ri and mini0 show you the values of these parameters. $\endgroup$
    – Ulrich Neumann
    Commented Oct 27, 2023 at 7:02
  • $\begingroup$ @Ulrich sorry Ulrich Could you display the ri and [Rho]i because i get a lot of negative values Thanks in advance.... $\endgroup$
    – Rémy Galli
    Commented Oct 27, 2023 at 8:49
  • $\begingroup$ @RémyGalli I modified my answer with displayed ri, rhoi ! Still I don't understand, why my code (with fresh kernel) doesn't evaluate in your Mathematica session. $\endgroup$
    – Ulrich Neumann
    Commented Oct 27, 2023 at 8:58

Your Answer

By clicking “Post Your Answer”, you agree to our terms of service and acknowledge you have read our privacy policy.

Not the answer you're looking for? Browse other questions tagged or ask your own question.