Commit c0284a5d authored by Håkon Harnes's avatar Håkon Harnes

completed assignment 5

parent 61da078a
% task a
figure(1)
orbit([0 100], [2, 0.2, 2, -0.2, 0, 0, 0, 0, -2, -0.2, -2, 0.2], 10000, 5)
% task b
figure(2)
orbit([0 100], [2, 0.20001, 2, -0.2, 0, 0, 0, 0, -2, -0.2, -2, 0.2], 10000, 5)
\ No newline at end of file
orbit([0 100], [5, 5, 5], 100000, 5)
\ No newline at end of file
function y = eulerstep(t, x, h)
y = x + h * ydot(t, x);
end
\ No newline at end of file
......@@ -3,51 +3,42 @@ function orbit(int, ic, n, p)
h = (int(2) - int(1)) / n;
% initial conditons
x1 = ic(1); vx1 = ic(2); y1 = ic(3); vy1 = ic(4);
x2 = ic(5); vx2 = ic(6); y2 = ic(7); vy2 = ic(8);
x3 = ic(9); vx3 = ic(10); y3 = ic(11); vy3 = ic(12);
x0 = ic(1); vx0 = ic(2); y0 = ic(3);
% y-vector
y(1, :) = [x1 vx1 y1 vy1 x2 vx2 y2 vy2 x3 vx3 y3 vy3];
y(1, :) = [x0 vx0 y0];
t(1) = int(1);
% sets up plot
set(gca, 'XLim', [-25 25], 'YLim', [-25 25], 'XTick', [-25 0 25], 'YTick', ...
[-25 0 25], 'Visible', 'on');
set(gca, 'XLim', [-25 25], 'YLim', [0 50], 'XTick', [-25 0 25], 'YTick', ...
[0 25 50], 'Visible', 'on');
cla;
% heads
head1=line('color', 'r', 'Marker', '.', 'markersize', 25, ...
'xdata', [], 'ydata', []);
head2=line('color', 'r', 'Marker', '.', 'markersize', 25, ...
'xdata', [], 'ydata', []);
head3=line('color', 'r', 'Marker', '.', 'markersize', 25, ...
head = line('color', 'r', 'Marker', '.', 'markersize', 25, ...
'xdata', [], 'ydata', []);
% tails
l1 = animatedline('Color', 'b', 'LineWidth', 3);
l2 = animatedline('Color', 'b', 'LineWidth', 3);
l3 = animatedline('Color', 'b', 'LineWidth', 3);
l = animatedline('Color', 'b', 'LineWidth', 3);
% calculates data
for k = 1:n/p
for i = 1:p
t(i+1) = t(i) + h;
y(i+1, :) = eulerstep(t(i), y(i, :), h);
y(i+1, :) = rkstep(t(i), y(i, :), h);
end
y(1, :) = y(p+1, :);
t(1) = t(p+1);
% sets heads
set(head1,'xdata', y(1, 1), 'ydata', y(1, 3))
set(head2,'xdata', y(1, 5), 'ydata', y(1, 7))
set(head3,'xdata', y(1, 9), 'ydata', y(1, 11))
% sets head
set(head,'xdata', y(1, 1), 'ydata', y(1, 3))
% sets tail
addpoints(l, y(2:p, 1), y(2:p, 3));
% sets tails
addpoints(l1, y(2:p, 1), y(2:p, 3));
addpoints(l2, y(2:p, 5), y(2:p, 7));
addpoints(l3, y(2:p, 9), y(2:p, 11));
% draws data
drawnow;
......
function y = rkstep(t, w, h)
s1 = ydot(t, w);
s2 = ydot(t + (h/2), w + (h/2)*s1);
s3 = ydot(t + (h/2), w + (h/2)*s2);
s4 = ydot(t + h, w + h*s3);
y = w + (h/6) * (s1 + 2*s2 + 2*s3 + s4);
end
function z = ydot(~, x)
g = 1;
m1 = 0.3;
m2 = 0.03;
m3 = 0.03;
mg1 = m1*g;
mg2 = m2*g;
mg3 = m3*g;
px1 = x(1); py1 = x(3); vx1 = x(2); vy1 = x(4);
px2 = x(5); py2 = x(7); vx2 = x(6); vy2 = x(8);
px3 = x(9); py3 = x(11); vx3 = x(10); vy3 = x(12);
dist1 = sqrt((px2-px1)^2 + (py2-py1)^2);
dist2 = sqrt((px2-px3)^2 + (py2-py3)^2);
dist3 = sqrt((px3-px1)^2 + (py3-py1)^2);
z = zeros(1, 8);
z(1) = vx1;
z(2) = (mg2*(px2-px1))/(dist1^3)+(mg3*(px3-px1))/(dist3^3);
z(3) = vy1;
z(4) = (mg2*(py2-py1))/(dist1^3)+(mg3*(py3-py1))/(dist3^3);
z(5) = vx2;
z(6) = (mg1*(px1-px2))/(dist1^3)+(mg3*(px3-px2))/(dist2^3);
z(7) = vy2;
z(8) = (mg1*(py1-py2))/(dist1^3)+(mg3*(py3-py2))/(dist2^3);
z(9) = vx3;
z(10) = (mg1*(px1-px3))/(dist3^3)+(mg2*(px2-px3))/(dist2^3);
z(11) = vy3;
z(12) = (mg1*(py1-py3))/(dist3^3)+(mg2*(py2-py3))/(dist2^3);
function z = ydot(~, y)
s = 10;
r = 28;
b = 8/3;
z(1) = -s*y(1) + s*y(2);
z(2) = -y(1)*y(3) + r*y(1) - y(2);
z(3) = y(1)*y(2) - b*y(3);
end
\ No newline at end of file
Markdown is supported
0%
or
You are about to add 0 people to the discussion. Proceed with caution.
Finish editing this message first!
Please register or to comment