Part 1 – Elliptic Integrals

Two of the most fundamental inductance formulae are those of Lorenz and Maxwell. Both of these formulae are expressed in terms of elliptic integrals, and both can and have been used as the basis for other inductance formulae. We will start with Lorenz.

A) Lorenz, Nagaoka and Solenoidal Current Sheets

Lorenz's formula is:

where K(k) and E(k) are, respectively, the complete elliptic integrals of the first and second kind, of modulus k where:

r being the radius,

ℓ being the length, and

N being the number of turns.

This formula is typical of those from the late 19th century. The inductance unit, the Henry, was not yet in common use. In order to use the formula with modern SI units, the term outside of the square brackets should include the factor µ0/4π which is exactly equal to 10-7. Then the result will be in Henries, when the radius and length are given in meters. Converted to modern SI units formula (1) becomes:

which now works in units of meters and Henries. And where,

µ0 = 4π×10-7 (The definition of the permeability of free space.)

While we are defining the modulus k, we will also define the complementary modulus k' which regularly shows up in elliptic integral calculations:

(The notation used here treats the elliptic integrals as standard functions, in that the argument is shown in parentheses immediately following the function name. It's very common, especially in older texts, simply to write K and E, without showing the argument, often leaving the reader to puzzle over what the argument should be.)

The formula as written does not look overly complicated, except that we haven't discussed the calculation of the elliptic integrals yet, and this is where things can get a bit messy. In the past, this formula has been considered impractical to calculate due to E and K, and scientists have gone to great lengths to simplify the formula in order to eliminate their calculation. However, it will be shown that the calculation of elliptic integrals is not terribly painful, and as a result, will allow us to avoid the approximations of the simplified formulae.

In 1909, Nagaoka [3] rearranged Lorenz's formula like so:

Note that in some of the early literature, uppercase N is total turns while lowercase n is the turns per unit length. In this discussion, all of the formulae have been converted to total turns. Also note that the factor of µ0/4π = 10-7 needs to be included outside the square brackets, in order to work in meters and Henries.

Formula (2) can be written as:

where kL is now known as Nagaoka's coefficient, and as is evident from above, is given by the last half of formula (2):

If we put the above formula in terms of the ratio of diameter to length, 2r/ℓ, which we will define as u, we are now dealing in relative rather than absolute dimensions (essentially a shape factor), and the formula for kL becomes:

The result of all this is that kL can now be calculated once and for all in terms of the coil shape factor u. Nagaoka [3] did this and tabulated the results. Thus, when calculating the inductance of a coil, a person only needed to calculate the diameter to length ratio u, use it to look up the value of kL from the table, and then plug the kL value into the simple formula (3) to determine the inductance. Of course, we are trying to avoid the use of tables, so our intent is to find a quick way to calculate the elliptic integrals. Having done this, we could use any of formulae (1), (2) or (3). However formula (3) will be the choice because it is in a form which makes it convenient to compare with other inductance formulae. This will be useful later when we want to evaluate the accuracy of various empirical formulae.

My first attempt at coding a BASIC program to calculate elliptic integrals began by referring to Dwight [4] (Art. 773.1-774.2) who gives several series formulae for elliptic integrals. Using some care, and proper optimization, these can provide accurate and fast results. I've described this in some detail in [5] and won't repeat it here, as I have since found a better method which will be described below. Incidentally, Dwight's book is typical of mathematical references written before the days of computers. It provides formulae in closed analytical form or as series expansions, but not in an algorithmic form which may sometimes be more useful for computer application. Judging from many of the technical papers on inductance calculations written since computers have become ubiquitous, one can see that many have fallen into the trap of assuming there is no better method of calculation. However, in 1921, Louis V. King [6] proposed a very efficient technique using the arithmetic-geometric mean method. Grover applied his method in a 1928 paper [7]. Since then, the method appears to have been overlooked by many, but has resurfaced a few times, most recently in papers by Miller, 1987 [8] and de Queiroz, 2003 [9]. It was in the latter paper where I first encountered the method. However, neither Miller nor de Queiroz appear to have been aware of King's work. This is hardly surprising considering how long ago it was published. However there are useful details in King's paper which make for very efficient computer implementation.

To describe the arithmetic-geometric mean method, first we need to explain what an arithmetic-geometric mean is. Everyone knows that the arithmetic mean of two numbers, aka the average, is simply the numbers added together and then divided by two. The geometric mean is found by multiplying the two numbers together and then taking the square root.

If we start with two numbers a0 and b0, then we can calculate the arithmetic-geometric mean by repetitively calculating both the arithmetic mean a1, and the geometric mean b1, and repeating this operation until the resulting numbers converge to the same value. In mathematical terms:

an=(an-1+bn-1)/2

bn=√(an-1 × bn-1)

Hence, an is the arithmetic mean of previous values an-1 and bn-1, while bn is the geometric mean of previous values an-1 and  bn-1. After very few iterations, an and bn will converge to the same value, which is defined as the arithmetic-geometric mean (AGM). We will define the limiting values of an and bn as a and b.

If we start with:

a0=1

and

b0 =k'    (the complementary modulus)

then the resulting AGM value a is related to the elliptic integral K(k) as follows:

The elliptic integral E(k) is almost as easy, but we first need to define another term, cn:

c0=k    (the modulus)

cn=(an-1-bn-1)/2

Then take the sum

Then from King [6]:

Therefore

and also

Formula (8) may seem trivial since we already have K and E from formulae (5) and (7) and can simply subtract to get this difference, but the significance will become apparent as when we consider calculation precision and roundoff error.

Following is the source code for a Javascript function to calculate Nagaoka's coefficient using formula (4) and the elliptic integral formulae (5), (7) and (8):

function Nagaoka(u) {

// returns Nagaoka's coefficient,

// aka: field non-uniformity coefficient

// where u=diameter/length

if (u==0) return(1);

else {

  uu=u*u;

  m=uu/(1+uu);

  m2=4*Math.sqrt(1+uu);

  a=1;

  b= Math.sqrt(1-m);

  c=a-b;

  ci=1;

  cs=c*c/2+m;

  do {

    ao=(a+b)/2;

    b=Math.sqrt(a*b);

    a=ao;

    co=c;

    c=(a-b);

    cs+=ci*c*c;

    ci*=2;

  } while (c<co);

  cs/=2;

  K=Math.PI/(a+a);

  KmE=K*cs;

  E=K*(1-cs);

  return(1/(3*Math.PI)*(m2/uu*(KmE)+m2*E-4*u));

  }

}

In this code, the ratio of diameter to length is represented by the variable u. The square of the modulus is m. (The modulus k, is not used directly.) A few squares, square roots and constants cancel out, and have been eliminated for efficiency. The code has been optimized primarily for accuracy, and also for speed; not necessarily for compactness, though it is quite compact. The variables a and b are the same as described in the text above. They both converge to the AGM. Variable c is the same as cn described above, except without the divisor of 2. Variable cs is the sum cs. Variables a and b are initialized to 1 and k' respectively, as previously described, and then the arithmetic and geometric means are iterated in the loop until they converge to the same value. The resulting values of a, b and cs are then used to calculate K, E, and K-E (represented by variable KmE). Finally Nagaoka's coefficient kL is calculated and returned in the last line. In testing it was found that for u<10 it requires no more than 7 iterations. For 10<u<1000, no more than 9 iterations are required. 

Now to explain a couple of subtleties.

When u is small, the modulus k is small, and both K and E approach the same value π/2. Since formula (4) subtracts E from K, round-off problems arise, and can result in lost precision. In fact, during testing it was found that as u gets smaller and smaller, the result of the function correctly approaches the value 1, then gets a bit noisy in the last couple of decimal places, then suddenly jumps to a value of 2/3. This appeared to be largely due to the subtraction of K-E when u is small. Fortunately, we have formula (8) which calculates the principal terms for K-E directly, minimizing round-off error. Hence, the variable KmE is used in formula (4) in place of the subtraction K-E. Miller [8] does not account for this in his derivation; he uses the separate values of K and E. I should point out though, the error that occurs when using K-E is generally not significant when using double precision and calculating inductance of coils with practical diameter/length ratios. However, if this function is implemented using single precision arithmetic, the error is certain to manifest itself at larger diameter/length ratios and could be a significant problem. In any event there is peace of mind to be had, knowing that the function will produce meaningful results over the whole continuous range of inputs.

Next, is the convergence test. If we were evaluating a typical infinite series where each term gets progressively smaller at a rate which accelerates with each iteration, we could set the convergence test almost arbitrarily, to see when the term is smaller than some insignificant value such as 10-50. However, in the case of the AGM, we don't have this situation. The two values a and b are converging on a common value, and we have to be aware of the precision with which they are represented. If the variables only have 15 digits of precision, we may find that setting the convergence test to exit on values 10-16 or smaller, may fail to terminate. This is because there is no guarantee that a and b will ever be equal in the least significant digit, again because of round-off error. However, a simple way to ensure that we are getting the highest accuracy, without the worry of the loop not terminating, is to monitor the value c=a-b, and terminate when c stops getting smaller. For this, the variable co is used to save the previous value of c. Consequently, in this routine the loop exits when the condition (c<co) is no longer true.

Now that we have coded a function for Nagaoka's coefficient, we now have an accurate and convenient formula for inductance:

Converted to SI units:

Where,

µ0 = 4π×10-7 (The definition of the permeability of free space.)
r and ℓ are given in meters, and then LS will be in Henries.

The same function coded in Open Office BASIC follows:

Function Nagaoka(ByVal u as Double) As Double

' returns Nagaoka's coefficient,

' aka: field non-uniformity coefficient

' where u=diameter/length

if u=0 then

  Nagaoka=1

else

  uu=u*u

  m=uu/(1+uu)

  m2=4*sqr(1+uu)

  a=1

  b=sqr(1-m)

  c=a-b

  ci=1

  cs=c*c/2+m

  do

    ao=(a+b)/2

    b=sqr(a*b)

    a=ao

    co=c

    c=a-b

    cs=cs+ci*c*c

    ci=2*ci

  loop while c<co

  cs=(cs)/2

  K=pi()/(a+a)

  KmE=K*cs

  E=K*(1-cs)

  Nagaoka=1/(3*pi())*(m2/uu*KmE+m2*E-4*u)

end if

end Function

This can be pasted into Open Office's macro editor, and can then be used just like any other spreadsheet function. The only thing to be aware of is that when Open Office encounters a name in a formula, it first checks built-in functions, then named cells, and finally macros. So, if you have a cell containing the same text as the name of the macro, it will be fail to find the macro.

Results

Both the BASIC and Javascript versions of the function were tested and compared with Nagaoka's six decimal place tabulated results. The results agreed ±1 in the last decimal place (ie., 1 part per million) for u<10. For u=10, which is the highest value that Nagaoka calculated, his value is 9×10-6 lower than that produced by this function. However, this is in a region of u values where round-off error should not be significant, and furthermore, the function agrees to the last decimal place with Lundin's Handbook Formula [10] for this and larger values of u. Lundin's formula is asymptotically correct for large values. Therefore, it appears that the error is likely in Nagaoka's value, due possibly to limited precision in his calculation, accumulated error, or human error. As a further check on round-off error for small input values, I calculated the differences between values and plotted them. The resulting curve was very smooth, so I concluded that round-off is not a problem. While this is not conclusive, I would expect that if the calculation has been compromised due to round-off error, the difference curve would be noisy.

Update on accuracy (July 2012)

After this work was completed in 2010, the Nagaoka() function has been used, in a few cases, as a standard against which some empirical formulae have been compared. Because of this, and because of some discovered anomalies for extremely short coils, it became apparent that some further testing of the Nagaoka() function was called for, to determine the limits of applicability. The earlier statement, that the function is correct over the range which Nagaoka calculated his data, still stands. However, this function has recently been used with input arguments several orders of magnitude beyond this range, both larger and smaller.

Because the theory upon which the Nagaoka() function is based, contains no approximations, the only errors which may be introduced, are due to the finite numeric precision of the computer performing the calculations. This is generally referred to as round-off error.

In testing the Nagaoka() function, double precision arithmetic operations were used. This is approximately 15 decimal digits of precision. Two of the most common ways for significant round-off error to occur are:

1. •An extremely large number of mathematical operations being performed on the data such that the small round-off error that occurs for each operation eventually accumulates into a large error.
   

2. •Two nearly equal terms, in a formula, are subtracted which significantly amplifies the existing round-off error.
   

The first situation is rarely a problem unless a calculation involves millions of mathematical operations. Even then, the error from a mathematical operation may be either positive or negative with equal probability, and the errors tend to cancel rather than accumulate.

The second situation is of more concern here. Although it was stated earlier that the direct calculation of K-E eliminates a significant source of round-off error, the AGM method still includes some subtraction operations. The most problematic is the calculation of the square of the modulus, m, and the subsequent calculation of the complementary modulus. For large input arguments (extremely short coils), the value of

m=u^2/(1+u^2)

approaches one, which is then subtracted from one to calculate the complementary modulus. Thus, two nearly equal values are subtracted, greatly amplifying round-off error before the AGM calculation loop even begins.

Inside the loop, variable c is calculated by subtracting the arithmetic mean a from the geometric mean b. However, this subtraction is done merely to determine when a=b within the available precision of the computer arithmetic routines, and then terminate the loop. Hence, it should not introduce any round-off error.

After the loop terminates, there is one more subtraction that needs to be examined. That is in the calculation of

E=K*(1-cs)

In examining the values of cs, the closest that it ever gets to 1, is about 0.95, which occurs for extremely short coils (u>107), but this is not enough introduce any significant round-off error.

In the case of long coils, when calculating the complementary modulus, the subtraction will involve operands of significantly different values, and so round-off errors are not likely to result.

Having examined the function to see where we can expect round-off error to occur, how do we know exactly when the error becomes so great that we can no longer trust the output data?

Two possibilities are:

1. •Compare the output of the function with numbers from another source which are known to be accurate.
   

2. •Look at the smoothness of the output data.
   

In the first case we need to find a source of known accurate data. All of these calculations involve elliptic integrals in one way or another, and therefore, are subject to the same problems. However, series expansions of the elliptic integral based formulae have been done, and using the principal low order terms we can use these to set some asymptotes. For long coils, the asymptote is simply the value 1. For very short coils, it can be shown that kL is approximately given by the following expression:

kL=2/(πu)*[Ln(8/π+4u)-1/2]

Since these are asymptotes, there won't be exact agreement between these and the Nagaoka() function, but the difference should be very small, and that difference should be smoothly decreasing as as the input argument decreases (or increases as the case may be). Though, not a perfect comparison, this test can be combined with the smootness test.

In the case where we simply look at the smoothness of the output data, we are relying on the fact that round-off error looks like a noisy signal, and one way to accentuate noise is to differentiate the signal. Or in terms of finite data points, we simply take the differences between successive points, plot them, and look at how smooth the resulting curve is.

Long Coils:

Using separate E-K causes error for u< 3e-5 fluctuating slightly above and below 1.5e-8 it jumps to 0.6666667

For AGM method, the graph of differences shows smooth behaviour down to 1e-16 where value goes from .999999999999999 to 1

Short Coils

Both E-K and AGM fail at around 9e7 but start to get noisy at around 1e5.

Conclusion

Before moving on, it's important to reiterate: This formula assumes that the coil is wound, not with round wire, but with an infinitesimally thin conducting tape with essentially no gap between turns. A coil wound with round wire will have a slightly different inductance than what this formula produces. There are accurate corrections to the current sheet formula to account for round wire. The calculation of these are discussed in Section 2a.

Continue to Part 1b

Back to:

Numerical Methods Introduction

Radio Theory

Home

This page last updated: February 25, 2021

Copyright 2010, 2015 Robert Weaver