Numerical method to solve a trigonometric (cotangent) function - transient heat transfer problem

Question

I was trying to develop a mathematical model for transient one-dimensional heat conduction of spheres using approximate analytical solution as mentioned in Cengel{refer page number 229 in that document}. To obtain the temperature gradient I will have to find the first few roots of the Eigen function (refer Table 4.1 in the above document). For the sake of the others who do not want to download the above PDF I am writing down the equation below.

$f(\lambda) = 1-\lambda_{n} cot \lambda_{n}-Bi=0$, where $Bi\geq 1$ (Biots number) and $\lambda>1$.

The issue with the above equation is that for various values of $Bi$ this equation will have roots and singularities over a range. A typical example for $Bi=5$ and $0\leq\lambda\leq10$ is shown in Recktenwald, 2006 {refer Figure 2 in page number 3 in this document}. Please keep in mind that $\lambda$ - the eigen value myself and Cengel have used is called $\zeta$ in the above figure taken from Recktenwald, 2006 linked above.

Can you please suggest me an efficient numerical way of finding roots which also avoids singularity. I am planning to implement this in C. As the program given in the appendix of Recktenwald, 2006 is in Matlab, I couldn't understand the underlying logic as I can't handle it well. I will be doing my code with gcc. Please help to write the algorithm.

Answer 1

If you look at the plots, you need $\cot \lambda_n$ to be pretty small because the other terms are pretty small, but $\lambda_n$ can be large. The singularities come at $\lambda_n=k\pi$, the zero crossings for $\cot$ are at $(k+\frac 12)\pi$. This suggests starting your root finding at $(k+\frac 12)\pi$. As the function is less than zero at those points, pick a value near $(k+1)\pi$ and you will have bracketed the root. Now bisection, the secant method, or Brent's method should have no trouble at all. Newton-Raphson may have trouble because the function is almost horizontal near $(k+\frac 12)\pi$, particularly for the lowest root, but once you get close you could switch to it if you want.

Answer 2

My advice (I am not going to write your algorithm) is to break the interval of interest up into subintervals between the poles of the function $\cot{x}$, which are at values of $x=k \pi$ for integer $k$. Within each subinterval, there is a single root of the function. To find that root in, say, the $k$th subinterval, you may simply use a bisection routine as follows.

Begin with bracket $[(k-1) \pi + \epsilon, (k-1/2) \pi, k \pi - \epsilon]$ for some $\epsilon = 0.01$ or something like that (you may want to scale it with, say, $Bi$. Is $(k-1/2) \pi$ the same sign as the left endpoint or the right? Replace that endpoint with $(k-1/2) \pi$, determine the midpoint of the new bracket, and iterate until you find the root to some desired tolerance. That gets you the root for that subinterval. Then move onto the next one, etc.