how to replace the ode45 method with the runge-kutta in this matlab?

Question

I tried everything and looked everywhere but can't find any solution for my question.

clc
clear all 


%% Solving the Ordinary Differential Equation 
G = 6.67408e-11; %Gravitational constant 
M = 10; %Mass of the fixed object 
r = 1; %Distance between the objects 

tspan = [0 100000]; %Time Progression from 0 to 100000s 
conditions = [1;0]; %y0= 1m apart, v0=0 m/s

F=@(t,y)var_r(y,G,M,r);

[t,y]=ode45(F,tspan,conditions); %ODE solver algorithm

%%part1: Plotting the Graph 
% plot(t,y(:,1)); %Plotting the Graph 
% xlabel('time (s)') 
% ylabel('distance (m)')

%% part2: Animation of Results 
plot(0,0,'b.','MarkerSize', 40); 
hold on    %to keep the first graph 
for i=1:length(t) 
k = plot(y(i,1),0,'r.','MarkerSize', 12); 
pause(0.05); 
axis([-1 2 -2 2]) %Defining the Axis 
xlabel('X-axis') %X-Axis Label 
ylabel('Y-axis') %Y-Axis Label 
delete(k)
end 

function yd=var_r(y,G,M,r) %function of variable r 
g = (G*M)/(r + y(1))^2; 
yd = [y(2); -g]; 
end 

this is the code where I'm trying to replace the ode45 with the runge kutta method but its giving me errors. my runge kutta function:

function y = Runge_Kutta(f,x0,xf,y0,h)

n= (xf-x0)/h;
y=zeros(n+1,1);
x=(x0:h:xf);
y(1) = y0;


for i=1:n
  k1 = f(x(i),y(i));
  k2= f(x(i)+ h/2 , y(i) +h*(k1)/2);
  y(i+1) = y(i)+(h*k2);
end

plot(x,y,'-.M')
legend('RKM')
title ('solution of y(x)');
xlabel('x');
ylabel('y(x)')
hold on
end 

Answer 1

y(1) is the first scalar element in the data structure of y (this counts in column-first order). You want to generate in y a list of column vectors, as your ODE is a system with state dimension 2. Thus you need to generate y with that format, y=zeros(length(x0),n+1); and then address the list entries as matrix columns y(:,1)=x0 and the same modification in every place where you extract or assign a list entry.

---

Matlab introduce various short-cuts that, if used consequently, lead to contradictions (I think the script-hater rant (german) is still valid in large parts). Essentially, unlike in other systems, Matlab gives direct access to the underlying data structure of matrices. y(k) is the element of the underlying flat array (that is interpreted column-first in Matlab like in Fortran, unlike, e.g., Numpy where it is row-first).

Only the two-index access is to the matrix with its dimensions. So y(:,k) is the k-th matrix column and y(k,:) the k-th matrix row. The single-index access is nice for row or column vectors, but leads immediately to problems when collecting such vectors in lists, as these lists are automatically matrices.

Answer 2

Before converting your ode45( ) solution to manually written RK scheme, it doesn't even look like your ode45( ) solution is correct. It appears you have a gravitational problem set up where the initial velocity is 0 so a small object will simply fall into a large mass M on a line (rectilinear motion), and that is why you have scalar position and velocity.

Going with this assumption, r is something you should be calculating on the fly, not using as a fixed input to the derivative function. E.g., I would have expected something like this:

F=@(t,y)var_r(y,G,M); % get rid of r
    :

function yd=var_r(y,G,M) % function of current position y(1) and velocity y(2) 
g = (G*M)/y(1)^2; % gravity accel based on current position
yd = [y(2); -g]; % assumes y(1) is positive, so acceleration is negative
end 

The small object must start with a positive initial position for the derivative code to be valid as you have it written. As the small object falls into the large mass M, the above will only hold until it hits the surface or atmosphere of M. Or if you model M as a point mass, then this scheme will become increasingly difficult to integrate correctly because the acceleration becomes large without bound as the small mass gets very close to the point mass M. You would definitely need a variable step size approach in this case. The solution becomes invalid if it goes "through" mass M. In fact, once the speed gets too large the whole setup becomes invalid because of relativistic effects.

Maybe you could explain in more detail if your system is supposed to be set up this way, and what the purpose of the integration is. If it is really supposed to be a 2D or 3D problem, then more states need to be added.

For your manual Runge-Kutta code, you completely forgot to integrate the velocity so this is going to fail miserably. You need to carry a 2-element state from step to step, not a scalar as you are currently doing. E.g., something like this:

y=zeros(2,n+1); % 2-element state as columns of the y variable
x=(x0:h:xf);
y(:,1) = y0; % initial state is the first 2-element column
% change all the scalar y(i) to column y(:,i)
for i=1:n
  k1 = f(x(i),y(:,i));
  k2= f(x(i)+ h/2 , y(:,i) +h*(k1)/2);
  y(:,i+1) = y(:,i)+(h*k2);
end

plot(x,y(1,:),'-.M') % plot the position part of the solution

This is all assuming the f that gets passed in is the same F you have in your original code.