Commit 61da078a authored by Håkon Harnes's avatar Håkon Harnes

completed CP10, 11, 1, 3

parent b45116f6
% 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
% task a
figure(
orbit([0 100], [2, 0.2, 2, -0.2, 0, 0, 0, 0, -2, -0.2, -2, 0.2], 0.01, 10)
% task b
orbit([0 100], [2, 0.2, 2, -0.2, 0, 0, 0, 0, -2, -0.2, -2, 0.2], 0.01, 10)
\ 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
function orbit(int, ic, n, p)
% plots n points
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);
% y-vector
y(1, :) = [x1 vx1 y1 vy1 x2 vx2 y2 vy2 x3 vx3 y3 vy3];
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');
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, ...
'xdata', [], 'ydata', []);
% tails
l1 = animatedline('Color', 'b', 'LineWidth', 3);
l2 = animatedline('Color', 'b', 'LineWidth', 3);
l3 = 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);
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 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;
end
end
\ No newline at end of file
% Program 6.4 Plotting program for one-body problem
% Inputs: time interval inter, initial conditions
% ic = [x0 vx0 y0 vy0], x position, x velocity, y pos, y vel,
% number of steps n, steps per point plotted p
% Calls a one-step method such as trapstep.m
% Example usage: orbit([0 100], [0 1 2 0], 10000, 5)
function z = orbit(int, ic, h, p)
% plots n points
n = round((int(2) - int(1)) / (p*h));
% 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);
% y-vector
y(1, :) = [x1 vx1 y1 vy1 x2 vx2 y2 vy2 x3 vx3 y3 vy3];
t(1) = int(1);
% plotting
set(gca, 'XLim', [-50 50], 'YLim', [-50 50], 'XTick', [-50 0 50], 'YTick', ...
[-5 0 5], 'Drawmode', 'fast', 'Visible', 'on');
cla;
%{
% dot in middle
line('color', 'r', 'Marker', '.', 'markersize', 25, ...
'xdata', 0, 'ydata', 0);
drawnow;
%}
% heads
head1=line('color', 'r', 'Marker', '.', 'markersize', 25, ...
'erase', 'xor', 'xdata', [], 'ydata', []);
head2=line('color', 'r', 'Marker', '.', 'markersize', 25, ...
'erase', 'xor', 'xdata', [], 'ydata', []);
head3=line('color', 'r', 'Marker', '.', 'markersize', 25, ...
'erase', 'xor', 'xdata', [], 'ydata', []);
% tails
l1 = animatedline('Color', 'b', 'LineWidth', 2);
l2 = animatedline('Color', 'b', 'LineWidth', 2);
l3 = animatedline('Color', 'b', 'LineWidth', 2);
for k = 1:n/p
for i = 1:p
t(i+1) = t(i) + h;
y(i+1, :) = eulerstep(t(i), y(i, :), h);
end
y(1, :) = y(p+1, :);
t(1) = t(p+1);
set(head,'xdata', y(1, 1), 'ydata', y(1, 3))
addpoints(l, y(2:p, 1), y(2:p, 3));
drawnow;
end
%}
end
\ No newline at end of file
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);
end
\ No newline at end of file
function z = ydot(t, x)
g = 1;
m1 = 0.3;
m2 = 0.03;
m3 = 0.03;
mg1 = m1*g;
mg2 = m2*g;
mg3 = m3*g;
mg2 = m2*g;
px2 = 0;
py2 = 0;
px1 = x(1); py1 = x(3); vx1 = x(2); vy1 = x(4);
px2 = x(5); py1 = x(7); vx1 = x(6); vy1 = x(8);
px3 = x(9); py1 = x(3); vx1 = x(2); vy1 = x(4);
dist = sqrt((px2-px1)^2+(py2-py1)^2);
z = zeros(1,4);
z(1) = vx1;
z(2) = (mg2*(px2-px1))/(dist^3);
z(3) = vy1;
z(4) = (mg2*(py2-py1))/(dist^3);
end
\ No newline at end of file
% task a
figure('Name', 'Task a')
orbit([0 100], [-0.970, -0.466, 0.243, -0.433, 0.970, -0.466, -0.243, -0.433, 0, 0.932, 0, 0.866], 100000, 5)
% task b
for k = 1:5
figure('Name', 'Task b')
orbit([0 100], [-0.970, -0.466, 0.243, -0.433, 0.970, -0.466, -0.243, -0.433, 0, 0.932 + (10^-k), 0, 0.866], 100000, 5)
end
\ 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
function orbit(int, ic, n, p)
% plots n points
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);
% y-vector
y(1, :) = [x1 vx1 y1 vy1 x2 vx2 y2 vy2 x3 vx3 y3 vy3];
t(1) = int(1);
% sets up plot
set(gca, 'XLim', [-5 5], 'YLim', [-5 5], 'XTick', [-5 0 5], 'YTick', ...
[-5 0 5], '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, ...
'xdata', [], 'ydata', []);
% tails
l1 = animatedline('Color', 'b', 'LineWidth', 3);
l2 = animatedline('Color', 'b', 'LineWidth', 3);
l3 = 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);
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 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;
end
end
\ No newline at end of file
function z = ydot(~, x)
g = 1;
m1 = 1;
m2 = 1;
m3 = 1;
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);
end
\ No newline at end of file
[t, y, e] = midpointmethod([0 1], 1, 10);
table(t', y', e', 'VariableNames', {'t', 'y', 'e'})
% exact solution of ydot
function y = exact(t)
y = exp(t^3 / 3);
end
% performs the trapezoid method on function ydot
function [t, y, e] = midpointmethod(inter, y0, n)
t(1) = inter(1);
y(1) = y0;
e(1) = 0;
h = (inter(2) - inter(1)) / n;
for i = 1:n
t(i+1) = t(i) + h;
y(i+1) = midpointstep(t(i), y(i), h);
e(i+1) = exact(t(i+1)) - y(i+1);
end
end
function y = midpointstep(t, y, h)
y = y + h * ydot(t + (h/2), y + (h/2) * ydot(t, y));
end
% right side of the ODE
function z = ydot(t, y)
z = t^2 * y;
end
% 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
function y = eulerstep(t, x, h)
y = x + h * ydot(t, x);
end
\ No newline at end of file
function orbit(int, ic, n, p)
% plots n points
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);
% y-vector
y(1, :) = [x1 vx1 y1 vy1 x2 vx2 y2 vy2 x3 vx3 y3 vy3];
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');
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, ...
'xdata', [], 'ydata', []);
% tails
l1 = animatedline('Color', 'b', 'LineWidth', 3);
l2 = animatedline('Color', 'b', 'LineWidth', 3);
l3 = 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);
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 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;
end
end
\ No newline at end of file
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);
end
\ No newline at end of file
% h = 0.1
figure('Name', 'h = 0.1');
[t, w, y, e] = rkmethod([0 1], 1, 10);
table(t', w', y', e', 'VariableNames', {'t', 'w', 'y', 'e'})
% h = 0.05
figure('Name', 'h = 0.05');
[t, w, y, e] = rkmethod([0 1], 1, 20);
table(t', w', y', e', 'VariableNames', {'t', 'w', 'y', 'e'})
% h = 0.025
figure('Name', 'h = 0.025');
[t, w, y, e] = rkmethod([0 1], 1, 40);
table(t', w', y', e', 'VariableNames', {'t', 'w', 'y', 'e'})
% exact solution of ydot
function y = exact(t)
y = exp(t^3 / 3);
end
% performs the trapezoid method on function ydot
function [t, w, y, e] = rkmethod(inter, y0, n)
t(1) = inter(1);
w(1) = y0;
y(1) = y0;
e(1) = 0;
h = (inter(2) - inter(1)) / n;
for i = 1:n t(i+1) = t(i) + h;
w(i+1) = rkstep(t(i), w(i), h);
y(i+1) = exact(t(i+1));
e(i+1) = y(i+1) - w(i+1);
end
plot(t, w, 'LineWidth', 2);
end
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
% right side of the ODE
function z = ydot(t, y)
z = t^2 * y;
end
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