Numerical Integration for Magnetic Field of a Loop of Wire

Rapier
Messages
83
Reaction score
0

Homework Statement


Calculate the magnetic field of a current loop. Compare your numerical results with exact solution above the center of the loop. Investigate the effect of the grid size based on this comparison.

Homework Equations


dB = u0*I/4pi * (dL * R) / (R^2 + Z^2)^3/2
Bz = u0*I*R^2/ (2 * (R^2 + Z^2)^3/2 )

The Attempt at a Solution


I've tried to paste this into my numerical integration program from a couple of assignments ago (as per instructions). I'm using the trapezoid rule to do my integration. This is my code:

Code:
// Initialise Values
            a = 0;
            b = 2 * pi;
            dl = (b - a) / n;
            dz = 1 / n;
            B[0] = 0; // +c
            dB[0] = (u0 * current * radius * dl) / (4 * pi *pow(radius, 3));

            // Integration
            for (int loop = 0; loop < n; loop++)
                {
                    if ((loop == 0) || (loop == n )) weight = 0.5;
                        else weight = 1.0;
                       
                    z[loop + 1] = z[loop] + dz;
                    dB[loop + 1] = (u0 * current * radius * dl) / (4 * pi * pow((pow(radius, 2) + pow(x, 2)), 1.5));
                    B[loop + 1] = B[loop] + weight * 2 * pi *(dB[loop + 1] + dB[loop]) / 2 * dl;
       
                    printf("Z = %.2f", z[loop + 1]) << printf(", dB = %.6f", dB[loop + 1]) << printf(", B = %.6f\n", B[loop + 1]);
                }

The intention is to draw this diagram:
BField.jpg


But there is something not right. This is the data I am getting:

Z = 0.05, dB = 0.008944, B = 0.053762
Z = 0.10, dB = 0.008944, B = 0.071417
Z = 0.15, dB = 0.008944, B = 0.089072
Z = 0.20, dB = 0.008944, B = 0.106728
Z = 0.25, dB = 0.008944, B = 0.124383
Z = 0.30, dB = 0.008944, B = 0.142038
Z = 0.35, dB = 0.008944, B = 0.159693
Z = 0.40, dB = 0.008944, B = 0.177349
Z = 0.45, dB = 0.008944, B = 0.195004
Z = 0.50, dB = 0.008944, B = 0.212659
Z = 0.55, dB = 0.008944, B = 0.230314
Z = 0.60, dB = 0.008944, B = 0.247970
Z = 0.65, dB = 0.008944, B = 0.265625
Z = 0.70, dB = 0.008944, B = 0.283280
Z = 0.75, dB = 0.008944, B = 0.300935
Z = 0.80, dB = 0.008944, B = 0.318591
Z = 0.85, dB = 0.008944, B = 0.336246
Z = 0.90, dB = 0.008944, B = 0.353901
Z = 0.95, dB = 0.008944, B = 0.371556
Z = 1.00, dB = 0.008944, B = 0.389212

This is a rather long ongoing saga of an incompetent instructor. I've been pounding my head against this particular wall for over a week now and I'm just stuck. Any help would be greatly appreciated.
 
Physics news on Phys.org
Check your first equation against e.g. 9.1.11 here.
##\vec B## is a vector, so ##\vec {dB}## is a vector. The factor ##(dl * R)## should be ##\vec {dl}\times \vec r##

Note that you are not integrating over z (z is a constant). Only over ##\vec {dl}##.

And, unlike this 9.1.11 you cannot use the simplification for ##\vec r = \vec r_p - \vec {r'}## that they use: you want ##\vec B## for the general ##\vec r_p##, not just for P on the z axis.
 
Last edited:
For starters, thanks for your response. My Z axis is from the center of the loop straight up so because of symmetry I only have a z component for the magnetic field at point X (which is on the z axis).

I did realize you were correct in that I needed another radius term. Maybe it's just because I've been frustrated and making tiny adjustments constantly for the past week. I just don't get it. I have looked over and over the problem and the equations and I think this should be working:

Code:
// Initialise Values
            a = 0;
            b = 2 * pi;
            dl = (b - a) / n;
            dz = 1 / n;
            B[0] = 0; // +c
            dB[0] = (u0 * current * radius * dl) / (4 * pi *pow(radius, 3));

            // Integration
            for (int loop = 0; loop < n; loop++)
                {
                    if ((loop == 0) || (loop == n )) weight = 0.5;
                        else weight = 1.0;
                       
                    z[loop + 1] = z[loop] + dz;
                    dB[loop + 1] = (u0 * current * pow(radius,2) * dl) / (4 * pi * pow((pow(radius, 2) + pow(x, 2)), 1.5));
                    B[loop + 1] = B[loop] + weight * (dB[loop + 1] + dB[loop]) / 2 * dz;
       
                    printf("Z = %.2f", z[loop + 1]) << printf(", dB = %.6f", dB[loop + 1]) << printf(", B = %.6f\n", B[loop + 1]);
                }

I've tried plotting them to see kind of where I am and this is what I'm getting:
c1.jpg
. And the red line isn't even close the blue line like it should be.
 
I read the problem statement differently: they want you to calculate generally. The restriction to points on the z-axis is only for the comparison, because there you have an analytical solution to compare with. Link section 9.8 can help with the general case.

Again: z is a constant. It shouldn't change inside the integration loop. You only integrate over ##dl = R d\phi##.

You also confuse me by using the name dl for dφ . It's not illegal and you do treat it properly (except in B[0] !), but for me ##dl = R d\phi## is easier on the eyes.

Your weight thing should disappear. (you have loop < n, so loop == n doesn't happen -- right ?)

[edit] sorry, it comes from the trapezoid rule. But netto: just adding the dB should be the same thing, isn't
it ?

[edit2] more sorry. z is a contant in your calculation. Only you call it x. The printing during the loop isn't relevant: it shows how the result builds up. All that counts is the result. Vary x to compare with the analytical expression 9.1.15 (where it's called z).
 
Last edited:
z is in the loop not for the purpose of integration, but so that I have an array for the axis range. This projects tend to get messy as I started by adapting the trapezoidal integration code (which integrated a sin function) to calculate the magnetic field over a straight wire. Then I adapted that code to try and form the wire into a loop. In this particular case, the timing forced me to abandon the straight wire code before it was finished. So, the long and short of it is that I had some legacy variable names. You are quite correct in it being confusing, so I've changed it to dtheta. I've also changed the loop == (n-1) on the last spot for the trapezoidal rule. It was originally n-1 but that was one of the desperate changes I made trying to get it to work.
 
Rapier said:

Homework Statement


Calculate the magnetic field of a current loop. Compare your numerical results with exact solution above the center of the loop. Investigate the effect of the grid size based on this comparison.

Homework Equations


dB = u0*I/4pi * (dL * R) / (R^2 + Z^2)^3/2
Bz = u0*I*R^2/ (2 * (R^2 + Z^2)^3/2 )

The Attempt at a Solution


I've tried to paste this into my numerical integration program from a couple of assignments ago (as per instructions). I'm using the trapezoid rule to do my integration. This is my code:

Code:
// Initialise Values
            a = 0;
            b = 2 * pi;
            dl = (b - a) / n;
            dz = 1 / n;
            B[0] = 0; // +c
            dB[0] = (u0 * current * radius * dl) / (4 * pi *pow(radius, 3));

            // Integration
            for (int loop = 0; loop < n; loop++)
                {
                    if ((loop == 0) || (loop == n )) weight = 0.5;
                        else weight = 1.0;
                      
                    z[loop + 1] = z[loop] + dz;
                    dB[loop + 1] = (u0 * current * radius * dl) / (4 * pi * pow((pow(radius, 2) + pow(x, 2)), 1.5));
                    B[loop + 1] = B[loop] + weight * 2 * pi *(dB[loop + 1] + dB[loop]) / 2 * dl;
      
                    printf("Z = %.2f", z[loop + 1]) << printf(", dB = %.6f", dB[loop + 1]) << printf(", B = %.6f\n", B[loop + 1]);
                }

The intention is to draw this diagram: View attachment 80514

But there is something not right. This is the data I am getting:

Z = 0.05, dB = 0.008944, B = 0.053762
Z = 0.10, dB = 0.008944, B = 0.071417
Z = 0.15, dB = 0.008944, B = 0.089072
Z = 0.20, dB = 0.008944, B = 0.106728
Z = 0.25, dB = 0.008944, B = 0.124383
Z = 0.30, dB = 0.008944, B = 0.142038
Z = 0.35, dB = 0.008944, B = 0.159693
Z = 0.40, dB = 0.008944, B = 0.177349
Z = 0.45, dB = 0.008944, B = 0.195004
Z = 0.50, dB = 0.008944, B = 0.212659
Z = 0.55, dB = 0.008944, B = 0.230314
Z = 0.60, dB = 0.008944, B = 0.247970
Z = 0.65, dB = 0.008944, B = 0.265625
Z = 0.70, dB = 0.008944, B = 0.283280
Z = 0.75, dB = 0.008944, B = 0.300935
Z = 0.80, dB = 0.008944, B = 0.318591
Z = 0.85, dB = 0.008944, B = 0.336246
Z = 0.90, dB = 0.008944, B = 0.353901
Z = 0.95, dB = 0.008944, B = 0.371556
Z = 1.00, dB = 0.008944, B = 0.389212

This is a rather long ongoing saga of an incompetent instructor. I've been pounding my head against this particular wall for over a week now and I'm just stuck. Any help would be greatly appreciated.

It's not clear from this code snippet where the value of 'n' comes from.

Writing a program to do trapezoidal integration is a bit of overkill. A spreadsheet should be able to calculate the value of the integral using the trapezoidal rule quite nicely, IMO. :wink:
 
That would be nice, unfortunately it's a computational physics class and writing the code is kind of the entire point. 8P
 
Rapier said:
That would be nice, unfortunately it's a computational physics class and writing the code is kind of the entire point. 8P
Still, you want to check the results from the program. :wink:

The source of the value of 'n' is still not known.
 
Good. Now, to reconstruct the comparison plot: what are the values you use ? It is improbable that B is 1 Tesla at z = 0

SteamKing: n is the number of steps. "big enough". Doesn't matter anyway: each step has same dB ... :smile:
 
  • #10
n is the number of divisions for the trapezoidal rule and in this case, 20.
 
  • #11
BvU said:
Good. Now, to reconstruct the comparison plot: what are the values you use ? It is improbable that B is 1 Tesla at z = 0

SteamKing: n is the number of steps. "big enough". Doesn't matter anyway: each step has same dB ... :smile:

Yes, I know n is the number of steps. However, n is not defined in the code snippet, so how does the loop terminate?
 
  • #12
Rapier said:
n is the number of divisions for the trapezoidal rule and in this case, 20.

Is that value of n in the code snipped you attached to the OP? Is it user input?
 
  • #13
dl appears twice in the integrand and as integration variable. I don't trust the 2 pi there either !


dB[loop + 1] = (u0 * current * radius * dl) / (4 * pi * pow((pow(radius, 2) + pow(x, 2)), 1.5));
B[loop + 1] = B[loop] + weight * 2 * pi *(dB[loop + 1] + dB[loop]) / 2 * dl;

You want to remove the red stuff...

[edit] sorry, copied from the quote of post #1.

dB[loop + 1] = (u0 * current * pow(radius,2) * dl) / (4 * pi * pow((pow(radius, 2) + pow(x, 2)), 1.5));
B[loop + 1] = B[loop] + weight * (dB[loop + 1] + dB[loop]) / 2 * dz;

So only the dz is suspicious. Factor 20 nevertheless.
 
Last edited:
  • #14
z is in the loop not for the purpose of integration, but so that I have an array for the axis range
I'm still afraid that what you plot is not what you think you plot. I also don't understand it. If your table has constant dB, then how can the green line drop like a brick ?

Again: during the integration loop you don't want to save or plot anything. You want to vary x/radius and plot B(x/radius) / B(x=0)
 
Last edited:
Back
Top