Week 3 - Solving second order ODEs

Aim:

To write a ode program in matlab for a pendulum occilation animation.

Code:

ode main code:

clear all
close all
clc

b=.05;
l=1;
g=9.81;
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);

figure(1)
plot(t,results(:,1))
hold on
plot(t,results(:,2))
ylabel('plot')
xlabel('time')

angu_dis=results(:,1);

ct=1
for i=1:length(angu_dis)
    ang_dis=angu_dis(i);
    x0=0;
    y0=0;
    
    y1=-l*cos(ang_dis);
    x1=l*sin(ang_dis);
    
    %plotting
    figure(2)
    plot([x0 x1],[y0 y1],'linewidth',2)
    hold on
    plot(x1,y1,'o','markersize',15,'markerfacecolor','r')
    hold off
    axis([-1 1 -1.5 0])
    pause(.01)
    M(ct)=getframe(gcf);
    ct=ct+1;
end

movie(M)
videofile=VideoWriter('pendulem.avi','Uncompressed AVI');
open(videofile)
writeVideo(videofile,M)
close(videofile)

ode_func code:

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

Code explanation:

The given derivation is break downed as dɵ1/dt =ɵ2 , dɵ2/dt=-(b/m)ɵ2-(g/l)sinɵ1. The given values has been written and the initial condition with angular displacement & angular velocity and the time span with intervals has been written.

To solve the ode a function called as ode_func is created and called. This ode_func function contains the above given expression with stating and taking all the given input value and calculates ɵ values for the given expression.

In the main program after entering all the given datas this ode_func function is called to calculed ode using the ode45 inbuilt function and saves the answer in t,results variables.

The graph for angular displacement and angular velocity is plotted as figure(1) with respect to time and plot.

To animate the pendulum the angular displacement is taken from the results and saved in angu_dis, counter(ct) is used to make individual plots. To get different values to plot the animaton we uses for loop to get the values, at first an variable i is set to 1 to length of angular displacement, then the angular displacement of i is set to another variable called ang_dis.

As one end of the pendulum is fixed x0,y0 is set to 0 then for the osscilation end the x1,y1 is set to the multiplication of the pendulum length and sin/cos of ang_dis respectivily.

In animation it should be like pendulum, so line and marker command is used to create the pendulum structure. Then the animation is saved movie in avi format.

Output/Results:

Animation link:

https://drive.google.com/file/d/15Z_YPszj3DEzWyQFFzOxufGPU2fL548x/view?usp=sharing