MATLAB Examples

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