# Homework Help: Numerical Integration for Magnetic Field of a Loop of Wire

1. Mar 17, 2015

### Rapier

1. The problem statement, all variables and given/known data
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.

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

3. 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 (Text):
// 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.

2. Mar 18, 2015

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

Last edited: Mar 18, 2015
3. Mar 18, 2015

### 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 realise 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 (Text):
// 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.

4. Mar 18, 2015

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

 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: Mar 18, 2015
5. Mar 18, 2015

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

6. Mar 18, 2015

### SteamKing

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

7. Mar 18, 2015

### Rapier

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

8. Mar 18, 2015

### SteamKing

Staff Emeritus
Still, you want to check the results from the program.

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

9. Mar 18, 2015

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

10. Mar 18, 2015

### Rapier

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

11. Mar 18, 2015

### SteamKing

Staff Emeritus
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. Mar 18, 2015

### SteamKing

Staff Emeritus
Is that value of n in the code snipped you attached to the OP? Is it user input?

13. Mar 18, 2015

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

 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: Mar 18, 2015
14. Mar 18, 2015

### BvU

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: Mar 18, 2015