# Solving and Plotting Infinite sums over solutions to transcendental equation involving Bessel functions

I am trying to use Mathematica to plot the following equation $E(q,\Delta)$ vs $qa$, for set values of $D \Delta/a^2$ and $\bar\rho a/D$. I have looked at other stackexchange posts which deal with transcendental equations and Bessel functions, but none have explained exactly what I am trying to do.

Where $\beta _{nk}$ is defined by the transcendental equation:

$\Large \beta_{nk}\frac{J'_{n}(\beta_{nk})}{J_{n}(\beta_{nk})}=-\frac{\bar \rho a}{D}$

Here is the code for the first equation:

E1[q_, \[CapitalDelta]_, rho_, k_, D_, a_, n_] := (4*Exp[(-((B0k[k])^2)*D*\[CapitalDelta])/a^2]*(B0k[k])^2/(((rho*a)/D)^2 + (B0k[k])^2)*((q*a) (-BesselJ[0, q*a]) + rho*BesselJ[0, q*a])^2/(q*a^2 - B0k^2)^2) + (8*Exp[(-((B[n, k])^2)*D*\[CapitalDelta])/a^2]*(B[n, k])^2/(((rho*a)/D)^2 + (B[n, k])^2 - n^2)*((q*a)*BesselJ[n, q*a] + (rho*a)/D*BesselJ[n, q*a])^2/((q*a)^2 - B[n, k]^2)^2)


I can plot the transcendental equation, but can't figure out how to find the solutions in such a way that I could incorporate them into the other equation.

Plot[B*(-BesselJ[1, B])/BesselJ[0, B], {B, -10, 10}]


• Please post copyable code (not images of equations). – David G. Stork May 26 '16 at 22:00
• At the point where you want this computation, do you have values for $\overline{\rho}$, $a$, and $D$, or are you attempting to do the (likely) impossible and evaluate this symbolically? – Eric Towers May 27 '16 at 5:43
• Hi Thank you for your help. I am not trying to evaluate this symbolically. Eventually I would like to fit data to this equation is possible but first would just like to plot the function for a known range of $\frac{rho a}{D}$. I have posted the requested code but it is obviously not complete because I am not sure how to code. – vwith May 28 '16 at 0:28
• Could you indicate (either in the picture of your equation, or otherwise) what are the lower bounds of your sums? – Eric Towers Jul 27 '16 at 5:36

This is two sums, which are similar, but the first is good practice for the second (and the $\beta_{0k}$s are a little different from the ones with $n>0$).

Throughout, we make the following changes in notation, consistent with your intention to treat $qa$ as the ordinate for plotting and the other two combinations as parameters:

• $qa \rightarrow x \text{ and }$ x,
• $D\Delta / a^2 \rightarrow \delta \text{ and }$ DDeltaOverASquared, and
• $\overline{\rho} a / D \rightarrow r \text{ and }$ RhoAOverD.

In this notation, the original sum is \begin{align} \hat{E}(x; \delta, r) = &\sum_{k=1}^{\infty} 4 \mathrm{e}^{-\beta_{0k}\delta}\frac{\beta_{0k}^2}{r^2 + \beta_{0k}^2 \phantom{{}-n^2}} \frac{(x J_0'(x) + r J_0(x))^2}{(x^2 - \beta_{0k}^2)^2} \\ + \sum_{n=1}^\infty & \sum_{k=1}^\infty 8 \mathrm{e}^{-\beta_{nk}\delta}\frac{\beta_{nk}^2}{r^2 + \beta_{nk}^2 - n^2} \frac{(x J_n'(x) + r J_n(x))^2}{(x^2 - \beta_{nk}^2)^2} \end{align}

# Facts about $\beta_{nk}$

For each $n$, $\beta_{nk}$ is a strictly monotonically increasing function of $k$. For very large $x$, $J_n(x) \approx \sqrt{\frac{2}{\pi x}} \cos \left( x - \frac{n \pi}{2} - \frac{\pi}{4} \right)$.

For $n=0$ and large $\beta$, we consider $f(\beta) = \frac{J_0'(\beta)}{J_0(\beta)} = -\beta \frac{J_1(\beta)}{J_0(\beta)}$ using the cosine formula above, and see that $f(\beta)$ has a zero at local minima and maxima of $J_0(\beta)$ and has a vertical asymptote at zeroes of $J_0(\beta)$, and the zeroes and asymptotes are interleaved. With some effort, we can show that this function is strictly monotonically decreasing between the vertical asymptotes on $(0,\infty)$. Since $0 \cdot f(0) = 0$, the zeroes of $f$ (equivalently, the zeroes of $J_1$) partition $[0,\infty)$ into intervals each of which is injectively mapped to $(\Bbb{R} \setminus \{0\})$ (ignoring the asymptote, since we assume $r$ is finite). The first few zeroes of $\beta f(\beta)$ are $\beta = 0$ and

NSolve[{-β (BesselJ[1, β]/BesselJ[0, β]) == 0, 1 < β < 14}, β]
(*  {{β -> 3.83171}, {β -> 7.01559}, {β -> 10.1735}, {β -> 13.3237}}  *)


the last two of which already differ by only $3.15...$ and the difference between subsequent zeroes gets closer and closer to $\pi$, the spacing between the zeroes of sine (the derivative of the cosine in the asymptotic formula). Our takeaway from this is $$(k-1)\pi \leq \beta_{0k}$$ and we will use this estimate to bound the error in truncating the first sum. We package the above into a machine which can find these $\beta$s:

Clear[B];
B[0, 1, 0, opts___] = 0;
B[0, k_ /; k > 1, 0, opts___] := BesselJZero[1, k - 1];
B[0, k_ /; k >= 1, RhoAOverD_ /; RhoAOverD > 0, opts___] :=
z /. FindRoot[
-z^2 BesselJ[1, z]/BesselJ[0, z] == -RhoAOverD,
{z,
If[k > 1, BesselJZero[1, k - 1], 0] + 1,
If[k > 1, BesselJZero[1, k - 1], 0],
BesselJZero[1, k]
}, opts];
B[0, k_ /; k >= 1, RhoAOverD_ /; RhoAOverD < 0, opts___] :=
z /. FindRoot[
-z^2 BesselJ[1, z]/BesselJ[0, z] == -RhoAOverD,
{z,
BesselJZero[1, k] - 1,
If[k > 1, BesselJZero[1, k - 1], 0],
BesselJZero[1, k]
}, opts];


We handle the first zero specially because FindRoot[] seems to have some difficulty locating the double root at zero. Also, since every one of our partitions of $(0,\infty)$ has an asymptote in it, we search for positive values of $r$ on one side of the asymptote and for negative values on the other side. The partitions are at least $\pi$ units wide, so backing into the partition by $1$ unit from the end (on the side of the asymptote where we expect to find the solution) is safe. Since the zeroes of the function are on the borders of two partitions, they're called out specially too (also since the extra zero at zero causes the zeroes of this function to be offset by one from the zeroes of $J_1$). Finally, at some point, we may wish to improve the accuracy/precision of a $\beta$, so we can pass relevant options in to FindRoot[].

B[0, 1, 0]
B[0, 2, 0]
B[0, 1, 1]
B[0, 1, -1]
B[0, 1, 1, WorkingPrecision -> 30]
B[0, 1, -1, WorkingPrecision -> 35]

(*  0  *)
(*  BesselJZero[1,1]  *)
(*  1.17726  *)
(*  3.76181  *)
(*  1.17725649964807633454475810267  *)
(*  3.7618061255679366126466190513942089  *)


For $n>0$ and large $\beta$, we consider $\beta f_n(\beta) = \beta \frac{J_n'(\beta)}{J_n(\beta)} = -n + \beta\frac{J_{n-1}(\beta)}{J_n(\beta)}$ and repeat our reasoning using the cosine formula above. This analysis is largely unchanged, except everything is shifted by $n$ units and there is a zero-free gap in $J_{n-1}(\beta)$, so the first partition is wider. Our estimates are \begin{align} 0 &\leq \beta_{n1} & n &\geq 1 \\ \frac{(n - 1) \pi}{2} + (k - 1)\pi = (2k + n -3)\frac{\pi}{2} &< \beta_{nk} & n &\geq 1, k > 1 \text{.} \end{align}

And we use these observations to put together the rest of the gadgets we need to find the $\beta$s. The only significant change is the shift by $n$ and having to move the condition between RhoAOverD and n to the end of the definition (because n cannot be used in a Condition at the point that RhoAOverD is declared).

B[n_ /; n > 0, 1, RhoAOverD_, opts___] := 0  /; RhoAOverD == -n ;
B[n_ /; n > 0, k_ /; k > 1, RhoAOverD_, opts___] :=
BesselJZero[1, k - 1]  /; RhoAOverD == -n ;
B[n_ /; n > 0, k_ /; k >= 1, RhoAOverD_, opts___] := (
z /.
FindRoot[-n + (z BesselJ[-1 + n, z])/BesselJ[n, z] == -RhoAOverD,
{z,
If[k > 1, BesselJZero[n - 1, k - 1], 0] + 1,
If[k > 1, BesselJZero[n - 1, k - 1], 0],
BesselJZero[n - 1, k]},
opts])  /; RhoAOverD < -n ;
B[n_ /; n > 0, k_ /; k >= 1, RhoAOverD_, opts___] := (
z /.
FindRoot[-n + (z BesselJ[-1 + n, z])/BesselJ[n, z] == -RhoAOverD,
{z,
BesselJZero[n - 1, k] - 1,
If[k > 1, BesselJZero[n - 1, k - 1], 0],
BesselJZero[n - 1, k]},
opts])  /; RhoAOverD > -n ;

B[1, 1, 0]
B[1, 2, 0]
B[1, 1, 1]
B[1, 1, -1]
B[1, 1, 1, WorkingPrecision -> 30]
B[1, 1, -1, WorkingPrecision -> 35]

(*  1.84118  *)
(*  5.33144  *)
(*  2.40483  *)
(*  0  *)
(*  2.40482555769577276862163187933  *)
(*  0  *)


(Note that 0 is already an infinite precision zero.)

As a future direction for improvement, starting the FindRoot searches one unit from the end of the interval containing the $\beta$s is not the best way to go. The asymptotic estimate for $J_n$ suggests we should use something vaguely like $\beta \tan \beta$, tweaked so its zero and asymptote match in the interval we are searching to get a starting estimate. This will likely make the FindRoot much quicker, especially for large $k$. However, solving this tangent equation will cost significantly more time, so instrumentation and testing are in order.

# Aside

We actually have enough to plot something now.

summand0[x_, DDeltaOverASquared_, RhoAOverD_] =
4 E^(-B[0, k, RhoAOverD] DDeltaOverASquared)
B[0, k, RhoAOverD]^2/(RhoAOverD^2 + B[0, k, RhoAOverD]^2)
(x D[BesselJ[0, x], x] + RhoAOverD BesselJ[0, x])^2/(x^2 - B[0, k, RhoAOverD]^2)^2
summandN[x_, DDeltaOverASquared_, RhoAOverD_] =
8 E^(-B[n, k, RhoAOverD] DDeltaOverASquared)
B[n, k, RhoAOverD]^2/(RhoAOverD^2 + B[n, k, RhoAOverD]^2 - n^2)
(x D[BesselJ[n, x], x] + RhoAOverD BesselJ[n, x])^2/(x^2 - B[n, k, RhoAOverD]^2)^2

eHat10[x_, DDeltaOverASquared_, RhoAOverD_] =
Sum[summand0[x, DDeltaOverASquared, RhoAOverD], {k, 1, 10}] +
Sum[summandN[x, DDeltaOverASquared, RhoAOverD], {k, 1, 10}, {n, 1, 10}];
eHat11[x_, DDeltaOverASquared_, RhoAOverD_] =
Sum[summand0[x, DDeltaOverASquared, RhoAOverD], {k, 1, 11}] +
Sum[summandN[x, DDeltaOverASquared, RhoAOverD], {k, 1, 11}, {n, 1, 11}];

Plot[Evaluate[eHat10[x, 1, 1]], {x, 0, 10}]
Plot[Evaluate[eHat11[x, 1, 1]], {x, 0, 10}]
Plot[Evaluate[eHat11[x, 1, 1] - eHat10[x, 1, 1]], {x, 0, 10}]


These are just truncations to $110$ terms (eHat10) and $132$ terms (eHat11). They plot the same and differ by less than $\sim 10^{-9}$ over the plotted interval. So for parameters near $1$, the convergence is rather quick. We don't yet know that the rest of the terms don't eventually accumulate significant error, but we will find out that they do not for parameters near $1$.

Given the large range of values appearing in the plots, a LogPlot might be more informative.

LogPlot[Evaluate[eHat11[x, 1, 1]], {x, 0, 10}, PlotRange -> All]


We use Evaluate[] in the plots so that the finite sums will be expanded and the $\beta$s evaluated once and the resulting rational expressions used repeatedly. This is much faster than repeating the calculation of the $\beta$s every time Plot evaluates its first argument. (Don't try to FullSimplify in this Evaluate call. The kernel will spend more time trying to simplify, without much effect, than it does plotting without simplifying.)

This might be enough to tell you everything you want to know. But generally, we need confidence that the truncated sums are actually close to the infinite sums and have some bound on the discrepancy. Let's get those error estimates and make some functions that give us values with bounds.

# The first sum

Now we analyze $$S_0(K) = \sum_{k=K+1}^{\infty} 4 \mathrm{e}^{-\beta_{0k}\delta}\frac{\beta_{0k}^2}{r^2 + \beta_{0k}^2 \phantom{{}-n^2}} \frac{(x J_0'(x) + r J_0(x))^2}{(x^2 - \beta_{0k}^2)^2}$$ to estimate the error in truncating the first sum after $K$ terms. From the above, $$\mathrm{e}^{- \beta_{0k} \delta} \leq \mathrm{e}^{-(k-1)\pi \delta} \text{.}$$ We could fuss with the fraction in $\beta$, but it is sufficient to estimate it with $$\frac{\beta_{0k}^2}{r^2 + \beta_{0k}^2} \leq 1 \text{.}$$ Further, both $J_0$ and $J_0'$ are bounded by $1$, so for $\beta_{0k} > (k-1)\pi > \sqrt{2}x$, $$\frac{(x J_0'(x) + r J_0(x))^2}{(x^2 - \beta_{0k}^2)^2} = \frac{(x J_0'(x) + r J_0(x))^2}{(\beta_{0k}^2 - x^2)^2} \leq \frac{(x J_0'(x) + r J_0(x))^2}{x^4} \leq \frac{(|x| + |r|)^2}{x^4} \text{.}$$ This gives us our first indication of a lower bound on $K$ : $$K > \frac{\sqrt{2}}{\pi} x \text{.}$$ Putting these together, we have $$S_0(K) \leq \sum_{k=K+1}^{\infty} 4 \mathrm{e}^{-(k-1)\pi \delta}\cdot 1 \cdot \frac{(|x| + |r|)^2}{x^4} \qquad\qquad [ K > \frac{\sqrt{2}}{\pi} x ] \text{.}$$ Pulling out everything not depending on $k$, we are left with a simple geometric series: $$S_0(K) \leq 4 \frac{(|x| + |r|)^2}{x^4} \sum_{k=K+1}^{\infty} \mathrm{e}^{-(k-1)\pi \delta} = 4 \frac{(|x| + |r|)^2}{x^4} \frac{\mathrm{e}^{(1-K)\pi \delta}}{\mathrm{e}^{\pi \delta} - 1} \qquad\qquad [ K > \frac{\sqrt{2}}{\pi} x ] \text{.}$$ Of course, we can make this sum as small as we like by increasing $K$. (Happily, this error estimate gets smaller as $x$ increases, starting once $x \approx \sqrt{r}$, so we do not have to take longer and longer sums as $x$ increases.)

(We could do better. The same sorts of estimates for the $\beta$s could be turned into upper bounds, which would give us a lower bound for the sum $S_0(K)$. This lower bound would be of about the same shape as the estimate we have, it would just skip a few terms. Once the lower bound and upper bound agreed sufficiently, we could stop increasing $K$. This would happen earlier than the method we are currently using, taking $0$ as the lower bound.)

Let's capture this in code

summand0[x_, DDeltaOverASquared_, RhoAOverD_] =
4 E^(-B[0, k, RhoAOverD] DDeltaOverASquared)
B[0, k, RhoAOverD]^2/(RhoAOverD^2 + B[0, k, RhoAOverD]^2)
(x D[BesselJ[0, x], x] + RhoAOverD BesselJ[0, x])^2/(x^2 - B[0, k, RhoAOverD]^2)^2
sum0K[K_, x_, DDeltaOverASquared_, RhoAOverD_] :=
Sum[summand0[x, \[Delta], r], {k, 1, K}]
sum0Error[K_, x_, DDeltaOverASquared_, RhoAOverD_] :=
4 (Abs[x] + Abs[RhoAOverD])^2/x^4
E^((1 - K) Pi DDeltaOverASquared)/(E^(Pi DDeltaOverASquared) - 1)


And let's see what the error does versus $x$.

LogPlot[sum0Error[3, x, 1, 1], {x, 0, 10}]


... and we find the error is nicely shrinking as $x$ increases.

# The second sum

There's an easy way to find cutoffs for double sums and a tricky way that generally uses half as many terms. The easy way is to find upper bounds $N$ and $K$ and sum on the rectangle $(n,k) \in [1,N]\times[1,K]$. A slightly smarter way is to find a line in the $N$-$K$ plane and only sum the terms bounded by the axes and this line. As a gross generalization, this line typically runs along the diagonal of the rectangle. I'm a little lazy (although maybe that's not obvious given the length of this answer and the endless stream of parenthetical comments), so I'll go with the two cutoffs.

Sadly, it's not quite as simple as that. We don't have shrinking bounds on the $\beta_{n1}$s because they could be in the interval prior to the first asymptote, which means we only know they're greater than or equal to zero. We'll have to handle $k=1$ separately, but the above comments about cutoffs $N$ and $K$ hold for the $\sum_{k=2}^K \dots$ sums.

We wish to study these sums to estimate the truncation error. For $K >1$, \begin{align} S_K(N,K) &= \sum_{n=1}^{N} \sum_{k=K+1}^\infty 8 \mathrm{e}^{-\beta_{nk}\delta}\frac{\beta_{nk}^2}{r^2 + \beta_{nk}^2 - n^2} \frac{(x J_n'(x) + r J_n(x))^2}{(x^2 - \beta_{nk}^2)^2} \text{,} \\ S_N(N) &= \sum_{n=N+1}^\infty \sum_{k=2}^\infty 8 \mathrm{e}^{-\beta_{nk}\delta}\frac{\beta_{nk}^2}{r^2 + \beta_{nk}^2 - n^2} \frac{(x J_n'(x) + r J_n(x))^2}{(x^2 - \beta_{nk}^2)^2} \text{, and } \\ S_1(N_1) &= \sum_{n=N_1+1}^\infty 8 \mathrm{e}^{-\beta_{n1}\delta}\frac{\beta_{n1}^2}{r^2 + \beta_{n1}^2 - n^2} \frac{(x J_n'(x) + r J_n(x))^2}{(x^2 - \beta_{n1}^2)^2} \text{.} \end{align}

The first two we analyze as for the first sum: \begin{align} \mathrm{e}^{-\beta_{nk} \delta} &\leq \mathrm{e}^{(2k+n-3)\pi \delta/2} \text{,} \\ \frac{\beta_{nk}^2}{r^2 + \beta_{nk}^2 - n^2} &\leq \frac{\pi^2}{\pi^2 - 4} \text{, and} \\ \frac{(x J_n'(x) + r J_n(x))^2}{(x^2 - \beta_{nk}^2)^2} &\leq \frac{|x| + |r|}{x^4} & &[K \geq \frac{3-n}{2} + \frac{\sqrt{2}}{\pi} x \geq 1 + \frac{\sqrt{2}}{\pi} x ] \text{.} \end{align} As before, we used the lower bound on the $\beta$s for the exponential and the exact same estimating for the Bessel functions fraction. The fraction in $\beta$s is a little trickier. Note that as $n$ increases to $r$, the ratio increases to approximately $1$. Subsequently, the $\beta$s increase by a bit more than $\pi/2$ (asymptotically decreasing to $\pi/2$) every time $n$ increases by $1$ or by a bit more than $\pi$ (asymptotically decreasing ot $\pi$) every time $k$ increases by $1$. Looking at $\lim_{k \rightarrow \infty} \frac{(\beta + k \alpha \pi)^2}{r^2 + (\beta + k \alpha \pi)^2 - (r+k)^2} = \frac{\pi^2}{\pi^2 - (1/\alpha)^2}$, which is maximized for $\alpha \in [1/2,\infty)$ by $\alpha = 1/2$, we find that $\frac{\pi^2}{\pi^2 - 4}$ is a bound we can use. Then summing the geometric series, we have \begin{align} S_K(N,K) &\leq \frac{8 \pi^2}{\pi^2 - 4} \frac{|x| + |r|}{x^4} \frac{\mathrm{e}^{-(N+2K-3)\delta \pi/2} \left( \mathrm{e}^{N \delta \pi / 2} - 1 \right)}{ \left( \mathrm{e}^{\delta \pi / 2} - 1 \right)^2 \left( \mathrm{e}^{\delta \pi/2} + 1 \right)} & &[K \geq 1 + \frac{\sqrt{2}}{\pi} x] \\ S_N(N) &\leq \frac{8 \pi^2}{\pi^2 - 4} \frac{|x| + |r|}{x^4} \frac{\mathrm{e}^{-(N-1) \delta \pi / 2}}{\left( \mathrm{e}^{\delta \pi / 2} - 1 \right)^2 \left( \mathrm{e}^{\delta \pi/2} + 1 \right)} & &[K \geq 1 + \frac{\sqrt{2}}{\pi} x] \text{.} \end{align}

This is not straightforward. You have some hope of stumbling into approximately correct results because of the aggressively decreasing exponentials in your sum, but I make no promises that you can get any particular degree of precision.

You should probably construct little machines for each $n$ that sequentially produce values of $\beta_{n,k}$ for $k=1,2,3, \dots$. As an example of how challenging you may find this, for $n \in [1,5]$, and $\dfrac{\overline{\rho} a}{D} = 1$, the $\beta_{n,k}$ in the interval $(0,5\pi)$ are

Table[
With[{lhs = B D[BesselJ[n, B], B]/BesselJ[n, B]},
Solve[{
lhs == -1,
0 < B < 5 Pi, Element[B, Reals]},
B]],
{n, 1, 5}]

(* Much mess and warnings *)

N[%]

(*  {{{B -> 2.40483}, {B -> 5.52008}, {B -> 8.65373}, {B -> 11.7915}, {B -> 14.9309}},
{{B -> 3.51832}, {B -> 6.86626}, {B -> 10.073}, {B -> 13.2477}},
{{B -> 4.61256}, {B -> 8.15781}, {B -> 11.44}, {B -> 14.6571}},
{{B -> 5.69431}, {B -> 9.41276}, {B -> 12.7689}},
{{B -> 6.76717}, {B -> 10.6411}, {B -> 14.0687}}}  *)


So starting with a family of gadgets, one for each $n$, each of which spits out $\beta_{k,1}, \beta_{k,2}, \beta_{k,3}, \dots$, you will want to construct your sums. You'll want some way to bound the portion of each sum you're tracking that you have not summed and also bound all the sums you haven't started yet. (I have not worked through the details -- replacing $\beta_{n,k}$ with $k$ seems like a feasible strategy for getting awfully loose bounds. Alternatively, you could just make a partial sum function where you provide upper bounds on $n$ and $k$ and then call that from another function that decides whether increasing either bound seems relevant given your choices of AccuracyGoal and PrecisionGoal. It'd be great if there were a built-in function that magically could do this for black-box functions like the partial sum function I'm describing above, but we don't have magic.

• i am still working on this but was not successful in the execution of this method. – vwith Jul 15 '16 at 0:35