0
$\begingroup$

I am trying to obtain a plot between two variables related through an equation as follows :

solk[logn_?NumericQ] := 
 logk /. FindRoot[(1.9894350000000002`*^-10 10^
        logn (-1 + 10^logn) Cos[
        200 (10^-logk)*ArcSinh[5*10^2 (10^logk)]])/
     Sqrt[1 + (25*10^2 10^(2 logk))] == -(3/
       100000), {logk, -5, -200, 0}]

list = Table[{logn, Quiet@solk[logn]}, {logn, 1/100, 5, 1/100}];

p1 = ListLinePlot[list, PlotRange -> {{0, 5}, Full}]

As can be seen I am taking an initial guess of logk = -5 which gives me some plot. Now, I change the initial guess from -1 to -120 (yes, literally, one step at a time) and try to observe the behaviour. The behaviour is :

  1. for initial guess as <code>-3</code>
  2. for initial guess as <code>-5</code>
  3. for initial guess as <code>-5</code> with <code>DampingFactor->0.001</code>
  4. for initial guess as <code>-10</code> or less

Image 1 : for initial guess as -3

Image 2 : for initial guess as -5

Image 3 : for initial guess as -5 with DampingFactor->0.001

Image 4 : for initial guess as -10 or less

To be noted that for initial guess less than or equal to -10 (say the guess is x) the plot has exact same behaviour : the plot has flat logk value of x until logn is around 3 and then jumps immediately to the maximum value of logk range provided to FindRoot.

Can anyone please help me obtain the correct result here? Do I need to use some other algorithm here? From this plot I need to be able to say things like : "When logk is around -120 (or so) then the value of logn is something.

$\endgroup$
2
  • $\begingroup$ You need to link to the *.png itself, and not just the Imgur page for it. I have done it for you this time. $\endgroup$ Commented May 8, 2020 at 2:14
  • $\begingroup$ @J.M. Understood! Thank you very much. $\endgroup$
    – Nitin
    Commented May 8, 2020 at 2:19

1 Answer 1

2
$\begingroup$

ORIGINAL:

This isn't an answer, but I needed to add pictures and maybe it will help a bit anyways.

The function that you're working with is ghastly, and I don't know remotely enough to figure out how to solve for zeroes that are in the messy part, though solve for the most negative zero shouldn't be too difficult.

I usually like to start by plotting functions to get an idea of what they might look like, although I didn't think to do it when I answered your other question. Let's pretend logn = 5, and plot for logk from -6 to 2. Here the logk value is on the x-axis, and we want to know when the total equation is equal to 0 on the y-axis:

Plot[(1.9894350000000002`*^-10 10^logn (-1 + 10^logn) Cos[
       200 (10^-logk)*ArcSinh[5*10^2 (10^logk)]])/
    Sqrt[1 + (25*10^2 10^(2 logk))] + (3/100000) /. {logn -> 5},
 {logk, -6, 2},
 PlotPoints -> 100,
 MaxRecursion -> 10]

Plot of function with logn = 5.

How do you decide which zero to grab? I don't even know how to count how many zeroes there are, but there are a lot. If you just want the most negative one, that might not be too bad since there's a good amount of space there, but if the initial guess happens to be anywhere in the range of -4 to 2, or if the algorithm ends up searching somewhere inside that region, I suspect the outcome will be completely unstable.

We can take a look at the plot in 3D to get a better understanding of what's happening:

Plot3D[
 (1.9894350000000002`*^-10 10^logn (-1 + 10^logn) Cos[
      200 (10^-logk)*ArcSinh[5*10^2 (10^logk)]])/
   Sqrt[1 + (25*10^2 10^(2 logk))] + (3/100000),
 {logn, 0, 20},
 {logk, -6, 0},
 AxesLabel -> {"logn", "logk", "Value"},
 PlotPoints -> 100,
 PlotRange -> {-0.1, 0.1}
 ]

3DPlot of function for logn and logk.

I didn't bother asking it to plot more points, but in that messy region is where we have millions of zeroes.

This leads me to wonder what would happen if I took a look again at that first plot, but with a few different values of logn. The following is just cross-sections of the 3D graph:

Plot[Evaluate@
  Table[(1.9894350000000002`*^-10 10^logn (-1 + 10^logn) Cos[
        200 (10^-logk)*ArcSinh[5*10^2 (10^logk)]])/
     Sqrt[1 + (25*10^2 10^(2 logk))] + (3/100000), {logn, {5, 3, 
     0}}], {logk, -6, 0}, AxesLabel -> {"logk", "Value"}, 
 PlotPoints -> 100, PlotRange -> {-0.0004, 0.0004},
 PlotLegends -> {"logn = 5", "logn = 3", "logn = 0"}]

Plot of logk for various values of logn.

So it seems like when logn = 0, the equation doesn't have any roots. For logn = 3, the equation has tons of roots but the amplitude is fairly small. When logn = 5, the amplitude is gigantic and there are also tons of closely spaced roots.

If you're okay with just find any one root, you could always try limiting the search so that it only finds that one fairly nicely spaced root around -5. If you were hoping that the function only crossed zero one, then I think you're out of luck. It's kind of like solving for where Sin[x] == 0. It has infinitely many roots, so there is no unique solution. But if you just want to find one possible root, or to see if there's a root between -1 and 1, or something, those are well-posed tasks.

EDIT:

In case you're okay with just finding the most negative root, I'll put that information here. If you have some other criteria for selecting the root, that will be more difficult.

If we plot a few values of logn:

ns = {3, 2.59, 2.5, 2};
Plot[Evaluate@
  Table[(1.9894350000000002`*^-10 10^logn (-1 + 10^logn) Cos[
        200 (10^-logk)*ArcSinh[5*10^2 (10^logk)]])/
     Sqrt[1 + (25*10^2 10^(2 logk))] + (3/100000), {logn, 
    ns}], {logk, -10, 0}, AxesLabel -> {"logk", "Value"}, 
 PlotPoints -> 100, PlotRange -> {-0.00004, 0.00004}, 
 PlotLegends -> ("logn = " <> ToString[#] & /@ ns)]

enter image description here

We can see that it's not until logn equals about 2.59 that we even have any roots. We can also see that the first root as we increase logn begins around -5.7, and then rapidly converges to -4.71683. I've modified solk here to only search in the region of -10 to -4.6. If it goes beyond -4.6, it will get stuck in the nasty region.

solk[logn_?NumericQ] := 
 logk /. FindRoot[(1.9894350000000002`*^-10 10^
        logn (-1 + 10^logn) Cos[
        200 (10^-logk)*ArcSinh[5*10^2 (10^logk)]])/
     Sqrt[1 + (25*10^2 10^(2 logk))] == -(3/
       100000), {logk, -4.7, -10, -4.6}]
list2 = Table[{logn, solk[logn]}, {logn, 2.59, 5, 0.01}];
ListLinePlot[list2, AxesLabel -> {"logn", "logk"}, PlotRange -> Full]

enter image description here

$\endgroup$
7
  • $\begingroup$ This is an absolute gem for me. I am looking for the behaviour in the region when logk is between -130 & -110. To avoid the algorithm going into the messy region I could take DampingFactor->0.1 or so. I tried what you have provided to Plot for several values of logn & determined that in the region of logk I am interested in, the value of logn is about 2.5898936 for my equation to be zero upto 10 digits of precision. I can atleast make some statement. The only problem being : there are many values of logk for which logn is 2.5898936. Can you tell me a good way to present this? $\endgroup$
    – Nitin
    Commented May 8, 2020 at 16:11
  • $\begingroup$ @Nitin Hmmm, well there's a very small range of logn which give a logk in the range of -130 to -110 as the curve is very quickly going to $-\infty$ in that region. You'll definitely run into precision issues. The main issue is the $1.9894350000000002 \times 10^{-10}$ which is only machine precision. How well do you know that number? If that can be turned into an exact number, then there's a chance to solve it. By assuming it was an exact number, I found the difference between logn when logk = -130 and logn when logk = -110 to be on the order of $3 \cdot 10^{-213}$. $\endgroup$
    – MassDefect
    Commented May 8, 2020 at 17:32
  • $\begingroup$ Sorry, I don't understand this comment completely. Do you mean to ask how precisely I know that logn = 2.5898936? The value of the expression is 0 when the plot scale is of the order of 1*10^-10 or so. So I am guessing that upto 10 digits of precision logn = 2.5898936. I do agree with your last statement about the difference between logn. I am actually interested at logk = -122. Originally, I was hoping to get a peak around at logk = -122 for some logn (basically a refinement of plot 2 above but without the irregular part towards the left) . But that seems impossible now. $\endgroup$
    – Nitin
    Commented May 8, 2020 at 18:11
  • $\begingroup$ ...and if that is indeed impossible then I would have like to make the statement other way round : For logk = -122 value of logn = 2.5898936 but the problem is I would have preferred that logn be something different for other logk values but as you point out the difference is too less!! Now, I don't know what kind of statement I can make here and how to properly represent it through a plot. $\endgroup$
    – Nitin
    Commented May 8, 2020 at 18:15
  • $\begingroup$ @Nitin In your original equation to FindRoot you have a number like $1.989435 \times 10^{-10}$. All of the other numbers in the equation are exact, but this one number has a decimal place. This decimal place means that it is a machine precision number with at most 16 decimal places. This is far too few decimal places to perform the calculation. Do you know any more decimal places of that number? Is it okay to pretend it's the exact number $19894350000000002 \times 10^{-26}$? It might seem weird at first, but exact numbers allow Mathematica to calculate values much more carefully. $\endgroup$
    – MassDefect
    Commented May 8, 2020 at 18:17

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.