MATLAB Examples

Contents

Cubic anisotropy example

We will find the equilibrium states for a single-domain magnet with cubic anisotropy in zero field. The magnetization has a fixed magnitude, so it can be represented by a unit vector $ (m_1,m_2,m_3) $. The normalized energy of the magnet is

$$ E = K(1-m_3)^2 + m_1^2m_2^2+m_2^2m_3^2+m_3^2m_1^2. $$

We use a Lagrange multiplier to ensure that the magnetization remains a unit vector. Defining

$$ E' = E + \lambda(m_1^2+m_2^2+m_3^2-1), $$

In equilibrium, the derivatives of $E$ with respect to $m_1,m_2,m_3$ and $\lambda$ are all zero.

m = sym('m',[3 1]);
syms K lambda

E = K*(1-m(3)^2) + m(1)^2*m(2)^2+m(2)^2*m(3)^2+m(3)^2*m(1)^2  + lambda*(m(1)^2+m(2)^2+m(3)^2-1);
f = [diff(E,m(1)); diff(E,m(2)); diff(E,m(3)); diff(E,lambda)];
disp(f)
          2*m1*m2^2 + 2*m1*m3^2 + 2*lambda*m1
          2*m2*m1^2 + 2*m2*m3^2 + 2*lambda*m2
 2*m3*m1^2 + 2*m3*m2^2 - 2*K*m3 + 2*lambda*m3
                       m1^2 + m2^2 + m3^2 - 1
 

First try

We will solve for a trial value of $K$.

Kval = 0.2;
poly_system = BertiniLab('variable_group',[m; lambda],'function_def',f, ...
    'constant',[K Kval]);
poly_system = solve(poly_system);

That run had a lot of truncated infinite paths:

results = poly_system.solve_summary;
istart = strfind(results,'Paths Tracked');
disp(results(istart:end))
Paths Tracked: 54
 Truncated infinite paths: 28 - try adjusting SecurityMaxNorm or set SecurityLevel to 1 in the input file
   Please see 'failed_paths' for more information about these paths.


Improved approach

We follow the advice of Bertini and set the security level to 1.

poly_system.config = struct('SecurityLevel',1);
poly_system = solve(poly_system);

results = poly_system.solve_summary;
istart = strfind(results,'Finite Solution Summary');
disp(results(istart:end))
Finite Solution Summary

NOTE: nonsingular vs singular is based on condition number and identical endpoints

		| Number of real solns	|  Number of non-real solns	|  Total
------------------------------------------------------------------------------------------
Non-singular	|	26		|		0		|   26
Singular	|	0		|		0		|   0
------------------------------------------------------------------------------------------
Total		|	26		|		0		|   26


Finite Multiplicity Summary

  Multiplicity	|  Number of real solns	|  Number of non-real solns
------------------------------------------------------------------------------------------
	1	|	26		|	0


Infinite Solution Summary

NOTE: nonsingular vs singular is based on condition number and identical endpoints

		| Number of real solns	|  Number of non-real solns	|  Total
------------------------------------------------------------------------------------------
Non-singular	|	0		|		0		|   0
Singular	|	0		|		1		|   1
------------------------------------------------------------------------------------------
Total		|	0		|		1		|   1


Infinite Multiplicity Summary

  Multiplicity	|  Number of real solns	|  Number of non-real solns
------------------------------------------------------------------------------------------
	28	|	0		|	1


------------------------------------------------------------------------------------------
The following files may be of interest to you:

main_data:             A human-readable version of the solutions - main output file.
raw_solutions:         A list of the solutions with the corresponding path numbers.
raw_data:              Similar to the previous, but with the points in Bertini's homogeneous
                         coordinates along with more information about the solutions.
real_finite_solutions: A list of all real finite solutions.
finite_solutions:      A list of all finite solutions.
nonsingular_solutions: A list of all nonsingular solutions.
singular_solutions:    A list of all singular solutions.
------------------------------------------------------------------------------------------

Paths Tracked: 54

Here are the solutions for the magnetization in double precision. The number of paths is equal to the Bezout number for the system ($3^3\cdot 2 = 54 $). They have reflection symmetries in the XY, YZ and XZ planes as well as invariance under swapping of the x and y coordinates.

sols = poly_system.match_solutions('finite_solutions',polysym(m));
mag = real(double(sols.m));
disp(mag.')
    0.6325    0.6325    0.4472
    0.6325    0.6325   -0.4472
    0.7071    0.7071    0.0000
    0.6325   -0.6325    0.4472
    0.7071   -0.7071    0.0000
    0.7746    0.0000    0.6325
    0.7746   -0.0000   -0.6325
    1.0000    0.0000    0.0000
    0.6325   -0.6325   -0.4472
   -0.6325    0.6325    0.4472
   -0.7071    0.7071    0.0000
   -0.7071   -0.7071    0.0000
   -0.7746    0.0000    0.6325
   -0.7746   -0.0000   -0.6325
   -1.0000   -0.0000   -0.0000
    0.0000    0.7746    0.6325
   -0.0000    0.7746   -0.6325
    0.0000    1.0000    0.0000
   -0.6325    0.6325   -0.4472
   -0.0000   -0.7746    0.6325
   -0.0000   -0.7746   -0.6325
   -0.0000   -1.0000   -0.0000
    0.0000    0.0000    1.0000
   -0.6325   -0.6325    0.4472
   -0.0000   -0.0000   -1.0000
   -0.6325   -0.6325   -0.4472

Parameter homotopy: ab initio run

Now we will use parameter homotopy to obtain solutions for several values of H. First comes the ab initio run:

poly_system.config = struct('ParameterHomotopy',1);
poly_system.constant = polysym.empty(0,2);
poly_system.parameter = K;
poly_system = solve(poly_system,'paramotopy.input');

Parameter homotopy: using symmetries

Next the nonsingular solutions are used as starting points. First, however, we can use the symmetries to reduce the number of paths to 5:

ia = [];

% Sometimes this method obtains the wrong number of starting points.
while (numel(ia) ~= 5)
    vars = poly_system.read_solutions('nonsingular_solutions');
    v = double(vars);
    v = [sort(v(1:2,:)); v(3:4,:)];
    v = round(v*1e10);
    [~,ia] = unique(abs(v.'),'rows');
end
vars_reduced = vars(:,ia(:));

Set up poly_system for parameter homotopy:

poly_system.starting_points = vars_reduced;
poly_system.config.ParameterHomotopy = 2;

Parameter homotopy: an example run

Run again for a single value of K.

poly_system.final_parameters = Kval;
poly_system = solve(poly_system,'paramotopy.input');

results = poly_system.solve_summary;
istart = strfind(results,'Paths Tracked');
disp(results(istart:end))
Paths Tracked: 5

The number of paths tracked has been reduced from 54 to 5.

Now run for several values of K.

N = 200;
Kval = linspace(-1,1,N);
mag = zeros(3,5,N);
lam = zeros(1,5,N);
problemCases = struct.empty;
for ii=1:N
    poly_system.final_parameters = Kval(ii);

    poly_system = solve(poly_system,'paramotopy.input');
    sols = poly_system.match_solutions('raw_solutions',polysym(m),polysym(lambda));

    if poly_system.has_failed_paths || poly_system.has_path_crossings
        problemCases.number = ii;
        problemCases.diagnosis = poly_system.solve_summary;
    end

    mag(:,:,ii) = sols.m;
    lam(:,:,ii) = sols.lambda;
end
if ~isempty(problemCases)
    warning('There were some failed paths or path crossings. See problemCases for more information.')
end

Plot magnetization

map = pmkmp(N,'IsoL');
%map = gray(N);
for ii=1:N
    I = all(abs(imag(mag(:,:,ii)))<1e-12);
    mx = real(mag(1,I,ii)); my = real(mag(2,I,ii)); mz = real(mag(3,I,ii));
    m1 = [mx -mx mx -mx mx -mx mx -mx my -my my -my my -my my -my];
    m2 = [my my -my -my my my -my -my mx mx -mx -mx mx mx -mx -mx];
    m3 = [mz mz mz mz -mz -mz -mz -mz mz mz mz mz -mz -mz -mz -mz];
    line(m1,m2,m3,'Marker','.','LineStyle','none','Color',map(ii,:))
end
colorbar
caxis([min(Kval) max(Kval)])
colormap(map)
view(148,26)
axis off equal

% Add axes and an equator
line([0 0],[0 0],[0 1],'Color','k')
line([0 0],[0 1],[0 0],'Color','k')
line([0 1],[0 0],[0 0],'Color','k')
phi = linspace(0,2*pi,200);
x = cos(phi); y = sin(phi); z = zeros(size(phi));
line(x,y,z,'Color','k')

% Make the colorbar a bit smaller.
h_bar = findobj(gcf,'Tag','Colorbar');
initpos = get(h_bar,'Position');
set(h_bar, ...
   'Position',[initpos(1)+initpos(3)*0.15 initpos(2)+initpos(4)*0.15 ...
      initpos(3)*0.7 initpos(4)*0.7])

Witness sets

Now we treat K as another variable and see whether there are any higher-dimensional solutions.

poly_system = BertiniLab('variable_group',[m; lambda; K],'function_def',f, ...
    'config',struct('TrackType',1));
poly_system = solve(poly_system);

results = poly_system.solve_summary;
istart = strfind(results,'************** Decomposition');
disp(results(istart:end))


% %%
% % There are 10 degree 1 components (i.e., points that do not change as K
% % changes) and four degree 4 components (the curves in the figure). We will
% % see which curves go with which components by doing a membership test for
% % some points in the middle of the curves.
% v100 = [mag(:,:,100); lam(:,:,100)];
% idx = any(abs(v100.^2-1)<1e-10 | abs(v100.^2-0.5)<1e-10);
% v100 = v100(:,~idx);
%
% %%
% % Look for permutations of the first root.
% mx = v100(1,1); my = v100(2,1); mz = v100(3,1); l = lam(1,1); kv = Kval(100);
% m1 = [mx -mx mx -mx mx -mx mx -mx my -my my -my my -my my -my];
% m2 = [my my -my -my my my -my -my mx mx -mx -mx mx mx -mx -mx];
% m3 = [mz mz mz mz -mz -mz -mz -mz mz mz mz mz -mz -mz -mz -mz];
% v = [m1; m2; m3; l*ones(size(m3)); kv*ones(size(m3))];
%
% %%
% vd = round(v*1e10);
% [~,ia] = unique(vd.','rows');
% v = v(:,ia(:));
%
% %%
% poly_system.member_points = v;
% poly_system.config.TrackType = 3;
% poly_system = poly_system.solve;
************** Decomposition by Degree **************

Dimension 1: 14 classified components
-----------------------------------------------------
   degree 1: 10 components
   degree 4: 4 components

*****************************************************