RK4 orbit problem

5 views (last 30 days)
Lawrence Coleman
Lawrence Coleman on 26 Jul 2019
Edited: James Tursa on 30 Jul 2019
% 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.

Accepted Answer

Walter Roberson
Walter Roberson on 26 Jul 2019
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.
  3 Comments
Jim Riggs
Jim Riggs on 29 Jul 2019
Edited: Jim Riggs on 29 Jul 2019
"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
Lawrence Coleman on 29 Jul 2019
thanks. I forgot to add the minus signs. You are a saint.

Sign in to comment.

More Answers (1)

James Tursa
James Tursa on 29 Jul 2019
Edited: James Tursa on 30 Jul 2019
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.
  4 Comments
James Tursa
James Tursa on 29 Jul 2019
Edited: James Tursa on 29 Jul 2019
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
James Tursa on 29 Jul 2019
Edited: James Tursa on 29 Jul 2019
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).

Sign in to comment.

Categories

Find more on Particle & Nuclear Physics in Help Center and File Exchange

Community Treasure Hunt

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

Start Hunting!