Can't generate plots of orbit using RK4

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.

Answer 1

[James Tursa]

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.

[Kyle Broder]

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]

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