Message Boards Message Boards

GROUPS:

Solving an equation that involves truncated series

Posted 14 days ago
257 Views
|
9 Replies
|
10 Total Likes
|

Hi all, I don't know how to handle this problem. I'll report the full code below. In short, I have an equation (at the bottom) in which delta(k) doesn't seem to give problems (I plotted it, and it works), while phi(k) contains a truncated series evaluated numerically and which is itself a function of k, so I can only plot it using DiscretePlot[]. There is also n*Pi because I have to solve the equation for some values of n, from 0 up to a certain integer. So, how can I find a solution to this equation? I just inserted Reduce, for example. Thank you, guys.

P.s.: I noticed that phi[k] can be plotted continuously, defining it with "=" instead of ":=". So I'll try this way. Sorry, I don't know how to eliminate this question.

m = 10;
ll = 1;
mro = 1;
mpi = 1;
g = 5.5;
vl := Table[{i, j, l}, {i, -m, m}, {j, -m, m}, {l, -m, m}];
vll := Flatten[vl, 2];
s[v_, k_] := 1/(v^2 - (k*ll/(2*Pi))^2);
ss[k_] := Total[Map[s[Norm[#], k] &, vll]];
phi[k_] := ArcTan[-2*(Pi)^2*k*ll/(2*Pi*ss[k])];
omsq[k_] := ((mro)^2 + k)*4;
gamma[k_] := g^2*k^3/(6*Pi*omsq[k]);
h1[k_] := 
  g^2*k^3*Log[(Sqrt[omsq[k]] + 2*k)/(2*mpi)]/(3*(Pi)^2*
      Sqrt[omsq[k]]);
h2[k_] := 
  g^2*k^2*(1 + (1 + 2*(mpi)^2/omsq[k])*Sqrt[omsq[k]]*
       Log[(Sqrt[omsq[k]] + 2*k)/(2*mpi)]/k)/(6*(Pi)^2*
      Sqrt[omsq[k]]);
delta[k_] := 
  ArcCot[((mro)^2 - omsq[k] - 
      h1[mro] - (omsq[k] - (mro)^2)*h2[mro]/(2*mro) + h1[k])/(omsq[k]*
      gamma[k])];
n = 0;
Reduce[delta[k] + phi[k] == n*Pi && k >= 0, k, Reals]
POSTED BY: Davide Alfano
9 Replies

Transcendental equation can't be solve symbolically, Reduce can't solve. We can solve with FindRoot:

n = 0;
Plot[Evaluate[delta[k] + phi[k] == n*Pi], {k, 0, 15}](*We see on the Plot a good starting points to solve roots*)

FindRoot[delta[k] + phi[k] == n*Pi, {k, 1/10}]
FindRoot[delta[k] + phi[k] == n*Pi, {k, 1/2}]
FindRoot[delta[k] + phi[k] == n*Pi, {k, 11}]
FindRoot[delta[k] + phi[k] == n*Pi, {k, 13}]

Or with NSolve:

 NSolve[delta[k] + phi[k] == n*Pi && 0 <= k < 15, k, Reals]

 (*{{k -> 0.}, {k -> 0.415907}, {k -> 3.97402}, {k -> 6.40243}, {k -> 
    9.0614}, {k -> 10.9818}, {k -> 12.6251}, {k -> 14.2583}}*)

Regards M.I.

POSTED BY: Mariusz Iwaniuk
Posted 12 days ago

Thank you for your reply. I made a (in any case irrelevant) mistake in the definition of omsq[k]. I report the correct code:

m = 10;
ll = 1;
mro = 1;
mpi = 1;
g = 5.5;
vl := Table[{i, j, l}, {i, -m, m}, {j, -m, m}, {l, -m, m}];
vll := Flatten[vl, 2];
s[v_, k_] := 1/(v^2 - (k*ll/(2*Pi))^2);
ss[k_] := Total[Map[s[Norm[#], k] &, vll]]
phi[k_] = ArcTan[-2*(Pi)^2*k*ll/(2*Pi*ss[k])];
omsq[k_] := ((mro)^2 + k^2)*4;
gamma[k_] := g^2*k^3/(6*Pi*omsq[k]);
h1[k_] := 
  g^2*k^3*Log[(Sqrt[omsq[k]] + 2*k)/(2*mpi)]/(3*(Pi)^2*
      Sqrt[omsq[k]]);
h2[k_] := 
  g^2*k^2*(1 + (1 + 2*(mpi)^2/omsq[k])*Sqrt[omsq[k]]*
       Log[(Sqrt[omsq[k]] + 2*k)/(2*mpi)]/k)/(6*(Pi)^2*
      Sqrt[omsq[k]]);
a0 := h1[mro] - mro*h2[mro]/2 + g^2*(mpi)^2/(6*(Pi)^2)
a[k_] := 
  h1[mro] + (omsq[k] - (mro)^2)*h2[mro]/(2*mro) - h1[k] + 
   I*Sqrt[omsq[k]]*gamma[k];
fpi[k_] := ((mro)^2 - a0)/((mro)^2 - omsq[k] - a[k])
delta[k_] := 
  ArcCot[((mro)^2 - omsq[k] - 
      h1[mro] - (omsq[k] - (mro)^2)*h2[mro]/(2*mro) + h1[k])/(omsq[k]*
      gamma[k])];

However, I'm looking for an instruction that automatically finds a solution without the need to set a starting point, so NSolve seems to be my choice. And in fact, I already tried it before opening this question, but the point is that if I write

NSolve[delta[k] + phi[k] == 0 && k >= 0, k, Reals]

the process seems to be too long. I waited 5 or 10 minutes whereupon I aborted the process. Looking at the plot of the two functions, I see that all the solutions are in the range k<=120. If, instead I write

 NSolve[delta[k] + phi[k] == 0 && 0 <= k < 120, k, Reals]

similar to how you suggested, it finds in 1 minute all these solutions:

{{k -> 0.}, {k -> 60.2778}, {k -> 60.7}, {k -> 61.2595}, {k -> 
   61.6217}, {k -> 61.9575}, {k -> 62.3592}, {k -> 62.681}, {k -> 
   62.8904}, {k -> 63.3799}, {k -> 63.7477}, {k -> 64.1951}, {k -> 
   64.5501}, {k -> 64.8526}, {k -> 65.1896}, {k -> 65.4012}, {k -> 
   65.7166}, {k -> 66.4009}, {k -> 66.93}, {k -> 67.2742}, {k -> 
   67.5128}, {k -> 67.8014}, {k -> 68.1651}, {k -> 68.6514}, {k -> 
   68.9567}, {k -> 69.2618}, {k -> 69.6373}, {k -> 69.9291}, {k -> 
   70.4181}, {k -> 70.9846}, {k -> 71.1348}, {k -> 71.5875}, {k -> 
   71.7683}, {k -> 72.0938}, {k -> 72.349}, {k -> 72.5804}, {k -> 
   73.1449}, {k -> 73.418}, {k -> 73.7302}, {k -> 73.9907}, {k -> 
   74.2564}, {k -> 74.5236}, {k -> 74.7941}, {k -> 75.2711}, {k -> 
   75.4883}, {k -> 75.7975}, {k -> 76.1651}, {k -> 76.4205}, {k -> 
   76.8743}, {k -> 77.3373}, {k -> 77.5906}, {k -> 77.9162}, {k -> 
   78.1684}, {k -> 78.8594}, {k -> 79.4396}, {k -> 79.9125}, {k -> 
   80.1565}, {k -> 80.3452}, {k -> 80.5801}, {k -> 80.8826}, {k -> 
   81.3629}, {k -> 81.8248}, {k -> 82.1206}, {k -> 82.3513}, {k -> 
   82.5568}, {k -> 82.8177}, {k -> 83.5153}, {k -> 83.7607}, {k -> 
   83.9889}, {k -> 84.1856}, {k -> 84.4299}, {k -> 84.7101}, {k -> 
   85.2929}, {k -> 85.8794}, {k -> 86.2753}, {k -> 86.5458}, {k -> 
   87.0464}, {k -> 87.4104}, {k -> 88.0525}, {k -> 88.3412}, {k -> 
   88.7546}, {k -> 89.0501}, {k -> 89.6893}, {k -> 90.0861}, {k -> 
   90.7352}, {k -> 91.2202}, {k -> 91.6101}, {k -> 92.2823}, {k -> 
   92.4872}, {k -> 94.1707}, {k -> 94.4087}, {k -> 94.8186}, {k -> 
   95.1993}, {k -> 96.4693}, {k -> 97.9171}, {k -> 98.2287}, {k -> 
   99.0925}, {k -> 101.624}, {k -> 102.032}, {k -> 105.249}, {k -> 
   108.799}}

but it missed all the solutions for k<60 !!!

I don't know why every time I insert an interval for k in which there are very many solutions, Mathematica misses all the ones in the first part of the interval. I suppose this problem can be overcome with special instruction. Maybe I have to request a finer calculation? I would insert such a command, but at the same time, I don't need all the solutions. So, to save time, I would ask Mathematica to find only the first, say, ten solutions. In this way, I don't even need to introduce an "ad hoc" upper limit for k.

How can I do this? Thank you.

POSTED BY: Davide Alfano

I don't have a solution for you; maybe you can find more help in here, here.

POSTED BY: Mariusz Iwaniuk
Posted 12 days ago

Thanks anyway. I'll take a look.

POSTED BY: Davide Alfano

You can try DiscretePlot with a short step:

DiscretePlot[phi[k] + delta[k], {k, 1, 10, 1/20}]
POSTED BY: Gianluca Gorni

The utililty RootSearch by Ted Ersek (see here) finds some more:

RootSearch[delta[k] + phi[k] == 0, {k, 1, 120}]

Unfortunately your function phi is discontinuous:

Plot[phi[k], {k, 0, 10}]

If you can foresee the singularities, then you may be able to better isolate the solutions.

POSTED BY: Gianluca Gorni
Posted 12 days ago

Thank you. I'm trying but i don't even know if I correctly loaded the package. This is the first time for me. I inserted the file RootSearch.m in the Packages folder and i wrote:

Needs["RootSearch`"]
RootSearch[delta[k] + phi[k] == 0, {k, 0, 120}]

After 7 error messages he found some solutions:

{{k -> 0.}, {k -> 9.0469}, {k -> 10.9726}, {k -> 12.6198}, {k -> 
   14.239}, {k -> 15.6005}, {k -> 17.8431}, {k -> 19.0224}, {k -> 
   20.0168}, {k -> 20.991}, {k -> 21.8113}, {k -> 22.7657}, {k -> 
   23.786}, {k -> 25.1573}, {k -> 26.0883}, {k -> 27.5029}, {k -> 
   28.1989}, {k -> 29.0041}, {k -> 30.8606}, {k -> 31.5041}, {k -> 
   32.2777}, {k -> 34.0496}, {k -> 34.6718}, {k -> 36.2231}, {k -> 
   37.3393}, {k -> 38.9894}, {k -> 39.7927}, {k -> 41.7476}, {k -> 
   44.0833}, {k -> 45.0265}, {k -> 47.1507}, {k -> 48.4782}, {k -> 
   49.2027}, {k -> 49.8007}, {k -> 50.8022}, {k -> 52.8209}, {k -> 
   53.3922}, {k -> 54.6042}, {k -> 57.4013}, {k -> 58.5957}, {k -> 
   60.2778}, {k -> 61.2595}, {k -> 62.681}, {k -> 65.1896}, {k -> 
   66.4009}, {k -> 67.5128}, {k -> 68.6514}, {k -> 69.9291}, {k -> 
   70.9846}, {k -> 71.5875}, {k -> 72.349}, {k -> 73.1449}, {k -> 
   74.2564}, {k -> 75.4883}, {k -> 76.4205}, {k -> 77.3373}, {k -> 
   77.9162}, {k -> 78.8594}, {k -> 79.9125}, {k -> 81.3629}, {k -> 
   82.1206}, {k -> 83.5153}, {k -> 85.2929}, {k -> 87.4104}, {k -> 
   88.0525}, {k -> 88.7546}, {k -> 89.6893}, {k -> 90.7352}, {k -> 
   91.6101}, {k -> 94.1707}, {k -> 95.1993}, {k -> 96.4693}, {k -> 
   97.9171}, {k -> 99.0925}, {k -> 101.624}, {k -> 105.249}, {k -> 
   108.799}}

It seems there are some solutions missing.

POSTED BY: Davide Alfano

With finer subdivisions it finds 177 solutions:

Union[Flatten@
  Table[RootSearch[delta[k] + phi[k] == 0, {k, a, a + 10}], {a, 0.01, 
    120, 10}]]
% // Length

You may try to transform your equation symbolically using some trigonometric identities, such as

ArcTan[a] + ArcTan[c] == ArcTan[( a + c)/(1 - a c)]

(in its domain of validity).

POSTED BY: Gianluca Gorni
Posted 9 days ago

Thank you. For now, it seems the most reliable method. Nevertheless, I'm not completely satisfied because this requires previously studying the equation to know how many solutions there are in the interval and then deciding how finer the subdivision must be. I was looking for an almost automatic procedure to implement. You have been helpful anyway. Thank you.

POSTED BY: Davide Alfano
Reply to this discussion
Community posts can be styled and formatted using the Markdown syntax.
Reply Preview
Attachments
Remove
or Discard

Group Abstract Group Abstract