0
$\begingroup$

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

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, enter image description here, 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 enter image description here

Thanks again for your help!

$\endgroup$
7
  • $\begingroup$ 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. $\endgroup$
    – Nasser
    Aug 18, 2023 at 4:56
  • $\begingroup$ @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. $\endgroup$ Aug 18, 2023 at 5:46
  • $\begingroup$ 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. $\endgroup$
    – Nasser
    Aug 18, 2023 at 6:00
  • $\begingroup$ @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. $\endgroup$ Aug 18, 2023 at 7:13
  • $\begingroup$ 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. $\endgroup$
    – Nasser
    Aug 18, 2023 at 11:49

1 Answer 1

2
$\begingroup$

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"]

Mathematica graphics

Alternative to the above is

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

Mathematica graphics

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"]

Mathematica graphics

And this is for a=1

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

Mathematica graphics

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.

$\endgroup$

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.