Asked by Lawrence Coleman
on 26 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.

Answer by Walter Roberson
on 26 Jul 2019

Accepted Answer

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
on 29 Jul 2019

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
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
on 29 Jul 2019

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

Sign in to comment.

Answer by James Tursa
on 29 Jul 2019

Edited by 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.

Lawrence Coleman
on 29 Jul 2019

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
on 29 Jul 2019

Sign in to comment.

Opportunities for recent engineering grads.

Apply Today
## 0 Comments

Sign in to comment.