Challenge

Week 3 - Solving second order ODEs

AIM
To write a programme that solves the following ODE.This ODE represents the equation of motion of a simple pendulam with damping.
(ⅆ^2 θ)/(ⅆt^2 )+b/m  ⅆθ/ⅆt+g/L  sin⁡θ =0
Where g=gravity in m/s^2
             L=length of pendulam=1meter
             m=mass of pendulam in kg=1
             b=damping coefficient=.05
solving the ODE
let θ=θ_1
so that ,    ⅆθ/ⅆt=(ⅆθ_1)/ⅆt=θ_2------------(1)
therefore              (ⅆ^2 θ)/(ⅆt^2 )=(ⅆ^2 θ_1)/(ⅆt^2 )=ⅆ/ⅆt ((dθ_1)/dt)=ⅆ/ⅆt (θ_2 )--------(2)
substituting values of equation 1 and 2 in the main ODE
(ⅆ^2 θ1)/(ⅆt^2 )+b/m  ⅆθ1/ⅆt+g/L  sin⁡θ1 =0
Equal to,  (ⅆ^ θ2)/ⅆt+b/m  θ2/+g/L  sin⁡θ1 =0
  So that,          (ⅆ^ θ2)/ⅆt=-b/m  θ2/-g/L  sin⁡θ1
So in matrix form we can represent the solution of above equations as below,
ⅆ/ⅆt [█(θ_1@θ2)]=[█(θ_2@-b/m θ_2-(g/l 〖)sin〗⁡θ )]
So making the above solving  to a matlab function 
function [dtheta_dt]=ode_func(t,theta,b,g,l,m)
 
 
theta1=theta(1);
theta2=theta(2);
dtheta1_dt=theta(2);
dtheta2_dt=-(b/m)*theta2-(g/l)*sin(theta1);
dtheta_dt=[dtheta1_dt;dtheta2_dt];
 end
 

Then the programme for simple pendulam animation is as follows using the above function


clear all
close all
clc
b=0.05;
g=9.81;
l=1;
m=1;
%initial condition
theta_0=[0;3];
%time points
t_span=linspace(0,20,500);
%solve ODE
[t,results]=ode45(@(t,theta)ode_func(t,theta,b,g,l,m),t_span,theta_0);
ang_disp=results(:,1);
ct=1;
for i=1:length(ang_disp)
    x0=0;                                                      
    y0=0;  
    x(i)=l*sin(ang_disp(i));
    y(i)=-l*cos(ang_disp(i));
%     plotting
    plot([x0,x(i)],[y0,y(i)],'linewidth',3);
    axis([-1.5 1.5 -1.2 .2]);
    ylabel('plot');
    xlabel('time');
    
    hold on
    plot(x(i),y(i),'o','MarkerFaceColor','red'); 
    hold off
    M(ct)=getframe(gcf);
    ct=ct+1;                                                                                                                                          
    pause(0.3);                                 
end
    movie(M)
videofile=VideoWriter('simple_pendulam01a.avi','uncompressed AVI');
open(videofile)
writeVideo(videofile,M)
close(videofile)
 




Your portfolio, next

Ship work like this
in your own course.

Talk to a career counsellor

Free 20-minute call. We'll map the right course to where you want to work.

By submitting you agree to our Terms and Privacy Policy.

Talk to a career counsellor

Free 20-minute call. We'll map the right course to where you want to work.

By submitting you agree to our Terms and Privacy Policy.