Matlab Tutorial for ODE solution: Pendulum Motion

This the tutorial associated with the computer assignment #4. In this tutorial, you will find the Matlab script necessary to solve the nonlinear equations of motion arising from the assigned problem ( You will need to make modifications to the program, however). The problem solved in this tutorial is the hanging pendulum problem (pictured below).

Matlab code

Unlike the vibration code for CA#2, you will need a main code, which sets up the Matlab solution and plotting parameters AND a subprogram which contains the information on the ODE to be solved. Therefore, you need to open up two files (e.g., main.m and subprogram.m):

For CA#4, you will need to change the equations of motion in the subprogram to match the assigned piston problem.

Description of Matlab Code

By modifying the codes above, you can solve the homework problem. In the following section, the pendulum problem is be solved completely, with comments as to the purpose of each command and a plot of the results. The main new concept is how to solve nonlinear (or linear for that matter) ODE's using Matlab.

How Do We Solve ODE's With Matlab?

The magic here is that we don't have to worry too much about solving the nonlinear ODE that governs the motion of the pendulum, we let Matlab do the heavy lifting. How do we solve Ordinary Differential Equations (ODE) using Matlab? Do we have to write lengthy codes? The answer is NO. Life is not that miserable due to the hard work of some numerical methods people. There are several ways in Matlab to solve ODE's, but we are going to discuss only about typical ODE solver in Matlab, ODE45, which is sufficient for our purpose.

So, How Do We Use ODE45 Solver?

ODE45 uses a 4th and 5th order Runge-Kutta method to integrate the equations of motion. Don't worry, you don't have to fully understand 4th and 5th order Runge- Kutta method to use ODE45 solver, check a numerical methods book if you are interested.

ODE45 has a typical form to use. 

tspan = [t_start, t_final];
X0 = [x0, xdot0];

options = odeset('Refine',6,'RelTol',1e-4,'AbsTol',1e-7); 
  
%  options 'Refine' was used to produce good-looking plots
%  default value of Refine was 4.
%  Refine doesn't apply if length(TSPAN) > 2.
%  scalar relative error tolerance 'RelTol' (1e-3 by default) and 
%  vector of absolute error tolerances 'AbsTol' (all components 1e-6 by
%  default).

[T, Y] = ODE45(F, tspan, X0, options);
T in the left bracket is for time, x is our main variables, in vector form. For example, if we want to know the displacement of a pendulum and use "x" as a rotation angle, "x" will be the main variable in the equations of motion. Now, F in the right side is a name of the subprogram for ODE45 solver (will be explained in detail a little bit later). tspan is the time period that you want to calculate, which consists of start time and end time. X0 is the initial conditions for ODE variable, x0 = x(0), xdot0 = xdot(0). options is the command that you can adjust the performance of ODE45 solver. There are many options you can choose from, type 'help odeset' at the Matlab command line. 'Refine' controls the number of data points for ODE solver. Higher the number after 'Refine' is more data points there are, but slower the calculation is. 'RelTol' controls the relative error tolerance and 'AbsTol' controls the absolute error tolerances.
Now, the tricky part is how to use a subprogram right. An example will be the best way to learn how to do it right.

Source Code

The complete source code to run this problem can be found by clicking  here.  Note that you must save the two parts separately and further, the sub-program must be saved as pendulum_eom.m in this example.  In your computer assignment you will have some different filename which must appear in the same three places that "pendulum_eom" is seen in the example codes.

Further Help

If you need some more information about commands, try "help". For example, 

       help plot