Mathcad to calculate magnetization

[ptc-5691311]

Hi, everyone!

I am new here and I hope I've found the right place to ask my question.

I've been trying to simulate the magnetization process of a ferrimagnet using Mathcad. The theory behind this is quite simple - one writes down an energy expansion, there're three terms that signify how two magnetic moments interact with each and with magnetic field. Mathematically, one needs to minimize this energy with respect to the angles the magnetic moments make with the field and some symmetry direction of the crystal lattice. The minimization should be done for each magnetic field value, e.g., B = 0,1,2,...,20 T. Then, the total magnetization, Mtotal, can be calculated using the obtaines values of the angles.

I've prepared a Mathcad document for the simulation (see attachment). The program works but gives me an incorrect unphysical result. Could you take a look and possibly exlpain to me what I've done wrong?

Thank you in advance.

Regards,

Den

[Werner_E]

I don't see any result at all in the sheet you have posted and also I don't see any attempt to minize a function. You have defined some range variables which are not used at all (resp. have no efffect) and I guess the function dE() you defined isn't working as you expected. It will only return the derivation of E after Phi for B=60 in every case - the range has no effect. What should dE() return in your opinion?

And which expression is it you want to minimize?

[ptc-5691311]

Yep, I didn't include any result because it takes too long for the program to arrive at it. Now I've included the plot I'm after - Mtotal(B). As you can see, it makes no sense.

In my opinion, dE() should return the energy minimized with respect to alpha and phi. I want to minimize the initial expression, that is, E(B,alpha,phi). I thought I minimized the E function by taking its derivatives and then finding their roots to be used for Mtotal... Should I have assumed a different approach? And how should I get the B range to have the desired effect?

Thank you.

[Werner_E]

Denis Gorbunov wrote:

Yep, I didn't include any result because it takes too long for the program to arrive at it. Now I've included the plot I'm after - Mtotal(B). As you can see, it makes no sense.

Thats no surprise given the way the functions are defined.

In my opinion, dE() should return the energy minimized with respect to alpha and phi.

In what form? a vector? a 61x2 matrix ? What did you had in mind?

I want to minimize the initial expression, that is, E(B,alpha,phi).

Its a function in 3 variables, so you would need the derivatives wrt all 3. Also, how yould you cope with multiple (local) minima? And if you want the calculation being done for a constant B, I guess both partial derivates should be zero at the same time. I would suggest using either a solve block with Find() (using the two partial derivative) or one with Minimize().

I also don't understand why your function Mtotal should still be dependent on alpha and Phi and why you had setup ranges for those angles. f() and g() should already return the "optimal" values for those angles, as I think. Also f() and g() would not be dependent on alpha and Phi.

[ptc-5691311]

> In what form? a vector? a 61x2 matrix ? What did you had in mind?

Exactly, a 61x2 matrix: field vs. angle.

I haven't yet thought of multiple (local) minima. Yes, both partial derivatives should be zero. Thank you for your suggestions, I'm going to try them.

The function Mtotal is the total magnetization that rotates under the action of applied magnetic field. How it rotates depends on the angles alpha and phi (unknowns), on the strength of the applied field and on the absolute value of the magnetic moments (all known). The ranges of the angles were chosen according to the geometry of the problem.

[AlanStevens · Accepted solution]

Is the attached any help? I'm not sure I've interpreted what you want correctly, but if I have, you need to note that the results are somewhat sensitive to the initial guesses.

Alan

[ptc-5691311]

Alan,

Thank you for your reply. However, I can't open the attachment in my version of MathCad (2001 Professional). Could you possibly save it in this, I guess, older version?

Regards,

Den

[AlanStevens]

Denis Gorbunov wrote:

Alan,

Thank you for your reply. However, I can't open the attachment in my version of MathCad (2001 Professional). Could you possibly save it in this, I guess, older version?

Regards,

Den

Ok. Here it is saved from Mathcad 11 (which is the oldest version I have) in 2001 and 2001i formats (I don't know if these are significantly different). I'm unable to check if it works in those formats though.

Alan

[Werner_E]

Compared to the original sheet for the calculation of Mtotal the sine of the sum of the angles alpha and Phi is used, mot the sum of the sines.

[AlanStevens]

Werner Exinger wrote:

Compared to the original sheet for the calculation of Mtotal the sine of the sum of the angles alpha and Phi is used, mot the sum of the sines.

Yes, a silly slip on my part!

Alan

[ptc-5691311]

Dear Alan,

I succeeded in opening your Mathcad file. Thank you!

I looked through your file and in the last expression, that for Mtotal, I changed (sin(alpha_i) + sin(phi_i) to (sin(alpha_i+phi_i)) and immediately got the expected result. It is correct in a sense that the low-field magnetization is (MEr-MCo) = 9 - 5 = 4 (magnetic moments are antiparallel), and the high field magnetization is (MEr + MCo) = 9 + 5 = 14 (magnetic moments are parallel). Moreover, the overall shape of the curve is correct: linear section -> slanted section -> linear section.

I've attached your file with the modified formula for Mtotal and an increased range of B values, 0...200, just to see the whole magnetization process.

From your file it is clear to me that I should learn to use solve blocks to do this sort of calculations. Well, at least I know what my problem is.

Thank you.

Regards,

Den

[AlanStevens]

1.

Denis Gorbunov wrote:

... It is correct in a sense that the low-field magnetization is (MEr-MCo) = 9 - 5 = 4 (magnetic moments are antiparallel), and the high field magnetization is (MEr + MCo) = 9 + 5 = 14 (magnetic moments are parallel).

Doesn't this mean you are dealing with an anti-ferrimagnet (as opposed to a ferrimagnet)?!!

2. I was puzzled by why the B=0 value of magnetization found by Mathcad is 0, not 4. It seems this is a numerical problem. If you set the initial guess value of phi to -pi/2 you do get 4 for the magnetization at B = 0.

Alan

[ptc-5691311]

Doesn't this mean you are dealing with an anti-ferrimagnet (as opposed to a ferrimagnet)?!!

No both, antiferromagnets and ferrimagnets, consist of at least two magnetic sublattices whose magnetic moments are antiparallel to each other. The difference between the two is that in the former the moments are equal in absolute value, whereas in the latter they are now. In our case, MEr = 9 and MCo = 5 so we have a ferrimagnet

2. I was puzzled by why the B=0 value of magnetization found by Mathcad is 0, not 4. It seems this is a numerical problem. If you set the initial guess value of phi to -pi/2 you do get 4 for the magnetization at B = 0.

This is exactly what I thought too. Luckily the problem is easily amended.

[AlanStevens]

Denis Gorbunov wrote:

2. I was puzzled by why the B=0 value of magnetization found by Mathcad is 0, not 4. It seems this is a numerical problem. If you set the initial guess value of phi to -pi/2 you do get 4 for the magnetization at B = 0.

This is exactly what I thought too. Luckily the problem is easily amended.

Thinking about it more carefully, I realise that when B is zero, the minimisation is independent of phi. This means phi will just retain its initial value.

Also, I notice that, for positive B (with initial guess for alpha of pi/2 and phi of -pi/2) the transition changes in type from the gradual increase you had previously to a step change at B=90. Which is the correct transition?

Also, shouldn't alpha and phi be 180deg out of phase for small B? Mathcad seems to insist on their being 90deg out of phase! Though I'm assumiong here that they are measured from the same origin, which might not be true.

Alan

[ptc-5691311]

You led me to the same conclusion: when B = 0, phi should stay zero.

I checked the situation for initial guess phi = -pi/2. Indeed, the type of the transition changes. The correct physical picture is when the transition is a monotonous increase in M rather than an abrupt jump. The reason is simple: the applied magnetic field competes against the exchange interaction between the two magnetic moments. They are antiaparallel initially, and they cannot abruptly assume the parallel orientation. So the competition results in a monotonous increase in M until finally the angle between the two moments is zero.

Now, about the angles... I was inattentive to them, unfortunately. As you pointed out, initially they should be 180 degrees out of phase, and for high fields the same should apply but with the opposite sign. Taking a closer look, I realize they are 90 degrees out of phase for both, low and high fields. I guess I must take a look at the initial picture of the problem. Despite the correct result!

Den

[ptc-5691311]

Here's a file where the angles are correct for low and high fields. I re-wrote the initial formula for E and changed alpha and phi intervals, in an attempt to simplify the initial model. As you pointed out before, the result is sensitive to the initial guess. Although now the result is correct, I guess it is more desirable to find a way to calculate M(H) without making an initial guess. I mean, the physics should be right regardless of the method used for the calculation, right?

[Werner_E]

I forgot to send my attempt yesterday. Basically it seems I had interpreted your attempts similary as Alan.

I have shown two different ways to get the minimum of E(). The firstis using the partial derivatives as you had tried but seems not appropriate given the shape of your function. I added a 3D-plot of the function we wish to minimize for B=7 and you can see that we have no useable minimum where both derivatives are zero. So I also passed the addition of a condition which would guarantee we get a minumum and not a maximum. Angles2() is pretty much the same as Alan had showed - the difference being that here the guess values are not constant but assigned as parameters of the function. Idea was that you could adopt those guesses dependent on B later. I didn't do that as I guess with the function you provided the minimum will always be at alpha=pi and Phi=pi/2. Maybe you should check your formulas.

Unfortunalety current versions of Mathcad (14/15) can only save down to MC11 format. We would need someone who has installed MC11 as well to save in the needed 2001 format.

So I attach my file and also pdf prints of Alan's and my file for your information.

EDIT: Had this reply open and worked intermediate on something, so I am late as Alan already sent his file in 2001 format. I'll delete the pdf of his file.

[Werner_E]

I had another look at my file and I noticed that Alan was right in that the parameter B, which is not optimized by minimize(), should be the last of the three, not the first.

I also added I third way to determine the minimum of E() (you already know that you can forget about the first one with the derivatives). Its a simple brute force algo which looks at each value in a grid of custom size and returns the angles which gives the minimum E-value. The results are (refrained from the first few values of B) pretty much the same as in Alan's sheet or with the second method in my sheet. We need a rather fine gridsize (grid>100) to avoid steps in the graph of Mtotal for higher B-values.

Anothe thing to mention: If we use minimize() we have to be very careful with constrains as those will give us wrong results with bad guess values.

I tend to call it a limitation of the numeric algoriothm rather than a bug.

[ptc-5691311]

Dear Werner,

Thank you for your Mathcad and .pdf files. I've made sure that using partial derivatives is not appropriate to my problem. The proof is that at some field values the magnetization decreases, whereas it should either increase or stay constant. However, when you used new guess values, the result you showed in the penultimate graph of Magnetization2_Werner file is definitely interesting. This is not what I expected but it is possible from the physics point of view. As I referred to Alan's example in the previous post, from this graph you get the parallel orientation of the magnetic moments in the lowest applied field.

In your example (magnetization3) where you use the same guess values for alpha and phi as Alan, you get the same expected result. Your third example where you use a grid algorithm yields the same result but seems less clear to me as I don't know how it works. But I am going to look into it.

Thank you for your help!

Regards,

Den

[Werner_E]

Your third example where you use a grid algorithm yields the same result but seems less clear to me as I don't know how it works. But I am going to look into it.

As its bruteforce its nothing special or complicated. E() is evaluated for a constant B (the argument of the function) for "all" values of alpha and Phi in the specified range (0..pi, -pi/2 .. pi/2) and the pair which yields the smallest E-value wins. This is done by two nested for-loops, but of course we cannot try ALL values so a grid of angle values is laid out, the denseness of that grid is controlled by the worksheet variable grid. If the grid is too coarse we may miss the real minimum and this seems to happen very easily for the higher values of B.

Find attached just for fun an animation showing the surface showing E depending on alpha and Phi as well as the minimum.

[ptc-5691311]

Thanks for the animation! At B = 40 where the dM/dH slope changes the point jumps away and starts moving! Awesome!

From your explanation it follows that a grid algorithm is really not complicated. Thank you for the explanation!

[Werner_E]

From your explanation it follows that a grid algorithm is really not complicated.

No, it isn't. But the drawback, as usual for brute force, is that its slow as of the many calculations necessary. A grid of 150x150 isn't that dense, but already requires 22500 times calculating E and comparing the result to the current minimum and replacing the latter if appropriate. And in graphing MTotal this procedure has to be done for every value of B which is plotted. So the conjugate gradient algorithm used by minimize() still is faster and given suitable guess values is delivering satisfying results. You may play around with the algorithms offered and some options by right clicking the word minimize (not sure if thats true for MC2001i, toom though) but usually auto select is doing a good job.

[ptc-5691311]

Thank you, I'll give it a try!

[ptc-5691311]

Dear Werner and Alan,

I've found a Mathcad file of one of my colleagues who I am no longer in touch with. The program was designed to treat pretty much the same problem - to calculate magnetization curves of ferrimagnets. It suits our case very well. However, it is much more complicated than what I am able to come up with so far. You might have noticed I can write just very simple programs. Could you please explain to me what individual steps in the program mean? I really don't get almost anything. The file is attached.

Thank you in advance.

Den

[Werner_E]

Any specific question?

I just looked at the routine AbsMin2(A) and it is really very tricky. It could easily be replaced by

but maybe match() was not available in MC2001.

The routine determines the minumum value of the matrix via the built-in function min() and then searches for the indices of that value. In the original routine this is done quite wily using just one sequential variable d (many of us would have used a pair of nested for loops, I guess).

Later a matrix of energy values is created with a stepwidth of 1deg for both angles and then the aformentioned routine used to determine the pair of angels yielding the minimum energy - kind of brute force. Its noticeable that the range for both angles is different from yours (180°..360° & 90°..270° vs yours 0°..180° & -90°..+90°).

Using the angles as matrix indices is a waste of memory and not as it should be done. I guess it was done for convenience so the scale in the countour plot is shown correct. Otherwise you would have to create a much more complex data structure for the contourplot. But then I may be wrong as I remember a thread where it was shown that there are differences in the behaviour of contour plots between the current version of MC and older versions.

[ptc-5691311]

Hehe I have no specific questions since I don't understand most of the algorithms used in the program I have to get a Mathcad manual and read carefully first.

Werner, I thank you for your helpful explanation. It'll be even more helpful when I understand the program better!

Den

P.S. As to the angles, this problem was designed for ferrimagnets with a different magnetic symmetry so the angles need be different too.

[AlanStevens]

I've attached an annotated copy of the Jump worksheet (in 2001 format). The comments relate mainly to what I think each routine is trying to do rather than to explain in detail the meaning of each Mathcad command.

I've also reduced the resolution of B in the magnetisation calculations to speed up the routine.

Hope this helps.

Alan

[ptc-5691311]

Alan,

Thank you for your comments! Of course they help! I'll try to digest them all.

Den