Finding and plotting roots with parameter

Question

I have the following function:

$$u(x,t)=\sum_{k=1}^\infty\frac{v_0 L}{\lambda_k c \sin\left(\lambda_k\right)}\sin\left(\frac{\lambda_kx}{L}\right)\sin\left(\frac{\lambda_kct}{L}\right)$$

and this equation whose solutions are the $\lambda_k$ needed for the function above:

$$\varepsilon\lambda=\cot\left(\lambda\right)$$

How can I feed the solutions of this (not analytically solvable) equation into my function? If possible, I'd like to create a function $u(x,t,\varepsilon)$ to create 3D plots. What I have found so far:

Clear[findRoots]
Options[findRoots] = Options[Reduce];
findRoots[gl_Equal, {x_, von_, bis_}, 
   prec : (_Integer?Positive | MachinePrecision | Infinity) : 
    MachinePrecision, wrap_: Identity, opts : OptionsPattern[]] := 
  Module[{work, glp, vonp, 
    bisp}, {glp, vonp, bisp} = {gl, von, bis} /. 
     r_Real :> SetPrecision[r, prec];
   work = wrap@Reduce[{glp, vonp <= x <= bisp}, opts];
   work = {ToRules[work]};
   If[prec === Infinity, work, N[work, prec]]];

uk         = (v0 L)/(lambdak c Sin[lambdak]) Sin[lambdak x/L] Sin[lambdak c t/L];

This function finds the roots in an interval manually defined since FindRoots or NSolve only finds one root. Removed the $\varepsilon$ for a test:

lambdalist = x /. findroots[x == Cot[x], {x, 1, 100}];

This gives me a list like this to use:

{0.860334, 3.42562...}

Which works with my function

u[x_,t_,L_,c_,v0_] = Sum[uk, {lambdak, lambdalist}];

My question is:

Is it possible to add $\varepsilon$ as a parameter or do I have to define it beforehand to find the values?

Answer 1

Can I create a function which efficiently returns a list of the small positive roots of $\varepsilon \lambda=\cot\lambda$?

Here is a routine I wrote quite a while back, when I was doing some research related to the square well potential. The procedure does the computation in two stages: an initial approximation computed through the Delves-Lyness contour integral method, and then a subsequent polishing with Halley's method.

SetAttributes[cotSol, {Listable, NHoldFirst, NumericFunction}];
cotSol[ε_] := cotSol[1, ε];
cotSol[k_Integer?Positive, ε_?InexactNumberQ] :=
Module[{prec = Precision[ε], r = π/4, h, w},
       h = π (k - UnitStep[ε]/2 - 1/4);
       w = Re[NIntegrate[# ((ε + 1) Sin[#] + ε # Cos[#])/
                         (ε # Sin[#] - Cos[#]) &[h + r Exp[I t]] Exp[I t],
                         {t, 0, 2 π}, AccuracyGoal -> 1,
                         Method -> {"Trapezoidal", "SymbolicProcessing" -> 0},
                         PrecisionGoal -> 2, WorkingPrecision -> prec]/8];
       w = FixedPoint[Block[{cw = Cos[#], sw = Sin[#], f, fp},
                            f = ε # sw - cw; fp = (ε + 1) sw + ε cw #;
                            # - 2 f fp/(2 fp^2 + f (ε # sw - (2 ε + 1) cw))] &,
                      SetPrecision[w, If[prec === MachinePrecision,
                                         prec, prec + 5]]];
       N[w, prec]]

Here's a plot of the first three solutions for different values of $\varepsilon$:

Plot[cotSol[Range[3], ε], {ε, -5, 5}, Evaluated -> True, 
     PlotStyle -> {RGBColor[7/19, 37/73, 22/31], RGBColor[59/67, 11/18, 1/7],
                   RGBColor[14/25, 9/13, 8/41]}]

---

As an example, let us plot a partial sum of the series given in the OP:

With[{ε = 1/10, c = 3, n = 18}, 
     Plot3D[Sum[With[{λ = N[cotSol[k, ε]]}, 
                     Sin[λ x] Sin[λ c t]/(λ c Sin[λ])], {k, n}] // Evaluate,
            {x, 0, 2 π}, {t, 0, 2 π}, Mesh -> False, PlotPoints -> 55]]

I'll leave the determination of how to sum the infinite series to full accuracy to somebody else.