- #1
aoeuDvorakhtn
- 2
- 0
1. For this problem, we are asked to design a MatLab/Octave code that will calculate and plot the orbit of a sattelite
2. We are given the following 2nd order equations
[itex] x''= \frac{-x}{( \sqrt{x^{2}+y^{2}})^{3}}[/itex]
[itex]y''= \frac{-y}{( \sqrt{x^{2}+y^{2}})^{3}}[/itex]
We can convert this into the following system of 4 1st order DEs by substitiuting z=x' and w=y'
[itex]x'=z[/itex]
[itex]y'=w[/itex]
[itex]z'= \frac{-x}{( \sqrt{x^{2}+y^{2}})^{3}}[/itex]
[itex]w'= \frac{-y}{( \sqrt{x^{2}+y^{2}})^{3}}[/itex]
We use the following initial vallues
[x0 y0 z0 w0] = [4 0 0 0.5]
I managed to build a program that works well enough that, after adding annotations, I could hand it in. However, as you can see below, the lines within the for loop could use a little help as far as elegance goes (I've attached the .m file if my copy/paste is too difficult to read). Ideally, I'd like to be able to make those lines a little more compact. Any ideas?
X0=[4 0 0 0.5];
t0=0;
tmax=50;
h=0.25;
SOL=[t0 X0;];
f1=inline('z','z');
f2=inline('w','w');
f3=inline('-x/((sqrt(x^2+y^2))^3)','x','y');
f4=inline('-y/((sqrt(x^2+y^2))^3)','x','y');
for i = 1:tmax/h
g1=[f1(X0(3)) f2(X0(4)) f3(X0(1),X0(2)) f4(X0(1),X0(2))];
g2=[f1(X0(3)+0.5*h*g1(3)) f2(X0(4)+0.5*h*g1(4)) f3(X0(1)+0.5*h*g1(1),X0(2)+0.5*h*g1(2)) f4(X0(1)+0.5*h*g1(1),X0(2)+0.5*h*g1(2))];
g3=[f1(X0(3)+0.5*h*g2(3)) f2(X0(4)+0.5*h*g2(4)) f3(X0(1)+0.5*h*g2(1),X0(2)+0.5*h*g2(2)) f4(X0(1)+0.5*h*g2(1),X0(2)+0.5*h*g2(2))];
g4=[f1(X0(3)+h*g3(3)) f2(X0(4)+h*g3(4)) f3(X0(1)+h*g3(1),X0(2)+h*g3(2)) f4(X0(1)+h*g3(1),X0(2)+h*g3(2))];
X1=X0+(h/6)*(g1+2*g2+2*g3+g4);
t1=t0+h;
SOL=[SOL;t1 X1];
X0=X1;
t0=t1;
end
plot(SOL(:,2),SOL(:,3))
axis('equal')
2. We are given the following 2nd order equations
[itex] x''= \frac{-x}{( \sqrt{x^{2}+y^{2}})^{3}}[/itex]
[itex]y''= \frac{-y}{( \sqrt{x^{2}+y^{2}})^{3}}[/itex]
We can convert this into the following system of 4 1st order DEs by substitiuting z=x' and w=y'
[itex]x'=z[/itex]
[itex]y'=w[/itex]
[itex]z'= \frac{-x}{( \sqrt{x^{2}+y^{2}})^{3}}[/itex]
[itex]w'= \frac{-y}{( \sqrt{x^{2}+y^{2}})^{3}}[/itex]
We use the following initial vallues
[x0 y0 z0 w0] = [4 0 0 0.5]
The Attempt at a Solution
I managed to build a program that works well enough that, after adding annotations, I could hand it in. However, as you can see below, the lines within the for loop could use a little help as far as elegance goes (I've attached the .m file if my copy/paste is too difficult to read). Ideally, I'd like to be able to make those lines a little more compact. Any ideas?
X0=[4 0 0 0.5];
t0=0;
tmax=50;
h=0.25;
SOL=[t0 X0;];
f1=inline('z','z');
f2=inline('w','w');
f3=inline('-x/((sqrt(x^2+y^2))^3)','x','y');
f4=inline('-y/((sqrt(x^2+y^2))^3)','x','y');
for i = 1:tmax/h
g1=[f1(X0(3)) f2(X0(4)) f3(X0(1),X0(2)) f4(X0(1),X0(2))];
g2=[f1(X0(3)+0.5*h*g1(3)) f2(X0(4)+0.5*h*g1(4)) f3(X0(1)+0.5*h*g1(1),X0(2)+0.5*h*g1(2)) f4(X0(1)+0.5*h*g1(1),X0(2)+0.5*h*g1(2))];
g3=[f1(X0(3)+0.5*h*g2(3)) f2(X0(4)+0.5*h*g2(4)) f3(X0(1)+0.5*h*g2(1),X0(2)+0.5*h*g2(2)) f4(X0(1)+0.5*h*g2(1),X0(2)+0.5*h*g2(2))];
g4=[f1(X0(3)+h*g3(3)) f2(X0(4)+h*g3(4)) f3(X0(1)+h*g3(1),X0(2)+h*g3(2)) f4(X0(1)+h*g3(1),X0(2)+h*g3(2))];
X1=X0+(h/6)*(g1+2*g2+2*g3+g4);
t1=t0+h;
SOL=[SOL;t1 X1];
X0=X1;
t0=t1;
end
plot(SOL(:,2),SOL(:,3))
axis('equal')
Attachments
Last edited: