RK4 orbit problem

% Parameters
h = 0.0027;      %step size
tfinal = 5;  % t=[0,tfinal]
GM = 4*pi^2 ;
% Initial Conditions
vx = zeros(1,1825);
vy = zeros(1,1825);
x = zeros(1,1825);
y = zeros(1,1825);

x(1) = 0;
y(1) = 1;
vy(1) = 0;
vx(1) = 1;
t = zeros(1,1825);


% Define the ODE function handle
dvxdt=@(t,x,vx) GM*x./(x.^2+y.^2).^3;
dvydt=@(t,y,vy) GM*y./(x.^2+y.^2).^3;
dxdt=@(t,x,vx) vx;
dydt=@(t,y,vy) vy;


%Runge Kutta Loop
for i=1:ceil(tfinal/h)
    t(i+1) = t(i) + h;
    
    
    k1vx = dvxdt(t(i)     ,x(i)     ,vx(i));
    k1vy = dvydt(t(i)     ,y(i)     ,vy(i));
    k1x = dxdt(t(i)     ,x(i)     ,vx(i));
    k1y = dydt(t(i)     ,y(i)     ,vy(i));
    
    k2vx = dvxdt(t(i)+0.5*h,x(i)+0.5*h*k1vx,vx(i)+0.5*h*k1vx);
    k2vy = dvydt(t(i)+0.5*h,y(i)+0.5*h*k1vy,vy(i)+0.5*h*k1vy);
    k2x = dxdt(t(i)+0.5*h,x(i)+0.5*h*k1x,vx(i)+0.5*h*k1x);
    k2y = dydt(t(i)+0.5*h,y(i)+0.5*h*k1y,vy(i)+0.5*h*k1y);
    
    k3vx = dvxdt(t(i)+0.5*h,x(i)+0.5*h*k2vx,vx(i)+0.5*h*k2vx);
    k3vy = dvydt(t(i)+0.5*h,y(i)+0.5*h*k2vy,vy(i)+0.5*h*k2vy);
    k3x = dxdt(t(i)+0.5*h,x(i)+0.5*h*k2x,vx(i)+0.5*h*k2x);
    k3y = dydt(t(i)+0.5*h,y(i)+0.5*h*k2y,vy(i)+0.5*h*k2y);
    
    k4vx = dvxdt(t(i)+    h,x(i)    +h*k3vx,vx(i)     +h*k3vx);
    k4vy = dvydt(t(i)+    h,y(i)    +h*k3vy,vy(i)     +h*k3vy);
    k4x = dxdt(t(i)+    h,x(i)    +h*k3x,vx(i)     +h*k3x);
    k4y = dydt(t(i)+    h,y(i)    +h*k3y,vy(i)     +h*k3y);
    
    vx(i+1)=vx(i)+h/6*(k1vx+2*k2vx+2*k3vx+k4vx);
    vy(i+1)=vy(i)+h/6*(k1vy+2*k2vy+2*k3vy+k4vy);
    x(i+1)=x(i)+h/6*(k1x+2*k2x+2*k3x+k4x);
    y(i+1)=y(i)+h/6*(k1y+2*k2y+2*k3y+k4y);
end
%plot the results
plot(x,y)

it says the left and right sides have a different number of elements and I dont know how to fix it.

Answer 1

[Walter Roberson]

dvxdt=@(t,x,vx) GM*x./(x.^2+y.^2).^3;

that accepts vx as a parameter but does not use it.

It does not have any parameter named y but uses y. In that case, it is going to use the variable value that the variable y had at the time the anonymous function was created, which comes from

y = zeros(1,1825);

This is a vector, and even though it only contributes 0 at all components, it makes the result of dvxdt be a vector. That becomes a problem because it means that

    k1vx = dvxdt(t(i)     ,x(i)     ,vx(i));

gives a vector result, so

vx(i+1)=vx(i)+h/6*(k1vx+2*k2vx+2*k3vx+k4vx);

has a vector on the right hand side, but designates a scalar output location.

[Lawrence Coleman]

I have edited the code and I now get an output. The output does not look like an orbit though. messing with the initial velocity did not help.

% Parameters
h = 0.0027;      %step size
tfinal = 5;  % t=[0,tfinal]
GM = 4*pi^2 ;
% Initial Conditions
vx = zeros(1,1825);
vy = zeros(1,1825);
x = zeros(1,1825);
y = zeros(1,1825);

x(1) = 0;
y(1) = 1;
vy(1) = 0;
vx(1) = 2;
t = zeros(1,1825);


% Define the ODE function handle
dvxdt=@(t,x,y) GM*x./(x.^2+y.^2).^3;
dvydt=@(t,y,x) GM*y./(x.^2+y.^2).^3;
dxdt=@(t,x,vx) vx;
dydt=@(t,y,vy) vy;


%Runge Kutta Loop
for i=1:ceil(tfinal/h)
    t(i+1) = t(i) + h;
    
    
    k1vx = dvxdt(t(i)     ,x(i)     ,y(i));
    k1vy = dvydt(t(i)     ,y(i)     ,x(i));
    k1x = dxdt(t(i)     ,x(i)     ,vx(i));
    k1y = dydt(t(i)     ,y(i)     ,vy(i));
    
    k2vx = dvxdt(t(i)+0.5*h,x(i)+0.5*h*k1vx,y(i)+0.5*h*k1vx);
    k2vy = dvydt(t(i)+0.5*h,y(i)+0.5*h*k1vy,x(i)+0.5*h*k1vy);
    k2x = dxdt(t(i)+0.5*h,x(i)+0.5*h*k1x,vx(i)+0.5*h*k1x);
    k2y = dydt(t(i)+0.5*h,y(i)+0.5*h*k1y,vy(i)+0.5*h*k1y);
    
    k3vx = dvxdt(t(i)+0.5*h,x(i)+0.5*h*k2vx,y(i)+0.5*h*k2vx);
    k3vy = dvydt(t(i)+0.5*h,y(i)+0.5*h*k2vy,x(i)+0.5*h*k2vy);
    k3x = dxdt(t(i)+0.5*h,x(i)+0.5*h*k2x,vx(i)+0.5*h*k2x);
    k3y = dydt(t(i)+0.5*h,y(i)+0.5*h*k2y,vy(i)+0.5*h*k2y);
    
    k4vx = dvxdt(t(i)+    h,x(i)    +h*k3vx,y(i)     +h*k3vx);
    k4vy = dvydt(t(i)+    h,y(i)    +h*k3vy,x(i)     +h*k3vy);
    k4x = dxdt(t(i)+    h,x(i)    +h*k3x,vx(i)     +h*k3x);
    k4y = dydt(t(i)+    h,y(i)    +h*k3y,vy(i)     +h*k3y);
    
    vx(i+1)=vx(i)+h/6*(k1vx+2*k2vx+2*k3vx+k4vx);
    vy(i+1)=vy(i)+h/6*(k1vy+2*k2vy+2*k3vy+k4vy);
    x(i+1)=x(i)+h/6*(k1x+2*k2x+2*k3x+k4x);
    y(i+1)=y(i)+h/6*(k1y+2*k2y+2*k3y+k4y);
end
%plot the results
plot(x,y)

[Jim Riggs]

"messing with the initial vlocity did not help"

The initial velocity is clearly way too large. Didi you try reducing the velocity?

Keep reducing it until the curve closes.

(Also, check the sign on the dvxdt and dvydt to make sure that they are pointed in the correct direction. Draw a picture - that always helps)

[Lawrence Coleman]

thanks. I forgot to add the minus signs. You are a saint.

Answer 2

[James Tursa]

Couple of things:

1) You have picked the most convoluted method of coding your equations. Having four separate variables for your state and switching the order of x and y in the arguments causes your code to bloat and makes it very hard to read ... so much that I am not taking the time to review your RK4 stuff for correctness at this time. You would be much better off defining a 4-element vector for your state with the elements being x,y,vx,vy. Then your code will need only one derivative function and you will need only one set of k's calculated (they will be vectors).

2) Not sure what kind of orbit problem has GM = 4*pi^2 as the central force parameter (where does this come from?), but if it is supposed to be a "proportional to inverse r^2" central force orbit problem then you need the equivalent of r^3 in the denominator of your derivative functions, not r^6.

E.g., this (effectively (r^2)^3 = r^6 in the denominator):

dvxdt=@(t,x,y) GM*x./(x.^2+y.^2).^3;
dvydt=@(t,y,x) GM*y./(x.^2+y.^2).^3;

should be this instead (effectively r^3 in the denominator):

dvxdt=@(t,x,y) -GM*x./sqrt(x.^2+y.^2).^3;
dvydt=@(t,y,x) -GM*y./sqrt(x.^2+y.^2).^3

I've added minus signs as well, since I am guessing you want this to be an attractive central force, not repulsive.

[Lawrence Coleman]

G = 4pi^2 when I choose the units of solar masses and AU for distance. Keeps the numbers from getting messy. Thanks for catching my mistake about the r^6.

[Lawrence Coleman]

also, can you elaborate on how I can make this work with only one derivative function? F=ma is a second order ode after all. Im pretty new to this MATLAB stuff so I am sorry if my methods do not seem efficient.

[James Tursa]

To align with MATLAB documentation on the ode methods, I am going to use y as the state vector with the following definitions (this is not code ... this is just definitions):

y(1) = x
y(2) = y
y(3) = vx = xdot
y(4) = vy = ydot

With that basic definition of the y vector, you define a single derivative function:

dy = @(t,y) [ y(3); y(4); -GM*y(1)/sqrt(y(1)^2+y(2)^2)^3; -GM*y(2)/sqrt(y(1)^2+y(2)^2)^3 ]

The (t,y) argument list matches the signature used by MATLAB's ode function suite, so it could easily be used in those calls as well.

Then the pre-allocation is simply this, where each column of y represents a single state vector at a specific time:

y = zeros(4,1825);

And the initial conditions are:

y0 = [ 0; 1; 2; 0 ];
y(:,1) = y0;

And inside your for-loop would be something like this for the k's and state vector update:

    k1 = dy( t(i)        , y(:,i)            );
    k2 = dy( t(i) + 0.5*h, y(:,i) + 0.5*h*k1 );
    k3 = dy( t(i) + 0.5*h, y(:,i) + 0.5*h*k2 );
    k4 = dy( t(i) +     h, y(:,i)      +h*k3 );
    y(:,i+1) = y(:,i) + h/6*( k1 + 2*k2 + 2*k3 + k4 );

See how much simpler this looks and how much easier it is to read? All of the individual equations involving x, y, vx, and vy are replaced with equations involving the single state vector y(:,i) etc.

For plotting you simply pick off the appropriate rows of y. E.g., all of the x values are y(1,:), all of the y values are y(2,:), etc.

[James Tursa]

P.S. I think you will find that for your initial conditions, the orbit is elliptical and comes close enough to the origin that your stepsize is too big. E.g., circular velocity would by sqrt(GM/r) = 2*pi but you are starting with a velocity of 2 which will result in an elliptical orbit. You will probably need to drop your stepsize by at least a factor of 10 to get good results (or have an adaptive stepsize).