Numerical Integration for Magnetic Field of a Loop of Wire

[Rapier]

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:

// 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:

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.

 

[BvU]

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.

 

[Rapier]

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:

// 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:

. And the red line isn't even close the blue line like it should be.

 

[BvU]

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).

 

[Rapier]

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.

 

[SteamKing]

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:

// 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.

 

[Rapier]

That would be nice, unfortunately it's a computational physics class and writing the code is kind of the entire point. 8P

 

[SteamKing]

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.

The source of the value of 'n' is still not known.

 

[BvU]

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

 

[Rapier]

n is the number of divisions for the trapezoidal rule and in this case, 20.

 

[SteamKing]

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

Yes, I know n is the number of steps. However, n is not defined in the code snippet, so how does the loop terminate?

 

[SteamKing]

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?

 

[BvU]

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.

 

[BvU]

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)