# Numerical Methods forInductance Calculation

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

In many cases when trying to analyze a real physical system, it's often found to be far too complex to yield a simple mathematical model. We then attempt to change the model into something that will behave similarly but will be much easier to model. In the case of a single layer solenoid coil, we find that if we replace the actual winding with an infinitesimally thin conducting tape with negligible gap between turns, similar to the illustration, then the magnetic field that is produced is virtually the same as that due to the actual winding. Lorenz [22][1] studied an even more simplified model (just a thin walled conducting tube) where the current flowed circularly around the solenoid, but did not have any axial current component. He gave an exact formula for calculating the inductance of this solenoidal current sheet. His formula is exact in this case, and is a very good approximation for coils wound with round wire, as long as they have many turns and the windings are closely spaced.

Lorenz's formula is:

 (1)

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

And where: r is the radius, is the length, and N is 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 equal to 10−7 H/m. An important implication is that because µ0 is in units of henries/meter, it converts the result of the formula from units of length into the electrical unit of inductance. (Early scientific formulae expressed both capacitance and inductance in centimeters.) In SI units, the result will be in henries, when the radius and length are given in meters. Converted to modern SI units formula (1) becomes:

 (1a)

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

µ0 = 4π × 10−7 H/m (The definition of the permeability of free space.)

While we are defining the modulus k, we will also define the complementary modulus 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 K and E, 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:

 (2)

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:

 (3)

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:

 (4)

the modulus k and complementary modulus are now defined as:

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—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 an1 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 simply as a and b.

a0 = 1

and

b0 =     (the complementary modulus)

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

 (5)

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

c0 = k    (the modulus)

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

Then take the sum

 (6)

Then from King [6]:

Therefore

 (7)

and also

 (8)

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 if (u > 100000) return (2/(Math.PI*u)*(Math.log(8/Math.PI+4*u)-0.5)); 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)); } } ```

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 elseif u>100000 then Nagaoka=2/(pi()*u)*(Log(8/pi()+4*u)-0.5) 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 ```

The BASIC version 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.

In these routines, 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 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 KE (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.

The loop convergence test deserves some comment. 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 to 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 due to 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 = ab, and terminate when c stops getting smaller. For this, the variable c0 is used to save the previous value of c. Consequently, in this routine the loop exits when the condition (c < c0) is no longer true.

###### Accuracy of the Subroutines

For very long coils, the value of 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 can arise and result in lost precision. During testing, it was found that as u gets small, the result of the function correctly approaches the value 1, but then develops unacceptable round-off error as u becomes extremely small. This appeared to be largely due to the subtraction of KE when u is small. Fortunately, we have formula (8) which calculates the value of KE directly, without subtraction, eliminating this source of error. Hence, the variable KmE is used in formula (4) in place of the subtraction KE. This significantly improves the accuracy for values of u right down to zero (infinitely long coils). However, for the limiting case where u is exactly equal to zero, the algorithm cannot be used due to a divide by zero condition, and so the if{} statement near the beginning of the code simply returns the correct value of 1.

For extremely short coils, the value of u is very large, and this can also create problems due to the subtraction of two nearly equal values when calculating the complementary modulus m. However, using a series expansion of the elliptic integral formulae and taking just the low order terms, we can get an approximate formula for kL that becomes increasingly accurate as u becomes larger. Thus, for very short coils it can be shown that kL is approximately given by the following expression:

kL = 2/(πu)×(Ln(8/π + 4u) − ½)

This formula is asymptotically correct for increasing values of u. The point where the round-off error of the exact formula exceeds the approximation error of the approximate expression is in the vicinity of u = 100,000. Therefore, the second conditional clause in program code substitutes the approximate formula when u > 100,000.

###### Results

Both the BASIC and Javascript versions of the function were tested and compared with Nagaoka's six decimal place table. The results agree within ±1 in the last decimal place (ie., 1 part per million) for < 10. For = 10, which is the highest u value that Nagaoka calculated, his value of kL is 9×10−6 lower than that produced by this function. However, this is in a region of u values where the program's round-off error is negligible. In addition, the program 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 can be concluded that this single anomaly is an error in Nagaoka's value, due to limited precision in his calculation, accumulated error, or human error. We can conclude that the program output is accurate to within 1 part per million.

###### Practical Use

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

 (9)

Where,
µ0 = 4π × 10−7   H/m
Radius and length r and are given in meters, and LS is in henries.

###### Revision Note - 2022-11-28

This page first went online in 2010. It included no discussion of the accuracy of the computer routines. In 2012, to remedy this, I appended a rather longwinded and rambling accuracy discussion to the end of the page. In the process of revamping the complete website in 2022, this page has been re-edited, and the discussion of accuracy has been merged in with the main body of the text.

In 2019, the base units used in the SI system were redefined. Planck’s constant, and the elementary charge were set as exact defined quantities. As a result, µ0 (which was previously defined as exactly 4π×10-7 H/m) is now an experimentally determined constant with the value 1.25663706212×10−6 H/m which varies from the previous value by one in the 10th significant digit. Hence, we can no longer say that µ0 is exactly 4π×10-7 H/m. However, for the degree of accuracy that is needed for normal engineering calculations, 4π×10-7 is accurate enough and easier to remember than 1.25663706212×10−6 H/m.

Back to:
Numerical Methods Introduction