MATLAB: Pendulum using Runge Kutta

errorsMATLABpendulumrunge kutta

Good Afternoon. I'm trying to run a code to simulate a simple pendulum but I am trying to run it using Runge Kutta, but I am getting errors. I have run it in the past using Euler and Verlet methods but the Runge Kutta method when put into code is confusing me. I want to chart theta vs time as well as multiple initial omega vs displacement. I don't know how to put Runge Kutta into code and what is in the code now is my best attempt. Thanks for the help.
%pendulrk - Simple Pendulum using Runge Kutta Method
clear all; help pendulrk;
%Initial Values
theta0 = input('Enter initial angle (in degrees): ');
theta = theta0*pi/180; %Convert angle to radians
omega = 0:2:20;
omega = 2*pi.*omega/60;
%Physical Constants
g_over_L = 1;   %The constant g/L
time = 0;       %Initial time
irev = 0;       %Used to count the number of reversals
tau = input('Enter time step: ');
%Loop over desired number of steps with given time step
N = input('Enter number of time steps: ');
for j = 1:length(omega)
    Om = omega(j);
for i=1:N
    %Record angle and time for plotting
    t_plot(i) = time;
    th_plot(i) = theta*180/pi; %Convert angle to degrees
    time = time + tau;
    
    %Record omega and theta for phase space plotting
    displace_plot(i,j) = theta;
    omega_plot(i,j) = Om;
    
    %Compute new position and velocity using Runge-Kutta Method
    accel = -g_over_L*sin(theta); %Gravitational Acceleration
    theta_old = theta;          %Save Previous Angle
    k1 = f(time(i),theta(i));
    k2 = f(time(i)+0.5*tau,theta(i)+0.5*tau*k1);
    k3 = f(time(i)+0.5*tau,theta(i)+0.5*tau*k2);
    k4 = f(time(i)+tau,theta(i)+k4*tau);
    theta(i+1)=theta(i)+(1/6)*tau*(k1+2*k2+2*k3+k4);
    
    theta = theta + tau*Om;     %Euler Method
    Om = Om + tau*accel;
end
clear Om;
end
%Graph the oscillations as theta versus time
clf; figure(gcf); %Clear and forward figure window
figure(1)
plot(t_plot,th_plot, '+');
grid on
xlabel('Time (s)'); ylabel('\theta (degrees)');
%Graph omega versus theta
figure(2)
plot(displace_plot,omega_plot);
grid on
xlim([-pi,pi]);
xlabel('\theta (rads)'); ylabel('\omega (rad/s)');

Best Answer

Yes, there are a few things wrong with this code.
James Tursa has identified one obvious error.
Another problem is with the loop indexing. If N is the number of steps, you need to pre-allocate all of your arrays to size N+1 (notice that you assign a value to theta(i+1) inside a loop that goes from 1 to N)
But the biggest problem is that the code contains a serious conceptual error. The Runge-Kutta method provides a numerical estimate for the integration of a function, but you are attempting to produce a double integration. To do this, you need to apply the RK method twice, just like you did with the Euler method : Om = Om + tau*accel and theta = theta + tau*Om. This represents two integrations. Om is the integral of accel, and theta is the integral of Om.
If function f() returns the pendulum acceleration, then the final result of this RK calculation will be omega, not theta.
Related Question