MATLAB Answers

Plotting a Trajectory using Euler's Method

50 ビュー (過去 30 日間)
Kyle Broder
Kyle Broder 2016 年 8 月 14 日
I'm solving the system of ODEs given by
R' = U
U' = Lz^2/R^2 +GM/R^2
using Euler's method. This is no problem, I have a functioning code given by
function euler
Tsim = 10; % simulate over Tsim seconds
h = 0.01; % step size
N = Tsim/h; % number of steps to simulate
x = zeros(2,N);
t = zeros(1,N);
G = 6.67*power(10,-11);
M = 100;
Lz = 20004444444;
R_1 = 10000;
U_1 = 20000;
x(1,1) = U_1; %IC_ODE_1
x(2,1) = (G*M/(power(R_1,2))+power(Lz,2)/(power(R_1,3))); %IC_ODE_2
for k=1:N-1
t(k+1)= t(k) + h;
x(:,k+1) = x(:,k) + h* Orbit(t(k),x(:,k));
end
%Compare with ode45
x0 = [1000;0];
[T,Y] = ode45(@Orbit,[0 Tsim],x0);
plot(t,x(1,:)) %Results from Euler algorithm
hold on
plot(T,Y(:,1),'^') %Results from ode45
legend('Euler','ode45')
end
function dx= Orbit(t,R)
G = 6.67*power(10,-11);
M = 100;
Lz = 200;
% Equations of motion
dx = zeros(2,1);
dx(1)=R(2);
dx(2)= (G*M/(power(R(1),2))+power(Lz,2)/(power(R(1),3)));
end
This generates plots fine. I want to plot the trajectory however. To do this, I need to be able to generate a value for two new variables. Labelling these by xvar and yvar, I essentially want to use the result of the Euler method code (which gives me a radius at time t) to give me trajectory in polar coordinates. Thus I require something like
theta(:,k+1) = theta(:,k) + h*Angle(t(k),theta(:,k));
xvar(:,k+1) = x(:,k+1)*cos(theta(:,k+1));
yvar(:,k+1) = x(:,k+1)*sin(theta(:,k+1));
Where theta is given by the ODE dtheta/dt = Lz^2/R^2. Therefore, I'm assuming I require something like this
function da = Angle(t,R)
Lz = 200;
da = Lz/(power(R(1),2))
end
But theta is given in terms of only one ODE, the one above, but that R(1) is obtained from solving the system above. I just want to be able to plot xvar against yvar with theta varying as it should.

  0 件のコメント

サインイン to comment.

回答 (1 件)

Swathik Kurella Janardhan
Swathik Kurella Janardhan 2016 年 8 月 16 日
You can calculate the initial theta(:,k) for an initial value of k by ODE
dtheta/dt = Lz^2/R^2
and then use this value for further calculations of theta(:,k+1) by
theta(:,k+1) = theta(:,k) + h*Angle(t(k),theta(:,k));

  0 件のコメント

サインイン to comment.

サインイン してこの質問に回答します。


Translated by