MATLAB TUTORIAL for the First Cource. Part 4: Pendulum

Prof. Vladimir A. Dobrushkin

This tutorial contains many matlab scripts.
You, as the user, are free to use all codes for your needs, and have the right to distribute this tutorial and refer to this tutorial as long as this tutorial is accredited appropriately. Any comments and/or contributions for this tutorial are welcome; you can send your remarks to

Email Vladimir Dobrushkin

Pendulum

function [t, y] = pendulum(alpha, epsilon, omega0, T, y0)

%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%
%   Usaging the built-in MATLAB function ode45 to solve the second-order 
%   ODE that describes the angular displacement y=theta of an oscillating 
%   pendulum with damping and two initial conditions
%
%   y''(t) + 2a y'(t) + e sign(y'(t)) (y'(t))^2 + omega0^2 sin(y(t)) = 0,
%                                                               0 < t < T
%   y(0) = ya
%   y'(0) = yb
%   
%   The constants in the equation are related to physical properties of the
%   pendulum, its mass and the damping coefficients:
%       damping force F = kappa*|v| + c*v^2
%       mass of bob m
%       length of pendulum string l
%       gravitational acceleration g
%
%   ode45 can only solve first-order problems, so the problem above should
%   be transformed into an equivalent system of two first-order ODEs. By
%   setting 
%       y1 = y      and     y2 = y'
%   we get the system
%       y1' = y2
%       y2' = -2a y2 - e sign(y2) y2^2 - omega0^2 sin(y1)
%   This is the system represented by function F below.
%___ ___ ___ ___ ___ ___ ___ ___ ___ ___ ___ ___ ___ ___ ___ ___ ___ ___ ___ ___ ___ ___ ___ ___ __
%   INPUT : (user must make sure that input values are valid; no check is
%                   performed here)
%           alpha : 0 < alpha = kappa/2m
%           epsilon : 0 < epsilon = cl/m
%           omega0 : 0 < omega0 = sqrt(g/l)
%           T : the maximum time we want to approximate the solution for
%           y0 : a column vector containing the initial values [ya; yb]
%   OUTPUT :
%           t, y : the vectors of timepoints and corresponding 
%               approximations of the solution of the IVP
%
%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%

    F = @(t, y) [y(2); -2*alpha*y(2) - epsilon*sign (y(2)).*(y(2)).^2 - omega0^2*sin(y(1))];
    
    [t, y] = ode45(F, [0, T], y0);
    
end

Next script shows how to use pendulum.m function:

close all
clear all

epsilon = 0;
omega0 = 1;
T = 100;
Colors = hsv(13);
alphas = [0.1, 1, 10];

figure();
for i = 1:3
    
    [t, y] = pendulum(alphas(i), epsilon, omega0, T, [0.7; 0]);
    plot(t, y(:, 1), 'Color', Colors(i, :), 'LineWidth', 1.5);
    hold on
    shg
    
end

legend('0.1', '1', '10', 'Location', 'best');

%--------------------------------------------------------

figure();
alphas = 0.1:0.1:1;
T = 10;
for i = 1:10
    
    [t, y] = pendulum(alphas(i), epsilon, omega0, T, [0.7; 0]);
    plot(t, y(:, 1), 'Color', Colors(i, :), 'LineWidth', 1.5);
    hold on
    shg
    
end

legend('0.1', '0.2', '0.3', '0.4', '0.5', '0.6', '0.7', '0.8', '0.9', '1.0', 'Location', 'best');

Example 1.1.1:

Example 1.1.2:

II. Plotting


                
Separable equations

Enter text here

Equations reducible to separable equations
Exact equations
Integrating Factors
Linear and Bernoulli equations
Riccati equation
Existence and Uniqueness
Qualitative analysis
Applications