Pendulum using Runge Kutta
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)');
4 Comments
James Tursa
on 3 Nov 2020
Why is this line still in your code:
theta = theta + tau*Om; %Euler Method
If you are using RK4, you should delete this.
Nicholas Davis
on 3 Nov 2020
James Tursa
on 3 Nov 2020
Well, delete it since it is screwing up your RK4 code.
Nicholas Davis
on 3 Nov 2020
Accepted 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.
13 Comments
Nicholas Davis
on 4 Nov 2020
For the looping, consider the case where you ask for only 1 step. You still have two data points. The first one is the initial condition, the second one is the output of your integration step. So if N=1, your arrays need to contain 2 values. (i=1) is the initial condition, and (i=2) is the result of the one integration step. The loop goes from 1 to N (for a single step, this would be from 1 to 1), the number of values in each array is N+1.
You preallocate your arrays using a function like:
time = zeros(N+1,1) % Creates an array N+1 x 1 (a column vector) filled with zeros
As far as evaluating the double integral, it is a bit more complicated using higher-order methods like RK4.
First, consider the system that you are modeling: a simple pendulum. This is a kinetic system that is defined by it's position and velocity.
The position (theta) and velocity (omega) are the STATES of the system. You want to propagate these system states in time, and the method for doing that is to apply a numerical estimate of the integral of the state derivatives. So, for each state, we define it's derivative and apply numerical integration over the time step to get the new state at the new time.
I like to think of the states as a vector, for example STATE(1) = theta, STATE(2) = omega.
So I need to define a state derivative vector that corresponds to STATE, I will call it DSTATE.
DSTATE(1) = omega, and DSTATE(2) = accel
Note that omega is both a state and a state derivative. The procedure for working this problem is to;
1) define the initial state values (theta and omega at time=0)
2) Define the state derivative function: (accel, omega = f(time, STATE)
3) Run the integration over the desired time span
% pre-allocate your arrays
time = zeros(N+1,1);
theta= zeros(N+1,1);
omega= zeros(N+1,1);
% Set Initial conditions
time(1) = 0;
theta(1) = theta0 * pi()/180; % radians
omega(1) = omega0 * pi()/180; % radians/sec (assume input in deg/sec)
% Each STATE value will need to have its own k value
% Index 1 {STATE(1)} is theta
% Index 2 {STATE(2)} is omega
k1 = zeros(2,1);
k2 = zeros(2,1);
k3 = zeros(2,1);
k4 = zeros(2,1);
% initialize the state vector
STATE(1) = theta(1);
STATE(2) = omega(1);
% -------------------- RK4 step -----------------------------------------v
for i=1:N
% Save the STATE vector and time at the start of the integration time-step
time_in = time(i)
STATE_in = STATE;
% === RK4 stage 1 ===
t = time_in
DSTATE = f(t,STATE); % calculate the state derivatives
k1(1) = tau*DSTATE(1); % theta
k1(2) = tau*DSTATE(2); % omega
% update the states
STATE(1) = STATE_in(1) + k1(1) / 2;
STATE(2) = STATE_in(2) + k1(2) / 2;
% === RK4 stage 2 ===
t = time_in + tau/2
DSTATE = f(t,STATE); % calculate the state derivatives
k2(1) = tau*DSTATE(1); % theta
k2(2) = tau*DSTATE(2); % omega
% update the states
STATE(1) = STATE_in(1) + k2(1) / 2;
STATE(2) = STATE_in(2) + k2(2) / 2;
% === RK4 stage 3 === (Notice that t doesn't change this time)
DSTATE = f(t,STATE); % calculate the state derivatives
k3(1) = tau*DSTATE(1); % theta
k3(2) = tau*DSTATE(2); % omega
% update the states
STATE(1) = STATE_in(1) + k3(1) / 2;
STATE(2) = STATE_in(2) + k3(2) / 2;
% === RK4 stage 4 ===
t = time_in + tau
DSTATE = f(t,STATE); % calculate the state derivatives
k4(1) = tau*DSTATE(1); % theta
k4(2) = tau*DSTATE(2); % omega
% update the states (Final)
STATE(1) = STATE_in(1) + (k1(1) +2*k2(1) + 2*k3(1) + k4(1))/6;
STATE(2) = STATE_in(2) + (k1(2) +2*k2(2) + 2*k3(2) + k4(2))/6;
% -------------------- RK4 step -----------------------------------------^
% Now we hve a completed time step, and we update the arrays
time(i+1) = time(i) + tau;
theta(i+1) = STATE(1);
omega(i+1) = STATE(2);
end % end of the time loop
The above is the basic technique for the general application of the 4th order Runge-Kutta method. You should be able to see how it can be generalized for systems that have many states by simply implementing loops, looping over the number of states in the STATE vector (in this case, it is just 2). I do a lot of work with 6 degree of freedom flight simulations, and they typically have 13 or more states; Position vector (3 components) velocity vector (3), Attitude quaternion (4), Angular rate vector (3)
Finaly, we need to define the function f() that defines the state derivatives. In general, the state derivatives are a function of time and the system states. In the case of the simple pendulum, there are no time-dependencies, so the system acceleration is simply a function of the theta state. But you should always consider writing code for the general case, so that you have flexibility for future changes.
Your state derivative function will look something like the following:
function DSTATE = f(time,STATE)
theta = STATE(1);
omega = STATE(2);
accel = -g_over_L*sin(theta);
DSTATE(1) = omega;
DSTATE(2) = accel;
return
end
(Note that "f" is probably not a good name for a function - mabe "derive" since it's purpose is to calculate the state derivatives)
This function has all the hooks that you might need to add additional features in the future. For example, suppose that you wanted to add air resistance to the pendulum. Now you can add a calculation for a torque force that is related to the velocity (omega) of the pendulum (this will change accel using a function of omega, so the new accel will be a function of both theta and omega). Since omega is in the state vector, it is already available for this purpose.
Nicholas Davis
on 4 Nov 2020
Jim Riggs
on 4 Nov 2020
"I do have an error with the code not recognizing 'deriv' in DSTATE = deriv(t,STATE);"
There are two ways to include a function in a script file.
1) Include the function code directly in the script.
2) Save the function in a separate file that is on the Matlab search path.
If you save the function to a file, the file name must match the function name, so function "deriv" must be saved in file "deriv.m"
Also, note that I do not have Matlab on my machine that I use to connect to the forum, so the code that I post is typed freehand (not copied) and is untested. I often will have typos in my posted code, most frequently I forget to include the semicolon at the end of the line. As should be obvious, I have not attempted to provide a complete set of code, but enough to explain the concept. For example, I did not include the definition of parameter "g_over_L" in the derivative function - you decide how it gets into the function. You can hard-code it, or add it to the argument list. So, make sure that you understand everything that I have done, and watch out for my typos.
Nicholas Davis
on 4 Nov 2020
Jim Riggs
on 4 Nov 2020
You can plot the results using the saved time history vectors "time", "theta", and "omega"; e.g. "plot(time,theta, ...)" for a plot of theta vs. time.
Nicholas Davis
on 4 Nov 2020
Jim Riggs
on 4 Nov 2020
I just noticed in your code that you posted that you initialize values before pre-allocating the arrays. When you pre-allocate, it wipes out the initial values, (fills the entire vector with zeros). Specifically, you set
theta(1) = theta0*pi/180; and then later
theta = zeros(N+1,1);
so that now theta(1) is zero. So watch out for the order of calculation of your variables. You need N to be defined before you can pre-allocate the vectors. And you need to pre-allocate vectors before you initialize the values.
Nicholas Davis
on 4 Nov 2020
Jim Riggs
on 4 Nov 2020
There is an error in stage 2:
you have:
k2(1) = STATE_in(1)+k2(1)/2;
k2(2) = STATE_in(2)+k2(2)/2;
it should be:
k2(1) = tau*DSTATE(1); %theta
k2(2) = tau*DSTATE(2); %omega
Jim Riggs
on 4 Nov 2020
One more error: In stage 3 you have
STATE(1) = STATE_in(1)+k3(1)/2;
STATE(2) = STATE_in(2)+k3(2)/2;
This should be:
STATE(1) = STATE_in(1)+k3(1);
STATE(2) = STATE_in(2)+k3(2);
Nicholas Davis
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)
More Answers (0)
Categories
Find more on Programming in Help Center and File Exchange
Community Treasure Hunt
Find the treasures in MATLAB Central and discover how the community can help you!
Start Hunting!
