Generate a Runge-Kutta Integrator solver for the orbital motion equation, when given the initial conditions. Execute for 500 seconds.
조회 수: 13 (최근 30일)
이전 댓글 표시
The Problem statement is as followed.
Using the constant for km/s. and the inital conditions
km
and
km/s
The orbital motion equation is:
integrate using runge kutta RK4 for an h =10 seconds and an end time of 500 seconds
So I am working on this problem and I am stuck on the state variable transformation to make it a first order system and subsequently applying the Runge Kutta integrator. I am not sure how to work through it seeing that we are dealing with vectors instead of a scalar . I did the Ode45 approach as well as an attempt as well but i am not certain if i am correct. Any pointers would be greatly appreciated.
function xdot = exfunc(t,x)
mu = 398600.64;
xdot(1) = x(2);
xdot(2) = (-mu/((x(1))^(3)))*x(1);
xdot(3) = x(4);
xdot(4) = (-mu/((x(3))^(3)))*x(3);
xdot(5) = x(6);
xdot(6) = (-mu/((x(5))^(3)))*x(5);
xdot = xdot';
end
clear
clc
[t,x] = ode45('exfunc',[0 1],[5371.844151186472, 3.230645742388105,...
-4141.062031065303, 3.522344029484922, 460.1400917227622,...
-5.981911152962826])
plot(t,x(:,[1,3,5]))
R_end = x(end,[1,3,5])
V_end = x(end,[2,4,6])
댓글 수: 2
채택된 답변
James Tursa
2021년 3월 15일
편집: James Tursa
2021년 3월 15일
First, I find it easier to keep the position elements together, and the velocity elements together, in the state vector. So I would make this definition:
x(1) = x position
x(2) = y position
x(3) = z position
x(4) = x velocity
x(5) = y velocity
x(6) = z velocity
Then in your derivative function you can work with the vectors more easily. Be sure to use the full position vector for calculating r. And I would advise appending comments with units on all lines with constants. E.g.,
function xdot = exfunc(t,x)
mu = 398600.64; % km^3/s^2
r = x(1:3);
rdot = x(4:6);
rmag = norm(r);
rdotdot = -(mu/rmag^3)*r;
xdot = [rdot;rdotdot];
end
And of course that means you would need to rearrange your initial conditions that you feed to ode45:
r0 = [5371.844151186472; -4141.062031065303; 460.1400917227622]; % km
v0 = [3.230645742388105; 3.522344029484922; -5.981911152962826]; % km/s
t0 = 0; % s
tend = 500; % s
[t,x] = ode45(@exfunc,[t0 tend],[r0;v0]);
plot(t,x(:,1:3))
R_end = x(end,1:3)
V_end = x(end,4:6)
댓글 수: 0
추가 답변 (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!