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:


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:

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?

  • $\begingroup$ I guess the question is equivalent to "can I create a function which efficiently returns a list of the small positive roots of $\epsilon\lambda=\cot\left(\lambda\right)$?". If you're fine with generating a 3D cube of data (might need a lot of RAM?), then you can add together the tensor products of the sinusoids for each of the root frequencies, and then save the array to disk. By doing this, you generate slices of the 3D cube for each value of $\epsilon$. Finally, stitch them together and call Interpolation on it, giving an InterpolatingFunction $u(x,t,\epsilon)$ as you desired. $\endgroup$ Commented Mar 2, 2015 at 0:00
  • $\begingroup$ At the moment I don't have time to test this (probably won't be free until tomorrow evening), but this might be one possible avenue of approach. $\endgroup$ Commented Mar 2, 2015 at 0:00
  • $\begingroup$ I guess the crux of what I'm trying to say is that because each term in the summand is separable, it is not necessary to compute the value of the term over an entire $(x,t)$ grid, and you can use TensorProduct instead, which ought to reduce computation time by an order of magnitude or so. $\endgroup$ Commented Mar 2, 2015 at 0:13

1 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]}]

first three roots of ελ = cotλ

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

plot of partial sum

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

  • $\begingroup$ wow +1 learning a lot from your code :) $\endgroup$
    – ubpdqn
    Commented Oct 24, 2015 at 4:27

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.