Question about how solve 2 body problem in orbit using MATLAB ode45

[saeede-]

hi guys, I'm trying to write a program in MATLAB to solve and plot equation of motion of 2 body problem but it errors as i don't know what it says!
do you know how to help me? please!

The main equation is r(double dot)=-(mu)*r(hat)/r^3
The code that i wrote is:

clear all
close all
clc

R=6371e+3;
r0=3e+5;
p0=[R+r0;0;0;7000;0;0];

[t,p]=ode45(@sattelite,[0,3.15581e+07],p0);
xout=p(:,1);
yout=p(:,2);
zout=p(:,3);
plot3d(xout,yout,zout)

function dp=sattelite(t,p)
p=[x position;y position;z position;x velocity(vx);y velocity(vy);z velocity(vz)];
G=6.67e-11;
Me=5.972e+24;
m=50;
mu=(Me+m)*G;
R=6371e+3;
r0=3e+5;

r=sqrt(p(1)^2+p(2)^2+p(3)^2);
dp=zeros(6,1);
dp(1)=p(4);
dp(2)=p(5);
dp(3)=p(6);
dp(4)=-(mu/r^3)*p(1);
dp(5)=-(mu/r^3)*p(2);
dp(6)=-(mu/r^3)*p(3);
end

<Moderator's note: code tags added. Please use them.>

 

[DrClaude]

What is the error you get?

 

[saeede-]

What is the error you get?

all the errors are:

Unrecognized function or variable 'x'.

Error in T3>sattelite (line 16)
p=[x position;y position;z position;x velocity(vx);y velocity(vy);z velocity(vz)];

Error in odearguments (line 90)
f0 = feval(ode,t0,y0,args{:}); % ODE15I sets args{1} to yp0.

Error in ode45 (line 115)
odearguments(FcnHandlesUsed, solver_name, ode, tspan, y0, options, varargin);

Error in T3 (line 9)
[t,p]=ode45(@sattelite,[0,3.15581e+07],p0);

 

[DrClaude]

Line 16 looks like it should be a comment, not a statement. Fixing that will probably make all those errors go away.

 

[saeede-]

Line 16 looks like it should be a comment, not a statement. Fixing that will probably make all those errors go away.

Oh yes! thanks. all errors gone .although i wanted to plot a circle 3d but it shows a line :)

 

[pasmith]

Now that you have a program which runs, there are some other improvements you should consider.

hi guys, I'm trying to write a program in MATLAB to solve and plot equation of motion of 2 body problem but it errors as i don't know what it says!
do you know how to help me? please!

The main equation is r(double dot)=-(mu)*r(hat)/r^3
The code that i wrote is:

clear all
close all
clc

R=6371e+3;
r0=3e+5;
p0=[R+r0;0;0;7000;0;0];

You are seeing a line because your initial position is in the $x$ direction and your initial velocity is also in the $x$ direction. Thus the angular velocity is zero, and the satellite will remain on the $x$ axis.

[t,p]=ode45(@sattelite,[0,3.15581e+07],p0);
xout=p(:,1);
yout=p(:,2);
zout=p(:,3);
plot3d(xout,yout,zout)

function dp=sattelite(t,p)
p=[x position;y position;z position;x velocity(vx);y velocity(vy);z velocity(vz)];
G=6.67e-11;
Me=5.972e+24;
m=50;
mu=(Me+m)*G;
R=6371e+3;
r0=3e+5;

In terms of floating point arithmetic: $5.972 \times 10^{24} + 50 = 5.972 \times 10^{24}$. You might as well set $m = 0$.

I would suggest scaling distances with $r + R_0$ and times with $((r+R_0)^3/G(M_e+m))^{1/2}$ so that in scaled units $$\ddot{\mathbf{r}} = - \frac{1}{r^2} \hat{\mathbf{r}}$$ with $\mathbf{r}(0) = (1,0,0)$. Velocities then scale as the length unit divided by the time unit, so the initial scaled velocity is $$\dot{\mathbf{r}}(0) = \left(0, 7000 \left(\frac{R + r_0}{G(M_e + m)}\right)^{1/2}, 0\right) \approx (0, 0.889, 0).$$ (The problem is entirely two-dimensional, so you could get rid of the $z$ components of position and velocity and save the computer some work.)

In these units, a circular orbit with initial velocity (0,1,0) has a period of $2\pi$ scaled time units. So your maximum time doesn't need to go beyond 10 scaled units, rather than the $3.15 \times 10^7\,\mathrm{s}$ you have allowed for.

 

[saeede-]

Now that you have a program which runs, there are some other improvements you should consider.
You are seeing a line because your initial position is in the $x$ direction and your initial velocity is also in the $x$ direction. Thus the angular velocity is zero, and the satellite will remain on the $x$ axis.In terms of floating point arithmetic: $5.972 \times 10^{24} + 50 = 5.972 \times 10^{24}$. You might as well set $m = 0$.

I would suggest scaling distances with $r + R_0$ and times with $((r+R_0)^3/G(M_e+m))^{1/2}$ so that in scaled units $$\ddot{\mathbf{r}} = - \frac{1}{r^2} \hat{\mathbf{r}}$$ with $\mathbf{r}(0) = (1,0,0)$. Velocities then scale as the length unit divided by the time unit, so the initial scaled velocity is $$\dot{\mathbf{r}}(0) = \left(0, 7000 \left(\frac{R + r_0}{G(M_e + m)}\right)^{1/2}, 0\right) \approx (0, 0.889, 0).$$ (The problem is entirely two-dimensional, so you could get rid of the $z$ components of position and velocity and save the computer some work.)

In these units, a circular orbit with initial velocity (0,1,0) has a period of $2\pi$ scaled time units. So your maximum time doesn't need to go beyond 10 scaled units, rather than the $3.15 \times 10^7\,\mathrm{s}$ you have allowed for.

in my problem that i solve it , I would take different velocity and see how change orbit. and i should use m because it's mass of satellite.
I put the velocity that gave me in y component and become Ok !
Now , I don't know how to plot orbit around Earth . this said in different topic.