Numerical integration of the differential equation of motion of the two body problem
33 ビュー (過去 30 日間)
古いコメントを表示
I have to write a code that integrates the differential equation of motion of the 2-body problem numerically, starting from initial values of position and velocity in the three-dimensional space, using this equation:
Initial values of a Geostationary satellite.
I have created a model, integration and main code file.
Model: model.m
function xdot = model(x,u);
global mu r
xdot(1,1) = x(2);
xdot(2,1) = (-mu/r^3)*x(1)
Integration: rk4int.m
function x = rk4int(model, h, x, u)
k1 = h*feval(model, x, u);
k2 = h*feval(model, x+k1/2, u);
k3 = h*feval(model, x+k2/2, u);
k4 = h*feval(model, x+k3, u);
x = x + (k1 + 2*k2 + 2*k3 + k4)/6;
Main code:
Define Constants and initial conditions
global mu r v
mu = 398600; % km^3/s^2
R0 = [42164 0 0]
V0 = [3.0746 0 0]
r = norm(R0)
v = norm(V0)
t = 0
E = v*v/2 - mu/r
a = -mu/(2*E)
n = sqrt(mu/(a*a*a))
T = 2*pi/n %period
Define the parameters
stepsize = T/1000 % Integration step size
comminterval = 0.01 % Communications interval
EndTime = T % Duration of the simulation (final time)
i = 0 % Initialise counter for data storage
Initial states
u = 0
x = [42164,0,0]'
xdot = [3.0746,0,0]'
Time
for time = 0:stepsize:EndTime
if rem(time,comminterval)==0
i = i+1 % increment counter
tout(i) = time % store time
xout(i,:) = x % store states
xdout(i,:) = xdot % store state derivatives
end
xdot = model(x,u)
x = rk4int('model',stepsize,x,u)
end
However, when I try to run the main code, it stops at the ''x = rk4int('model',stepsize,x,u)'' line gives me the following error:
What is the issue? How can I fix it?
0 件のコメント
採用された回答
James Tursa
2020 年 11 月 24 日
編集済み: James Tursa
2020 年 11 月 24 日
Your biggest problem is that you don't carry enough states in your derivative function. Your ODE is a 3D 2nd order equation, so that means you will need 3*2 = 6 states to carry. But your derivative function only calculates two scalar derivatives when it should be calculating six scalar derivatives. Also, you need to calculate r from the current state vector ... not pass in a constant for this. So, first let's decide what goes into your six element state vector x:
x(1) = x position
x(2) = y position
x(3) = z position
x(4) = x velocity (i.e., xdot)
x(5) = y velocity (i.e., ydot)
x(6) = z velocity (i.e., zdot)
Given those definitions, your derivative function should look something like this instead:
function xdot = model(x,u);
global mu
xdot = zeroes(size(x));
xdot(1:3) = x(4:6);
r = norm(x(1:3));
xdot(4:6) = (-mu/r^3)*x(1:3);
end
Then for initial state you would simply have
R0 = [42164 0 0];
V0 = [0 3.0746 0]; % <-- changed! The velocity needs to be in the ydot element for circular orbit
x = [R0,V0]';
2 件のコメント
BHOOMIKA
2022 年 11 月 23 日
Hello ! I have the same question and I want to solve it using 0de45, could you help me with the code when we use that function.
その他の回答 (0 件)
参考
カテゴリ
Help Center および File Exchange で Numerical Integration and Differential Equations についてさらに検索
Community Treasure Hunt
Find the treasures in MATLAB Central and discover how the community can help you!
Start Hunting!