Week 4.1 - Solving second order ODEs

Code:

clc;

clear;

close all;

% Given parameters

g = 9.81; % Gravity (m/s^2)

L = 1; % Length of pendulum (m)

m = 1; % Mass of the ball (kg)

b = 0.05; % Damping coefficient

% Initial conditions

theta0 = 0; % Initial angular displacement (rad)

omega0 = 3; % Initial angular velocity (rad/s)

y0 = [theta0; omega0]; % Initial state vector

% Time span

tspan = [0 20];

% Solve the ODE using ode45

[t, y] = ode45(@(t, y) pendulumODE(t, y, g, L, b, m), tspan, y0);

% Extract solutions

theta = y(:, 1);

% Convert to Cartesian coordinates

x = L * sin(theta);

y = -L * cos(theta);

% Animation setup

figure;

hold on;

axis equal;

xlim([-1.2 1.2]);

ylim([-1.2 0.2]);

grid on;

xlabel('X Position (m)');

ylabel('Y Position (m)');

title('Damped Pendulum Motion');

% Plot pivot point

plot(0, 0, 'ko', 'MarkerSize', 8, 'MarkerFaceColor', 'k');

% Create line and marker objects for animation

h_line = line([0 x(1)], [0 y(1)], 'LineWidth', 2, 'Color', 'b');

h_ball = plot(x(1), y(1), 'ro', 'MarkerSize', 10, 'MarkerFaceColor', 'r');

% Animation loop

for i = 1:length(t)

set(h_line, 'XData', [0 x(i)], 'YData', [0 y(i)]);

set(h_ball, 'XData', x(i), 'YData', y(i));

pause(0.05);

end

% Function defining the pendulum ODEs

function dydt = pendulumODE(~, y, g, L, b, m)

theta = y(1);

omega = y(2);

dtheta_dt = omega;

domega_dt = - (b / (m * L)) * omega - (g / L) * sin(theta);

dydt = [dtheta_dt; domega_dt];

end

Step 1: Define the ODE

The equation of motion for a damped pendulum is:

$$mL\frac{d^2\theta}{dt^2} + b\frac{d\theta}{dt} + mg\sin\theta = 0$$

Dividing by $mL$:

$$\frac{d^2\theta}{dt^2} + \frac{b}{mL} \frac{d\theta}{dt} + \frac{g}{L} \sin\theta = 0$$

Define:

$$\omega = \frac{d\theta}{dt}$$

Rewriting as a system of first-order ODEs:

$$\frac{d\theta}{dt} = \omega$$ $$\frac{d\omega}{dt} = -\frac{b}{mL} \omega - \frac{g}{L} \sin\theta$$

Step 2: MATLAB/Octave Code

Here is the MATLAB/Octave program to solve the ODE and animate the motion:

clc;
clear;
close all;

% Given parameters
g = 9.81;  % Gravity (m/s^2)
L = 1;     % Length of pendulum (m)
m = 1;     % Mass of the ball (kg)
b = 0.05;  % Damping coefficient

% Initial conditions
theta0 = 0;       % Initial angular displacement (rad)
omega0 = 3;       % Initial angular velocity (rad/s)
y0 = [theta0; omega0]; % Initial state vector

% Time span
tspan = [0 20];

% Solve the ODE using ode45
[t, y] = ode45(@(t, y) pendulumODE(t, y, g, L, b, m), tspan, y0);

% Extract solutions
theta = y(:, 1);

% Convert to Cartesian coordinates
x = L * sin(theta);
y = -L * cos(theta);

% Animation setup
figure;
hold on;
axis equal;
xlim([-1.2 1.2]);
ylim([-1.2 0.2]);
grid on;
xlabel('X Position (m)');
ylabel('Y Position (m)');
title('Damped Pendulum Motion');

% Plot pivot point
plot(0, 0, 'ko', 'MarkerSize', 8, 'MarkerFaceColor', 'k');

% Create line and marker objects for animation
h_line = line([0 x(1)], [0 y(1)], 'LineWidth', 2, 'Color', 'b');
h_ball = plot(x(1), y(1), 'ro', 'MarkerSize', 10, 'MarkerFaceColor', 'r');

% Animation loop
for i = 1:length(t)
    set(h_line, 'XData', [0 x(i)], 'YData', [0 y(i)]);
    set(h_ball, 'XData', x(i), 'YData', y(i));
    pause(0.05);
end

% Function defining the pendulum ODEs
function dydt = pendulumODE(~, y, g, L, b, m)
    theta = y(1);
    omega = y(2);
    
    dtheta_dt = omega;
    domega_dt = - (b / (m * L)) * omega - (g / L) * sin(theta);
    
    dydt = [dtheta_dt; domega_dt];
end

Step 3: Explanation of Code

1. Define parameters: The gravitational acceleration, pendulum length, mass, and damping coefficient are set.

2. Set initial conditions: The initial angle and angular velocity are defined.

3. Solve ODE using ode45: The equation is rewritten in first-order form and solved using MATLAB’s ode45 function.

4. Convert to Cartesian coordinates: The solution in terms of $\theta$ is converted into x and y coordinates.

5. Animation loop: A figure is created, and the pendulum’s motion is animated using line and plot functions.

Step 4: Common Errors and Fixes

1. "Index exceeds matrix dimensions" error

   • Cause: If the loop index exceeds the number of solution points.

   • Fix: Ensure the loop iterates only up to length(t).

2. "Array indices must be positive integers" error

   • Cause: Using t(i) incorrectly when t is continuous.

   • Fix: Ensure t is indexed correctly with an integer loop.