# A problem of Mathematica function with different parameters in the same model

I'm simulating a graph of a function, and have made a preliminary attempt with mathematica, my test code is as follows::

Clear["*"]
M = 1/2; a = 1; s = 1; intc = 5; L = 0; b = 1;
f[r_] = 1 - L - (2 M)/Sqrt[r^2 + a^2];
rs[r_] = Integrate[1/f[r], r] - intc;
v[r_] = (1 - L - (2 M)/Sqrt[
r^2 + a^2])*((b*(b + 1))/(r^2 + a^2) + (1 -
s)*((1/(r^2 + a^2)^(5/2)) (2*M*r^2 +
a^2 (-2 M - (-1 + L) \[Sqrt](r^2 + a^2)))));
V[r_] = v@InverseFunction[rs][r];
Plot[V[r], {r, -200, 20}, PlotRange -> All]


My results show that when $$a=1$$, I achieve the second one of Figure 1, , but when I test $$a=0.99$$, I failed. In principle, the same model only has different parameters. I don't know whether it is a problem of function setting. I would appreciate it if you could provide some help.

Another method of the test codes as follows,

Clear["*"]
a = 0.99;
f[r_] = 1 - 1/Sqrt[r^2 + a^2];
ff[r_] = Integrate[1/f[r], r];
fid = OpenWrite["/Users/cc.csv"];
data = For[n = 1, n <= 2000, n++, rs = -250 + 0.25*n;
M = 1/2; a = 0.99; s = 1; intc = 5; L = 0; b = 1;
x = r /. NSolve[-intc + ff[r] - rs == 0, r, Reals][[1]];
v = (1 -
L - (2 M)/
Sqrt[x^2 + a^2])*((b*(b + 1))/(x^2 + a^2) + (1 -
s)*((1/(x^2 + a^2)^(5/2)) (2*M*x^2 +
a^2 (-2 M - (-1 + L) \[Sqrt](x^2 + a^2)))));
Export[fid, List[{rs, v}], "CSV", "FieldSeparators" -> ","]];
Close[fid];

data = Import["/Users/cc.csv"];
ListPlot[{data[[1 ;; All]], data[[All ;; 1]]}, PlotRange -> All,
Joined -> True]


and the result is

Thanks again for your help!

• I can't read your code because you are using l = 1; which makes it hard to read. Do not use lower case L for variable as it looks like 1. If you use different letter I can try to help. Aug 18, 2023 at 4:56
• @Nasser Thank you very much for your suggestion. I have replaced the lowercase $l$ with $b$, and added a picture of the original code, hoping to give you a good read. Aug 18, 2023 at 5:46
• btw, you have to watch having integrate command taking r as argument and then integrating w.r.t. to that argument. If r is numerical, this will fail. The code you wrote can be re-written to avoid this problem coming into play. But I assume this is working for you so I did not change it. But I would not write it this way. Aug 18, 2023 at 6:00
• @Nasser Thanks for the possible solution to the problem. Physically, the parameter $a$ represents different physical backgrounds. The parameter $a$ can be divided into two parts, the cases of $a>=1$ and $0<a<1$. Currently, the test code I wrote is only applicable to $a>=1$. I don't know the reason for this. Before this, I used Matlab to complete the case of $0<a<1$ (the idea is to first discretize $rs$, replace the numerical result of $rs$ integral with an analytical expression, then solve inversely $r$, that is, $r[rs]$. Finally, substitute $r[rs]$ into $V[r]$). But in Mathematica, I failed. Aug 18, 2023 at 7:13
• If you have Matlab code that works for $a<1$, then why not just translate that code as is to Mathematica? This way you know it will work the same way. It should be easy to do the translation assuming your did not use any special functions in Matlab which do not exist in Mathematica or something like this. I can only use the code you provided and that gives problems with $a<1$, with Mathematica, it says the result is complex. Aug 18, 2023 at 11:49

When using Plot, it uses real numbers. For a=0.99 the code you showed give complex numbers. Also it seems there are division by zero encountered. One option is to generate the data your self and use ListPlot instead and use Abs.

I used division of 1.0. Notice that Plot uses adaptive divisions instead. But this way, we have avoided that problem.

Clear["*"]
M = 1/2; a = 99/100; s = 1; intc = 5; L = 0; b = 1;
f[r_] = 1 - L - (2 M)/Sqrt[r^2 + a^2];
rs[r_] = Integrate[1/f[r], r] - intc
v[r_] = (1 -
L - (2 M)/
Sqrt[r^2 + a^2])*((b*(b + 1))/(r^2 + a^2) + (1 -
s)*((1/(r^2 + a^2)^(5/2)) (2*M*r^2 +
a^2 (-2 M - (-1 + L) \[Sqrt](r^2 + a^2)))));
V[r_?NumericQ] = v@InverseFunction[rs][r];


And now

data = Table[V[r], {r, -200, 20, .1}];
ListLinePlot[Abs@data, PlotRange -> All, PlotLabel -> "Using a=0.99"]


Alternative to the above is

Plot[Abs@V[N@r], {r, -200, 20}, PlotRange -> All, PlotLabel -> "Using a=0.99"]


This is using a=0.9 (need to change a at top)

Plot[Abs@V[N@r], {r, -200, 20}, PlotRange -> All, PlotLabel -> "Using a=0.9"]


And this is for a=1

Plot[Abs@V[N@r], {r, -200, 20}, PlotRange -> All, PlotLabel -> "Using a=1"]


Note that the above plot for a=0.99` does not match the one shown in the question. I do not know why, but that is what the code gives.