Take the 2-minute tour ×
Mathematics Stack Exchange is a question and answer site for people studying math at any level and professionals in related fields. It's 100% free, no registration required.

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.

share|improve this question
    
I think your equation to solve is $$f(x) = 1-x \cot{x} - Bi = 0$$ –  Ron Gordon May 29 '13 at 21:23

2 Answers 2

up vote 0 down vote accepted

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.

share|improve this answer

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.

share|improve this answer
    
Thanks @RonGordon. I was wondering whether Bisection method will be computationally efficient compared to other numerical methods in this specific scenario and also whether other methods like Newton-Raphson will give me a sure convergence. –  L'Unità May 30 '13 at 5:24
    
In this case, bisection is better because within a subinterval, the function is monotonic and the root simple. You really should be able to get your root in about $\log_2{1/\epsilon}$, where $\epsilon$ is the specified tolerance. –  Ron Gordon May 30 '13 at 12:52
    
I tried out secant method as recommended by @ross in the other answer. It worked out well as compared to the crude method I was following to find the root. I think I should do a comparative study for this particular case. –  L'Unità May 30 '13 at 15:42

Your Answer

 
discard

By posting your answer, you agree to the privacy policy and terms of service.

Not the answer you're looking for? Browse other questions tagged or ask your own question.