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

  • MATLAB
  • Thread starter saeede-
  • Start date
  • #1
6
0
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:
Matlab:
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.>
 
Last edited by a moderator:

Answers and Replies

  • #2
DrClaude
Mentor
7,584
3,946
What is the error you get?
 
  • #3
6
0
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);
 
  • #4
DrClaude
Mentor
7,584
3,946
Line 16 looks like it should be a comment, not a statement. Fixing that will probably make all those errors go away.
 
  • Like
Likes pasmith and saeede-
  • #5
6
0
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 :)
 
  • #6
pasmith
Homework Helper
2,004
637
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:
Matlab:
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 [itex]x[/itex] direction and your initial velocity is also in the [itex]x[/itex] direction. Thus the angular velocity is zero, and the satellite will remain on the [itex]x[/itex] axis.
Code:
[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: [itex]5.972 \times 10^{24} + 50 = 5.972 \times 10^{24}[/itex]. You might as well set [itex]m = 0[/itex].

I would suggest scaling distances with [itex]r + R_0[/itex] and times with [itex]((r+R_0)^3/G(M_e+m))^{1/2}[/itex] so that in scaled units [tex]
\ddot{\mathbf{r}} = - \frac{1}{r^2} \hat{\mathbf{r}}
[/tex] with [itex]\mathbf{r}(0) = (1,0,0)[/itex]. Velocities then scale as the length unit divided by the time unit, so the initial scaled velocity is [tex]\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).[/tex] (The problem is entirely two-dimensional, so you could get rid of the [itex]z[/itex] 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 [itex]2\pi[/itex] scaled time units. So your maximum time doesn't need to go beyond 10 scaled units, rather than the [itex]3.15 \times 10^7\,\mathrm{s}[/itex] you have allowed for.
 
  • #7
6
0
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 [itex]x[/itex] direction and your initial velocity is also in the [itex]x[/itex] direction. Thus the angular velocity is zero, and the satellite will remain on the [itex]x[/itex] axis.


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

I would suggest scaling distances with [itex]r + R_0[/itex] and times with [itex]((r+R_0)^3/G(M_e+m))^{1/2}[/itex] so that in scaled units [tex]
\ddot{\mathbf{r}} = - \frac{1}{r^2} \hat{\mathbf{r}}
[/tex] with [itex]\mathbf{r}(0) = (1,0,0)[/itex]. Velocities then scale as the length unit divided by the time unit, so the initial scaled velocity is [tex]\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).[/tex] (The problem is entirely two-dimensional, so you could get rid of the [itex]z[/itex] 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 [itex]2\pi[/itex] scaled time units. So your maximum time doesn't need to go beyond 10 scaled units, rather than the [itex]3.15 \times 10^7\,\mathrm{s}[/itex] 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.
 

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

  • Last Post
Replies
4
Views
6K
Replies
2
Views
24K
Replies
5
Views
8K
  • Last Post
Replies
1
Views
3K
  • Last Post
Replies
7
Views
6K
Replies
5
Views
2K
  • Last Post
Replies
0
Views
2K
Top