Ode45 producing NaN values only

Hello,

I am working on a matlab program to plot the x,y positions of 5 particles based on a differential equation describing their motion and the influence that each particle has on eachother's motions. I attempted to make the code scalable to any N number of particles. To run the program you call the following:

clc
clear

[t,s] = ode45('NthSwarmDynamics', [0 50], [1 1 0, 2 2 0, 3 3 0, 4 4 0, 5 5 0, 1 1 0, 1 1 0, 1 1 0, 1 1 0, 1 1 0]);

%plot X1_x and X1_y
plot3(s(:,1),s(:,2),t,':b');
hold on;

%plot X2_x and X2_y
plot3(s(:,4),s(:,5),t,':r');
hold on;

%plot X3_x and X3_y
plot3(s(:,7),s(:,8),t,':g');
hold on;

%plot X4_x and X4_y
plot3(s(:,10),s(:,11),t,':m');
hold on;

%plot X5_x and X5_y
plot3(s(:,13),s(:,14),t,':c');
hold on;

The input to this code is t, and s. s is further broken down such that the first half of the vector is the x,y,z positions for each particle (3 coordinates X 5 particles = 15 values, seen here at 1,1,0 2,2,0, etc) and the second half of the vector is the x,y,z velocities for each particle.

The NthSwarmDynamics file follows:

%SwarmDynamics: Models the dynamics of a N-agent swarm in 3D
%S is the position variable
%V is the velocuty variable

function dSdt = NthSwarmDynamics(t, Z)

dimZ = length(Z);
S = Z(1:dimZ/2);
V = Z(dimZ/2+1:dimZ);

N = length(S)/3;
U_x = zeros(N,N);
U_y = zeros(N,N);
U_z = zeros(N,N);
alpha = 1;
beta = 1;
mass = 1;
C_a = 1;
l_a = 2;
C_r = 2;
l_r = 1;

for i = 1:N
    for j = 1:N
        if i == j
            U_x(i,j) = 0;
            U_y(i,j) = 0;
            U_z(i,j) = 0;
        else
            r_x = S(i) -S(3*(j-1)+1);
            r_y = S(i+1) -S(3*(j-1)+2);
            r_z = S(i+2) -S(3*(j-1)+3);
            M = sqrt((S(i)-S(3*(j-1)+1))^2+(S(i+1)-S(3*(j-1)+2))^2+(S(i+2)-S(3*(j-1)+3))^2);

              U_x(i,j) = (C_a/l_a)*exp(-M/l_a)*(r_x/M)-(C_r/l_r)*exp(-M/l_r)*(r_x/M);
              U_y(i,j) = (C_a/l_a)*exp(-M/l_a)*(r_y/M)-(C_r/l_r)*exp(-M/l_r)*(r_y/M);
              U_z(i,j) = (C_a/l_a)*exp(-M/l_a)*(r_z/M)-(C_r/l_r)*exp(-M/l_r)*(r_z/M);
          end
      end
  end

dSdt = zeros(length(Z),1);

j=1;
u=1;

for i=1:6:length(dSdt)
    dSdt(i) = V(j);
    dSdt(i+1) = V(j+1);
    dSdt(i+2) = V(j+2);
    dSdt(i+3) = ((alpha - beta * ((sqrt(V(j)^2 + V(j+1)^2 + V(j+2)^2))^2)*V(j))- sum(U_x(u,1:N)))/mass;
    dSdt(i+4) = ((alpha - beta * ((sqrt(V(j)^2 + V(j+1)^2 + V(j+2)^2))^2)*V(j+1))- sum(U_y(u,1:N)))/mass;
    dSdt(i+5) = ((alpha - beta * ((sqrt(V(j)^2 + V(j+1)^2 + V(j+2)^2))^2)*V(j+2))- sum(U_z(u,1:N)))/mass;
    j = j+3;
    u = u+1;
end

When I run the swarm plotter function, I got an empty plot and upon inspecting my s vector, there is nothing but NaN values. Does anyone know where my flaw is?

Thank you so much for your time, Jeremy

Answer 1

[James Tursa]

I haven't looked over your code in detail, but this definitely looks funny:

for i = 1:N
    for j = 1:N
          :
          r_x = S(i) -S(3*(j-1)+1);
          r_y = S(i+1) -S(3*(j-1)+2);
          r_z = S(i+2) -S(3*(j-1)+3);
          M = sqrt((S(i)-S(3*(j-1)+1))^2+(S(i+1)-S(3*(j-1)+2))^2+(S(i+2)-S(3*(j-1)+3))^2);

If N is the number of particles, why is your indexing for the i'th particle i,i+1,i+2? Shouldn't it be the same pattern as the j'th particle using the factor of 3 since you have x,y,z positions for both? E.g., I would have expected something like this:

for i = 1:N
    I = 3*(i-1);
    for j = 1:N
          J = 3*(j-1);
          :
          r_x = S(I+1) - S(J+1);
          r_y = S(I+2) - S(J+2);
          r_z = S(I+3) - S(J+3);
          M = sqrt((S(I+1)-S(J+1))^2+(S(I+2)-S(J+2))^2+(S(I+3)-S(J+3))^2);

Of course, some of this could be simplified using vector operations. E.g.,

          M = norm(S(I+1:I+3)-S(J+1:J+3));

[Jeremy Friedman]

That's a good spot, James, thanks for that!

I've gone through some hand calculations and your suggestion does seem to produce the correct logic, though the matlab result is still a vector of NaN values.

If I run:

NthSwarmDynamics([0 100], [1 1 0 2 2 0 3 3 0 1 1 0 1 1 0 1 1 0])

I get the result:

ans =
      1.0000
      1.0000
           0
         NaN
         NaN
         NaN
      1.0000
      1.0000
           0
         NaN
         NaN
         NaN
      1.0000
      1.0000
           0
     -0.1322
     -0.1322
      0.2000

Which leads me to believe that the problem lies in the following section:

U_x(i,j) = (C_a/l_a)*exp(-M/l_a)*(r_x/M)-(C_r/l_r)*exp(-M/l_r)*(r_x/M);
U_y(i,j) = (C_a/l_a)*exp(-M/l_a)*(r_y/M)-(C_r/l_r)*exp(-M/l_r)*(r_y/M);
U_z(i,j) = (C_a/l_a)*exp(-M/l_a)*(r_z/M)-(C_r/l_r)*exp(-M/l_r)*(r_z/M);

But the code is not displaying the U_x, U_y, U_z matrices in the Workspace for me to verifty. Does anyone know how I can access the contents of these matrices?

[Jeremy Friedman]

After some further testing using

NthSwarmDynamics([0 100], [1 1 0 2 2 0 3 3 0 1 1 0 1 1 0 1 1 0])

I have found that removing the sum() term from

dSdt(i+3) = ((alpha - beta * ((sqrt(V(j)^2 + V(j+1)^2 + V(j+2)^2))^2)*V(j))- sum(U_x(u,1:N)))/mass;
dSdt(i+4) = ((alpha - beta * ((sqrt(V(j)^2 + V(j+1)^2 + V(j+2)^2))^2)*V(j+1))- sum(U_y(u,1:N)))/mass;
dSdt(i+5) = ((alpha - beta * ((sqrt(V(j)^2 + V(j+1)^2 + V(j+2)^2))^2)*V(j+2))- sum(U_z(u,1:N)))/mass;

generates values instead of NaN.

However, when running this using the separate file which calls ode45, I still get nothing but NaN.

[James Tursa]

"But the code is not displaying the U_x, U_y, U_z matrices in the Workspace for me to verifty. Does anyone know how I can access the contents of these matrices?"

The easiest way is to set a breakpoint in your function. When you run the code, it will pause at the breakpoint will all current workspace variables available for inspection/modification. To set a breakpoint, in the editor simply click on the dash that is just to the right of the line number and it will turn into a red dot. The code will pause at this line just prior to executing the line and pop you out to the command window. At the command window you will have direct access to the current function workspace variables. You can then use the Step, Continue, Etc buttons to manually direct the subsequent code execution. You can also set other breakpoints and clear breakpoints when you are paused.

[James Tursa]

M will be identically 0 if any of your two particles are at the same position. Then when you divide by M you risk getting either an inf or a NaN result. In your example input vector:

[1 1 0 2 2 0 3 3 0 1 1 0 1 1 0 1 1 0]

It appears to me that the 1st particle is at [1 1 0] and the last particle is also at [1 1 0]. This will lead to your observed problems. Your equations of motion are not valid if two particles are at the same position, so you can't have initial conditions that have this.

[Jeremy Friedman]

The last three sets of three values represent the initial velocity though, not the initial position. This is due to

dimZ = length(Z);
S = Z(1:dimZ/2);
V = Z(dimZ/2+1:dimZ);

So my assumption is that after ode45 runs, the resulting matrix is broken up into: the first half, which can be plotted as positions, and the second half, which can be plotted as velocities.

[James Tursa]

Ah yes, I see that now. But then I am confused by the following which led me to believe the r's and v's were interleaved:

    dSdt(i) = V(j);
    dSdt(i+1) = V(j+1);
    dSdt(i+2) = V(j+2);
    dSdt(i+3) = ((alpha - beta * ((sqrt(V(j)^2 + V(j+1)^2 + V(j+2)^2))^2)*V(j))- sum(U_x(u,1:N)))/mass;
    dSdt(i+4) = ((alpha - beta * ((sqrt(V(j)^2 + V(j+1)^2 + V(j+2)^2))^2)*V(j+1))- sum(U_y(u,1:N)))/mass;
    dSdt(i+5) = ((alpha - beta * ((sqrt(V(j)^2 + V(j+1)^2 + V(j+2)^2))^2)*V(j+2))- sum(U_z(u,1:N)))/mass;

The above code appears to interleave the V's and the accel's rather than having all of the V's in the first half of the return vector and all of the accel's in the last half of the return vector (which would seem to be required for your ordering description). I think you need to reorder these assignments.

[Jeremy Friedman]

In that case, the return vector would be organized such that every 6 values represent 1 particle, with the first 3 being position and the second 3 being velocity.

So to plot the position change with time I was thinking I could use something like:

%plot X1_x and X1_y
plot3(s(:,1),s(:,2),t,':b');
hold on;
%plot X2_x and X2_y
plot3(s(:,7),s(:,8),t,':r');
hold on;
%plot X3_x and X3_y
plot3(s(:,13),s(:,14),t,':g');
hold on;

which would just plot the x and y coordinates for three particles.

The good news is after implementing your indexing fix the result is no longer NaN!

However, the result now is nonsensical involving tight clumping of most particles and one particle shooting off to infinity, so it looks like I have a new problem. But thank you very much for helping me with the NaN issue!

[James Tursa]

Well, you just need to make a decision about how the states are organized and stick with it. Your posted code is a mixed bag. Either have the particles interleaved pos/vel/pos/vel/etc, or have all the pos followed by all the vel.