MATLAB Answers

Can't generate plots of orbit using RK4

2 ビュー (過去 30 日間)
Kyle Broder
Kyle Broder 2016 年 8 月 11 日
編集済み: James Tursa 2016 年 8 月 15 日
I'm considering a system of differential equations given by dR/dt = U and dU/dt = (L-GM)/R^2. I'm trying to use RK4 to solve this system and generate some plots. My code is given below, however I'm not generating any plots. Is therefore something wrong with my RK4 code? Note that the code does run.
function [x,t] =nma_RK4_4( )
delT = 0.001; %time grid
t = linspace(0,5,round(5/delT)); %time
N = length(t);
x = zeros(2,N);
x(1,1) = 100000; %initial conditions
x(2,1) = 100000;
G = 6.67*power(10,-11);
M = 1.5*1.989*power(10,41);
L = x(1,1)*x(2,2);
theta = -x(2,1)/(power(x(1,1),2));
xvar = x(1,1)*cos(theta);
yvar = x(1,1)*sin(theta);
for i = 1:(N-1)
k1 = delT * f( t(i) , x(:,i));
k2 = delT * f( t(i)+1/2*delT , x(:,i)+1/2*k1);
k3 = delT * f( t(i)+1/2*delT , x(:,i)+1/2*k2);
k4 = delT * f( t(i)+delT , x(:,i)+k3);
x(:,i+1) = x(:,i)+1/6*(k1+2*k2+2*k3+k4);
end
function r=f(t,x)
r=[x(2)
(L-G*M)/(power(x(1),2))];
end
figure()
plot(xvar,yvar)
end
Essentially, what I want to do is solve the above system of ODEs for R and then use that value obtained by using RK4 to then plot the result. I want to plot x against y however, and this is obtained from setting x = Rcos(theta) and y=Rsin(theta). As can be seen this is not working for me.
I don't want to plot just a point, I want something like this, https://arxiv.org/pdf/1503.05861.pdf, see the first graph on page 8. Note that the plot I've given as an example in that link does not look like what I have, given that it is a different system of ODEs, but the type of plot it what I want to generate.

  0 件のコメント

サインインしてコメントする。

回答 (1 件)

James Tursa
James Tursa 2016 年 8 月 11 日
編集済み: James Tursa 2016 年 8 月 11 日
Try moving your f code outside of your function code. E.g.,
:
figure()
plot(xvar,yvar)
end
function r=f(t,x)
r=[x(2)
(L-G*M)/(power(x(1),2))];
end
Also, you are plotting xvar and yvar, but these variables are never updated by your RK4 code.

  2 件のコメント

Kyle Broder
Kyle Broder 2016 年 8 月 11 日
Even if I change the code to
function [x,t] =nma_RK4_4( )
delT = 0.001; %time grid
t = linspace(0,5,round(5/delT)); %time
N = length(t);
x = zeros(2,N);
x(1,1) = 100000; %initial conditions
x(2,1) = 100000;
G = 6.67*power(10,-11);
M = 1.5*1.989*power(10,41);
L = x(1,1)*x(2,2);
theta = -x(2,1)/(power(x(1,1),2));
xvar = x*cos(theta);
yvar = x*sin(theta);
for i = 1:(N-1)
k1 = delT * f( t(i) , x(:,i));
k2 = delT * f( t(i)+1/2*delT , x(:,i)+1/2*k1);
k3 = delT * f( t(i)+1/2*delT , x(:,i)+1/2*k2);
k4 = delT * f( t(i)+delT , x(:,i)+k3);
x(:,i+1) = x(:,i)+1/6*(k1+2*k2+2*k3+k4);
figure()
plot(xvar,yvar)
end
function r=f(t,x)
r=[x(2)
(L-G*M)/(power(x(1),2))];
end
end
The plot still isn't generated.
James Tursa
James Tursa 2016 年 8 月 15 日
Try this code to get a plot based on the "x" result from your RK4 loop. Although even with this code there is an apparent issue based on the resulting plot. Still not sure what your xvar and yvar are supposed to be and how they are connected to the problem. But at least you have something you can step through and debug now and get a plot. (Your current code has the figure() and plot() inside your for loop, so it is attempting to generate some 5000 separate plots!)
function [x,t] =nma_RK4_4( )
delT = 0.001; %time grid
t = linspace(0,5,round(5/delT)); %time
N = length(t);
x = zeros(2,N);
x(1,1) = 100000; %initial conditions
x(2,1) = 100000;
G = 6.67*power(10,-11);
M = 1.5*1.989*power(10,41);
L = x(1,1)*x(2,2);
theta = -x(2,1)/(power(x(1,1),2));
xvar = x*cos(theta);
yvar = x*sin(theta);
for i = 1:(N-1)
k1 = delT * f( t(i) , x(:,i));
k2 = delT * f( t(i)+1/2*delT , x(:,i)+1/2*k1);
k3 = delT * f( t(i)+1/2*delT , x(:,i)+1/2*k2);
k4 = delT * f( t(i)+delT , x(:,i)+k3);
x(:,i+1) = x(:,i)+1/6*(k1+2*k2+2*k3+k4);
end
figure()
plot(x(1,:),x(2,:))
function r=f(t,x)
r=[x(2)
(L-G*M)/(power(x(1),2))];
end
end

サインインしてコメントする。

Community Treasure Hunt

Find the treasures in MATLAB Central and discover how the community can help you!

Start Hunting!

Translated by