Trajectorys of 4 objects on 3D plane not symmetrical

Hi all,

Sorry i'm a beginner at Matlab and am having issues with getting 4 spots on a 3D plane with the intitial conditions (0,0,0), (1,0,3), (1,2,0), (0,2,3) to symmetrically meet at (0.5,1,1.5).

My 2D simulation appeared like this:

But from my 3D simulation they appeared like this:

Im unsure what line of code I need to change to fix this from my code.

My 3D code is given below:

clc, clear, clear all

ti=0;
tf=1000;
tspan= linspace(ti,tf,1000000);
f0 = [0; 0; 0; 1; 2; 0; 1; 0; 3; 0; 2; 3];
[t, f] = ode45(@g, tspan, f0);

%Trajectories are done in 3 coordinates x,y,z
i_1=(f(:,1)); %Trajectory of bug 1
j_1=(f(:,5));
k_1=(f(:,9));

i_2=(f(:,2)); %Trajectory of bug 2
j_2=(f(:,6));
k_2=(f(:,10));

i_3=(f(:,3)); %Trajectory of bug 3
j_3=(f(:,7));
k_3=(f(:,11));

i_4=(f(:,4)); %Trajectory of bug 4
j_4=(f(:,8));
k_4=(f(:,12));

figure
plot3(i_1, j_1, k_1)
hold on
plot3(i_2, j_2, k_2)
plot3(i_3, j_3, k_3)
plot3(i_4, j_4, k_4)
hold off
title('test')
xlabel('x')
ylabel('y')
zlabel('z')

function dxdt = g(t,f)
a = sqrt((f(2)-f(1))^2 + (f(6)-f(5))^2 + (f(10)-f(9))^2);
b = sqrt((f(3)-f(2))^2 + (f(7)-f(6))^2 + (f(11)-f(10))^2);
c = sqrt((f(4)-f(3))^2 + (f(8)-f(7))^2 + (f(12)-f(11))^2);
d = sqrt((f(1)-f(4))^2 + (f(5)-f(8))^2 + (f(9) -f(12))^2);

dxdt = [ ...
    (f(2)  - f(1))  / a;
    (f(6)  - f(5))  / a;
    (f(10) - f(9))  / a;
    (f(3)  - f(2))  / b;
    (f(7)  - f(6))  / b;
    (f(11) - f(10)) / b;
    (f(4)  - f(3))  / c;
    (f(8)  - f(7))  / c;
    (f(12) - f(11)) / c;
    (f(1)  - f(4))  / d;
    (f(5)  - f(8))  / d;
    (f(9)  - f(12)) / d];
end

%My 2D code is here:

tinitial = 0;
tfinal = 1000;
tspan = linspace(tinitial, tfinal, 1000000);
f0 = [0; 0; 1; 1; 0; 2; 2; 0];
[t, f] = ode45(@g, tspan, f0);
j = (f(:, 1));
k = (f(:, 5));
a = (f(:, 2));
b = (f(:, 6));
d = (f(:, 3));
e = (f(:, 7));
h = (f(:, 4));
i = (f(:, 8));
figure
hold on
plot3(j, k, tspan)
plot3(a, b, tspan)
plot3(d, e, tspan)
plot3(h, i, tspan)
title('test')
xlabel('x')
ylabel('y')
zlabel('t')
function dfdt = g(t,f)
dfdt = [(f(2) - f(1))/(sqrt((f(2)-f(1))^2 + (f(6)-f(5))^2));
    (f(3)-f(2))/(sqrt((f(3)-f(2))^2 + (f(7)-f(6))^2));
    (f(4) - f(3))/(sqrt((f(4)-f(3))^2 + (f(8)-f(7))^2));
    (f(1) - f(4))/(sqrt((f(1)-f(4))^2 + (f(5)-f(8))^2));
    (f(6) - f(5))/(sqrt((f(2)-f(1))^2 + (f(6)-f(5))^2));
    (f(7) - f(6))/(sqrt((f(3)-f(2))^2 + (f(7)-f(6))^2));
    (f(8) - f(7))/(sqrt((f(4)-f(3))^2 + (f(8)-f(7))^2));
    (f(5) - f(8))/(sqrt((f(1)-f(4))^2 + (f(5)-f(8))^2))];
end

Cheers

[Torsten]

I you change the plot command in the 2d-case as you did in the 3d-case, are the curves still as expected ?

figure
plot3(j,k,tspan)
hold on
plot3(a,b,tspan)
plot3(d,e,tspan)
plot3(h,i,tspan)
hold off
title('test')
xlabel('x')
ylabel('y')
zlabel('t')

[Jonas Freiheit]

Yes

Answer 1

[Bjorn Gustavsson]

Your main problem (or rather: step to take) is that you do the calculations element-by-element in the different arrays. That is not the way to go through life. You will have a much easier time if you start working with the arrays using vector-algebra. If you do that correctly there is no real difference between the 2-D and 3-D case. For your problem I guess that at each time-step you want to move the points towards the next point, some fraction of the distance, or something like that. For this you should use something like (for moving point 2 towards point 1):

dr21 = r1 - r2;
l_21 = norm(dr21);
lstep = l_21/30;
r2_next = r2 + (lstep/l_21)*dr21;

Then you have to put this into your loop and do it for all 4 points.

This way of working with vector-algebra removes (almost) all the thinking you have to do for a problem like this. 5-10 minutes with pen and paper sketching arrows and scribling trivial additions and subtractions will cut down the time required by factors of 10. Not just a single factor of 10, factors...

HTH

[Jonas Freiheit]

Thats a genius idea to work through that problem. I was wondering what r1 and r2 are defined as?

[Bjorn Gustavsson]

Well Those should be the first elements of the bug-trajectories you extract so I guess something like this:

r1 = f0([1 5 9]); 
r2 = f0([1 5 9]+1);
r3 = f0([1 5 9]+2);
r4 = f0([1 5 9]+3);

(I'm not sure this qualifies as genious, but I'll take it - on the othe hand "vector algebra" certainly does qualify!)

[Jonas Freiheit]

is r2 = f0([1 5 9]+1) the same as f0([2 6 10]).

Also there is no loop in my program, is that talking about ode45?

[Jonas Freiheit]

Also what is the 30?

[Bjorn Gustavsson]

Ah, sorry about that. I just realised that I confused/combined your question with someone elses where the task was to compare the ODE-solution with a simple fixed-size stepping solution. For that "someonelses problem" I thought that a step-length of 1/30 might be sensible if starting at coordinates of single-digit magnitudes. Just cut out that constant - in your setting I guess it would simply be something like a scalling of time.

[Bjorn Gustavsson]

Yes

r2 = f0([1 5 9]+1);
r2 == f0([1 5 9]+1);

Matlab's indexing capabilities are ridiculously handy! It lets you do things like:

A = [1 2;3 4];
A = zeros(size(A));
A([1 3 4]) = rand(1,3);
C = zeros(size(A));   
C(:) = randn(4,1);

which can be increadibly useful - not that these examples are, but these capacities have been very handy for problems I've had to solve.

Glad that the suggestion helped. If you have time once you've solved the problem we'd be happy to see the result - everyone likes a nice figure.