Pendulum using Runge Kutta

66 views (last 30 days)
Nicholas Davis
Nicholas Davis on 3 Nov 2020
Edited: Jim Riggs on 6 Nov 2020
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)');
  4 Comments
James Tursa
James Tursa on 3 Nov 2020
Well, delete it since it is screwing up your RK4 code.
Nicholas Davis
Nicholas Davis on 3 Nov 2020
The error still remains, though I have a feeling I am missing something.

Sign in to comment.

Accepted Answer

Jim Riggs
Jim Riggs on 4 Nov 2020
Edited: Jim Riggs on 4 Nov 2020
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.
  13 Comments
Nicholas Davis
Nicholas Davis on 4 Nov 2020
That's it! Thank you Jim. I really appreciate it. I'm sure I could figure this out, but to save myself the trouble, how do I run the code using multiple initial omega, so omega = 0:2:20 instead of just 0? I get an error when I get to this line:
omega(1)=omega0*pi/180;
Jim Riggs
Jim Riggs on 4 Nov 2020
Edited: Jim Riggs on 6 Nov 2020
You have to think carefully about what you are doing. If you set omega = 0:2:20, you have assigned 11 values to omega. But omega is being used as the time history vector which is initialized to N+1 zero values, where N is the number of time steps. So, what you want to do is define a set of initial values omega0 = 0:2:20.
I supposethat this means that you want to run the same case 11 times, using a different omega0 each time.
This means that you will want to save 11 time histories of theta and omega. So let "no" be the number of omega initial conditions.
no = length(omega0);
By using the length function to find the number of omega0 values, you simply set the desired values of omega0 and the code will do the rest. For example, you can specify a list of the values you want, such as
omega0 = [-2, 0, 2, 5, 12]
and this will work too.
Now the vectors that you use to save the data will have to be 2-dimentional arrays, so you want to pre-allocate your vectors as:
time = zeros(N+1, no);
theta = zeros(N+1, no);
omega = zeros(N+1, no);
Add a loop around the code to run "no" times, using a different omega0 value each time
omega0 = [-2, 0, 2, 5, 12] ; % specify whatever you like. You can also make this a user input
no = length(omega0);
time = zeros(N+1, no);
theta = zeros(N+1, no);
omega = zeros(N+1, no);
for j=1:no % loop over multiple omega0 values
time(1,j)= 0;
theta(1,j) = theta0*pi/180;
omega(1,j) = omega0(j) % index j corresponds to omega0
...
STATE(1) = theta(1,j);
STATE(2) = omega(1,j);
% now run the i loop
for i=1:N
...
...
% results are now saved in a 2-dimentional array
time(i+1,j) = time(i) + tau;
theta(i+1,j) = STATE(1);
omega(i+1,j) = STATE(2);
end % of i loop
...
end % of j loop
After you have run this, you plot individual columns from the saved data, e.g.
plot( time(:,1), theta(:,1)...) % this means plot all rows (:) from column 1
hold on;
plot( time(:,2), theta(:,2)...) % plot all rows (:) from column 2
... (etc)

Sign in to comment.

More Answers (0)

Categories

Find more on General Applications in Help Center and File Exchange

Products


Release

R2020b

Community Treasure Hunt

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

Start Hunting!