Week 3 - Solving second order ODEs

Challenge 4, Solving second-order ODEs and animating a simple pendulum

Report:

The program starts with providing input data for the simple pendulum such as damping coefficient, the mass of the ball, acceleration due to gravity and length of the pendulum.

The following inputs were provided,

Damping coefficient, b = 0.05

Mass of the ball, m = 1 kg

Acceleration due to gravity, g = 9.81 m/s^2

Length of the pendulum, l = 1 m

The initial condition for the ODE was specified next,

Starting displacement = 0 rad

Initial angular velocity = 3 rad/s

The ODE was solved for a time span ranging from 0-20 seconds as specified in the input code for the variable t_span.

Next, the ODE was solved by calling the ode45 solver with the independent variables t and theta and by passing the function ‘ode_func’ with the input variables t, theta, b, n, g and l.

Furthermore, the integration span ‘t_span’ was specified with the initial condition stored in the variable ‘theta_0’.

Code for animation was specified next.

The fixed points for the pendulum were specified with the coordinates x_0 and y_0.

The theta values ranging from 0 to 20 were converted into degrees and stored in the variable named disp_values.

Each still of the plot was fed into a variable named M. subsequently a movie was made using the command Movie(M).

Lastly, a video file was created, opened, named pendulum and finally closed using the commands open(), writeVideo(), close()

Inside the function ‘ode_func’,

The second-order ODE was simplified into a system of first-order ODEs, solved and sent back to the main program.

Main Program

% challenge 4, ODE for second order equation
clear all;
close all;
clc;

%inputs
b=0.05; %damping
m=1; %mass in kg
g=9.81; %acclereation due to gravity m/s^2
l=1; %length of pendulum

%initial condition
theta_0=[0; 3];

%time span for integration
t_span= linspace(0,20,100);


%solving ODE
[t, results]= ode45(@(t, theta) ode_func(t, theta, b, m, g, l), t_span, theta_0);

%animation

%fixed co ordinate of the pendulum
x_0=0;
y_0=0;
disp_values=results(:,1).*(180/pi);%converting radians to degrees
ct=1;

for i=1:length(disp_values)
    current_theta= disp_values(i);
    
    %end point co ordinate of the pendulum
    x_1=-l*sind(current_theta);
    y_1=-l*cosd(current_theta);

    %plotting
    plot([x_0 x_1],[y_0 y_1], 'color', 'r','LineWidth', 3 )
    hold on
    plot(x_1,y_1,'.','MarkerSize',50,'color','g');
    axis([-1 1.2 -1 1.2]);
    grid on
    pause(0.05);
    M(ct)=getframe(gcf);
    ct=ct+1;
    hold off
end

%making animation video file
movie(M);%has default frame rate of 12
videofile = VideoWriter('pendulum','Uncompressed AVI');
videofile.FrameRate = 12;%reducing the default frame rate from 30 to 12 to mtch movie command
open(videofile);%once opened frame rate cannot be adjusted
writeVideo(videofile,M);
close(videofile);
    

Function

% function for solving ODE

function [dtheta_dt] = ode_func(t, theta, b, m, g, l)
theta_1=theta(1);
theta_2=theta(2);

%solving equations/ first order ODE's

dtheta_dt_1= theta(2);
dtheta_dt_2= -(b/m)*theta_2-(g/l)*sin(theta_1);
dtheta_dt= [dtheta_dt_1;dtheta_dt_2];
end

Animation Link : https://drive.google.com/file/d/1njjEvstL4cKVviCPphww87i9jKtmq1FY/view?usp=sharing

Errors faced

• An incorrect variable name was given to adjust the video framerate.

Correction: variable name was typed correctly

• The incorrect syntax used in plot command with marker indices

Correction: syntax was corrected

• The incorrect syntax used with axis command

Correction: syntax was corrected

• Incorrect variable name typo

Correction: variable name was named correctly

• Forgot to mention the initial condition for the ode45 solver

Correction: initial condition was specified