Numeric integration for systems of vector equations

I am trying to model the rotational motion of a system of particles under an applied magnetic field. To accomplish this, I have a system of differential equations that describes the motion of each particle based on the applied magnetic field and the magnetization of nearby particles. Each particle is represented by a unit vector for its direction, so for N particles, this is N vectors, and I use a 3 by N matrix to represent the system. Is there a more appropriate built-in MATLAB function that can handle the integration of this system? I attempted to use ODE45, but it is not giving the results that I would expect.

This is how my code is set up:

% Declare a number of particles
N = 27;
% This is just a constant used in later calculations
QQQ
% Create a set of random initial directions for the particles
D = rand(3, N);
% Convert the particles to unit vectors
D = D./(ones(3,1)*sum(D.^2).^(1/2));
% Initialize field in z-direction
H0x = 0;
H0y = 0;
H0z = 10^5;
% Make D into a column vector of initial directions for use in ODE45
init = (D(1:3*N))';
% Use ode45 to numerically integrate
[t, data] = ode45(@(t, params) odefcn(t, params, N, H0x, H0y, H0z, QQQ), [0 0.5], init);
% Reform D into a matrix of with N columns, each representing a direction vector
    for i = 1:3*N
        if rem(i - 1, 3) == 0
            D(1, floor(i/3) + 1) = data(i);
        elseif rem(i - 2, 3) == 0
            D(2, floor(i/3) + 1) = data(i);
        else
            D(3, floor(i/3)) = data(i);
        end
    end

And here is the function:
function [dydt] = odefcn(t, params, N, H0x, H0y, H0z, QQQ)
D = zeros(3, N);
for i = 1:3*N
    if rem(i - 1, 3) == 0
        D(1, floor(i/3) + 1) = params(i);
    elseif rem(i - 2, 3) == 0
        D(2, floor(i/3) + 1) = params(i);
    else
        D(3, floor(i/3)) = params(i);
    end
end
D = D./(ones(3,1)*sum(D.^2).^(1/2));

Dnx = ones(N,1)*D(1,:);
Dny = ones(N,1)*D(2,:);
Dnz = ones(N,1)*D(3,:);

Dmx = D(1,:)'*ones(1,N);
Dmy = D(2,:)'*ones(1,N);
Dmz = D(3,:)'*ones(1,N);

% rhatx, rhaty, and rhatz are matrices that relate the positions of each particle and are initiatialized elsewhere. This has been omitted for brevity
DmR = rhatx.*Dmx+rhaty.*Dmy+rhatz.*Dmz;

% Magnetic Field calculations (M is a scalar that represents the particle's magnetization, V a constant that represents particle volume, r is a matrix that represents the distance between each pair of particles. Also initiated elsewhere)
Hddx = (M*V/4/pi)*sum((1./r.^3-eye(N)).*(3*DmR.*rhatx-Dmx));
Hddy = (M*V/4/pi)*sum((1./r.^3-eye(N)).*(3*DmR.*rhaty-Dmy));
Hddz = (M*V/4/pi)*sum((1./r.^3-eye(N)).*(3*DmR.*rhatz-Dmz));

Hx = H0x*ones(1,N) + Hddx;
Hy = H0y*ones(1,N) + Hddy;
Hz = H0z*ones(1,N) + Hddz;

% Calculate the dot product of H and D
DnH = Dnx(1,:).*Hx + Dny(2,:).*Hy+ Dnz(2,:).*Hz;

Wx = (QQQ)*(Hx-Dnx(1,:).*DnH);
Wy = (QQQ)*(Hy-Dny(1,:).*DnH);
Wz = (QQQ)*(Hz-Dnz(1,:).*DnH);

T = [Wx;Wy;Wz];

dydt = T(1:3*N)';
end

Answer 1

[James Tursa]

I haven't looked through your code, but from your description it sounds like you are making a fundamental error in your state space. You have N particles in a magnetic field, so the state of each of these particles will involve both position and velocity in 3-space, and be governed by 2nd order differential equations according to the inverse square law and the local magnetic field. So your state space should be 6xN in size, where the 6 elements for each particle are 3D position and 3D velocity, and the derivative for each particle will be 3D velocity and 3D acceleration. You need to rewrite your entire code to account for this.

[James Zola]

Hi, thanks for taking the time to answer my question. I'm less concerned with the physics at this point and more with MATLAB's options for numeric integration. I have a more complete version of the code that maybe accounts for some of the concerns you raised here. I can include it if you like, this version makes some assumptions in order to simplify the presentation. The actual code is much more lengthy and accounts for angular velocity as you suggest. This version simply assumes angular acceleration to be zero and is governed by the below equation:

The expression is summed up by the constant QQQ in my code.

[James Tursa]

OK, so if you put my concerns aside and assume there is only a 1st order equation to integrate as you show, then I see no reason why a 3xN system wouldn't work with ode45() as long as you have the derivative function correctly coded. That part I can't comment on because you don't show the formula for Hn.

[James Zola]

Ok, thanks for your help!