1. Not finding help here? Sign up for a free 30min tutor trial with Chegg Tutors
    Dismiss Notice
Dismiss Notice
Join Physics Forums Today!
The friendliest, high quality science and math community on the planet! Everyone who loves science is here!

Numerical Integration for Magnetic Field of a Loop of Wire

  1. Mar 17, 2015 #1
    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: 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.
     
  2. jcsd
  3. Mar 18, 2015 #2

    BvU

    User Avatar
    Science Advisor
    Homework Helper
    Gold Member

    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
  4. Mar 18, 2015 #3
    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: c1.jpg . And the red line isn't even close the blue line like it should be.
     
  5. Mar 18, 2015 #4

    BvU

    User Avatar
    Science Advisor
    Homework Helper
    Gold Member

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

    SteamKing

    User Avatar
    Staff Emeritus
    Science Advisor
    Homework Helper

    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:
     
  8. Mar 18, 2015 #7
    That would be nice, unfortunately it's a computational physics class and writing the code is kind of the entire point. 8P
     
  9. Mar 18, 2015 #8

    SteamKing

    User Avatar
    Staff Emeritus
    Science Advisor
    Homework Helper

    Still, you want to check the results from the program. :wink:

    The source of the value of 'n' is still not known.
     
  10. Mar 18, 2015 #9

    BvU

    User Avatar
    Science Advisor
    Homework Helper
    Gold Member

    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:
     
  11. Mar 18, 2015 #10
    n is the number of divisions for the trapezoidal rule and in this case, 20.
     
  12. Mar 18, 2015 #11

    SteamKing

    User Avatar
    Staff Emeritus
    Science Advisor
    Homework Helper

    Yes, I know n is the number of steps. However, n is not defined in the code snippet, so how does the loop terminate?
     
  13. Mar 18, 2015 #12

    SteamKing

    User Avatar
    Staff Emeritus
    Science Advisor
    Homework Helper

    Is that value of n in the code snipped you attached to the OP? Is it user input?
     
  14. Mar 18, 2015 #13

    BvU

    User Avatar
    Science Advisor
    Homework Helper
    Gold Member

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

    BvU

    User Avatar
    Science Advisor
    Homework Helper
    Gold Member

    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
Know someone interested in this topic? Share this thread via Reddit, Google+, Twitter, or Facebook

Have something to add?
Draft saved Draft deleted



Similar Discussions: Numerical Integration for Magnetic Field of a Loop of Wire
Loading...