7 views (last 30 days)

Show older comments

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)');

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.

Jim Riggs
on 4 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)

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

Start Hunting!