Numerical Methods for
Inductance 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 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 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 kʹ
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 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
simply 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:
 |
(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−1
− bn−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 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.
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 = a−b,
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 K−E
when u is small. Fortunately, we have formula (8) which
calculates the value of K−E directly,
without subtraction, eliminating this source of error. Hence, the
variable KmE is used in formula (4) in place of the
subtraction K−E. 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
u < 10. For u = 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.
Continue to Part 1b
