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