Week 3 - Solving second order ODEs

Program to solve the ode for pendulum

To solve this I have used ode45 function of MATLAB. ode45 functions takes below parameters
1. ode function which breaks down the ode to first order equation
2. tspan which represent the start and end of integration limits and the size of each steps
3.initial condition

After calling the ode45 with the ode function , a for lopps iterates for the length of t and plots the position of the pendulum . Postion of pendulum is calulated using sine and cosine of theta angle. Pendulum bob is represented using a 'O' and fixed point of pendulum as '*'

Main program:
This program initializes variables , calls ode45 with ode function and do the plotting

%Program to solve the ODE function of pendulum
clear all
close all
clc

%initialize variables
xf = 0;
yf = 0;
b = .05;
m = 1;
l = 1;
g  = 9.81;
count = 1;

% This is initial [displacement, velocity]
theta_0 = [0;3];

%start and end of integration limits and the size of each steps 
tspan = linspace(0,20,250);

%call ode45 function with our ode function and intil parameters and tspan
%integrator
[t,results] = ode45(@(t,theta) Pendulum_odefun(t,theta,b,g,l,m), tspan,theta_0);
figure();
for i = 1 : length(t)
    x = l * sind(results(i));
    y = -l * cosd(results(i));
    
    %plot the fixed point of pendulum for cosmetics and hold it
    plot(0,0,'*','color','b')
    hold on
    
    %plot the pendulum bob movement and hold it
    plot(x,y,'O','color','r')
    hold on
    
    %plot the line joining the fixed point and pendulum bob
    plot([xf x], [yf y], 'linewidth',1);
    hold off  
    
    axis([-1 1 -1 1]) ;
    xlabel('Displacement');

    %capture the frames
    frames(count)  = getframe(gcf);
    count = count + 1;   
    
end

%create the movie
movie(frames);
videoFile = VideoWriter('ODE.avi','Uncompressed AVI');
open(videoFile);
writeVideo(videoFile, frames);
close(videoFile)

ode function
This is the function where ode of pendulum is converted to fir

function [dtheta_dt] = Pendulum_odefun(t,theta,b,g,l,m)
    theta1 = theta(1);
    theta2 = theta(2);
    dtheta1_dt = theta2;
    dtheta2_dt = -(b/m)*theta2 - (g/l)*sind(theta1);
    dtheta_dt = [dtheta1_dt; dtheta2_dt];
end

Result animation

Errors encountered

Sin was used instead of sind in ode function. With this error and initial condition of 0:3 was working fine but when larger displacement or velocity is used the pendulum animiation goes for a toss.