Week 3 - Solving second order ODEs

https://drive.google.com/file/d/1hfGQ8GIQXEsU3mWDS17-3fnxsWsO1ohw/view?usp=sharing

%function is defined
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

clear all
close all
clc
% input parameters
b=0.02;%damping coefficient in Ns/m
g=9.81;%accelaration off gravityon m/s^2
l=1;%length of wire or rod i m
m=0.1;% mass of the ball in kg

%parameters to plot graph of displacment , velocity, vs time 

theta_0=[0;3];% inttial displace ment and angular velocity
%linspace(X1, X2, N) generates N points between X1 and X2
t_span=linspace(0,20,150);  % time duration fore which the simulation is done 
% solving given problem using ode45 
[t,x]=ode45(@(t,theta) ode_func(t,theta,b,g,l,m),t_span,theta_0);
 % here t is time and c contain 2 column ie column 1 contain displacement and column
 %and column 2 is angular velocity which is defined in ode_func named
 %function 
 figure(1) % figure 1 to show the time vs displacement plot   
plot(t,x(:,1));%plot of time vs displacement here x(:,1) means data of x matrix of column 1 is extracted
hold on;% to hold the current plot
plot(t,x(:,2))%plot of time vs angular velocity here x(:,2) means data of x matrix of column 2 is extracted
legend('displacement','velocity')%to differntial b/w displacement plot and angular velocity plot
title('spring mass system')%for title
ylabel('displacemen & velocity')%y axis label
xlabel('time')%x axis label


% to generating animation of simple pendulum 
ct=1;%count
for i=1:length(t);
    x0=0;
    y0=0;
    
  p1=l*sin(x(i,1));%to locate the wire or rod end of pendulum of x axis
  p2=-l*cos(x(i,1));%to locate the wire or rod end of pendulum of y axis
  
%   to plot animation of pendulum
figure(2)
  plot([-0.5,0.5],[0,0],'linewidth',10,'color','c');%to plot ceiling of pendulum of line width of 10 and color of cyan
  
 hold on;% to hold current plot
 grid on;% create grid on current plot
 
  plot([x0 p1],[y0 p2],'linewidth',3,'color','b');%to plot wire or rod of pendulum
  
  
  plot(p1,p2,'o','MarkerSize',20,'markerfacecolor','r');%to create sphere
 v=[-1.5 1.5 -1.5 0.5];
 axis(v)%to fix the axis at vector v
 pause(0.001)
 M(ct)=getframe(gcf);%getframe returns a movie frame. The frame is a snapshot
%     of the current axes
 hold off%torelease helpd plot

 
 ct=ct+1;

end
movie(M);
videofile=VideoWriter('simple_pendulum.AVI');%to create video file named simple_animation
open(videofile);
writeVideo(videofile,M)%to write the video with frame M
close(videofile)

Q

solve this eqn using matlab ode function `(d^2(theta))/dt^2+(d(theta))/dt+g*sin(theta)/l=0`

soln:-