Week 3 - Solving second order ODEs

Step 1- Define function program which will be used in main program to solve the ODE

Step2- Basically a function is a group of statements that together perform a task. In Matlab, functions are defined in separate file. Save the function file in such a way that the name of the file & of the function should be same.

% function program which will be used in main program to solve the ODE
function [dtheta_dt]=ode_func(t,theta,b,g,l,m)
theta1=theta(1)
theta2=theta(2)
dtheta1_dt=theta2;
dtheta2_dt=-(b/m)*theta2-(g/l)*sin(theta1);
dtheta_dt=[dtheta1_dt; dtheta2_dt];
end

Step3- Before starting the code, it’s important to use clear all, close all & clc command.

Step4- Defined the inputs & variables i.e. damping coefficient, gravity, spring length, bob mass.

Step5- Defined initial conditions, like angular displacement, angular velocity & time span

Step6-Solve the ODE, ode45 was used, it is function in Matlab that solve ordinary differential partial equation usng a numeric method. The function ode45 take 3 inputs. The 3 inputs are the function to be integrated, the time span & initial condition. [t, results] means that the results that have been calculated with respect to time

Step7- Plot graph which content two plots. One for angular displacement with respect to time & one for angular velocity with respect to time

Step7-The title, legends & lables are defined.

Step8-Refer attached figure.

Step9-The second part of code with help to create the animation of the pendulum

1. Define ct value
2. Define length of results
3. Define x & start value
4. Calculate end value of both x & y
5. Plot the animation
6. Create the animation movie of pendulum

Refer below codes

% program to solve 2nd order ODE of a pendulam
close all
clear all
clc

% inputs for challenge 
% length for pendulum string = 1 meter
l=1;       
% mass of ball = 1 kg
m=1;        
% daping coefficient = 0.05
b=0.05;      
% gravitational acceleration =9.81 m/s^2
g=9.81;     

% initial condition
% for solving first ODE, required to provide value for disp & velocity
% initial displacement of pendulam = 0 m & initial velocity = 3 m/s
theta_0=[0 3];

% time interval or span for pendulam occilation between 0 to 20 second
t_span=linspace(0,20,100);

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

%ploting
figure(1)
hold on 
plot(t,results(:,1));
plot(t,results(:,2));
title('Damped Harmonic Motion of Pendulum');
legend('Angular Displacement','Angular Velocity');
xlabel('Time(sec)');
ylabel('Plot');

%for loop so before use the for loop counter value was defined
ct=1;
for i=1:length(results(:,1));
    x_start=0;
    y_start=0;
    x_end=l*sin(results(i,1));
    y_end=-1*cos(results(i,1));

    % plot the animation 
    figure(2)
    plot([-1 1],[0 0],'LineWidth',3);
    hold on 
    plot([x_start x_end],[y_start y_end],'LineWidth',3);
    hold on
    plot(x_end,y_end,'o','markersize',20,'markerFacecolor','r');
    title('Animation of Pendulam Movement');
    axis([-1.5 1.5 -1.5 1]);
    pause(0.03);
    hold off
    M(ct)=getframe(gcf);
    ct=ct+1;
end

% Animation of the Pendulam
movie(M);
videofile=VideoWriter('Pendulum Animation.avi','Uncompressed AVI');
open(videofile);
writeVideo(videofile,M);
close(videofile)

Refer Output Figure 1

Refer below you tube link for animation video

https://www.youtube.com/watch?v=ODRRmG-98d8&ab_channel=pravindixit

Error faced during coding - refer attached screen used Ct instead of ct. used capital C instead of small c in ct value representation. After making small c program easly run.