Trajectorys of 4 objects on 3D plane not symmetrical

조회 수: 1 (최근 30일)
Jonas Freiheit
Jonas Freiheit 2022년 5월 16일
댓글: Bjorn Gustavsson 2022년 5월 18일
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
  댓글 수: 2
Torsten
Torsten 2022년 5월 16일
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')

댓글을 달려면 로그인하십시오.

채택된 답변

Bjorn Gustavsson
Bjorn Gustavsson 2022년 5월 16일
편집: Bjorn Gustavsson 2022년 5월 16일
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
  댓글 수: 6
Bjorn Gustavsson
Bjorn Gustavsson 2022년 5월 18일
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
Bjorn Gustavsson 2022년 5월 18일
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.

댓글을 달려면 로그인하십시오.

추가 답변 (0개)

카테고리

Help CenterFile Exchange에서 Pulsed Waveforms에 대해 자세히 알아보기

제품

Community Treasure Hunt

Find the treasures in MATLAB Central and discover how the community can help you!

Start Hunting!

Translated by