Week 3 - Solving second order ODEs

clear all

close all

clc

% inputs

b=0.05; %Damping coefficent

l=1; %Length of the pendulum

m=1; %Mass of the pendulum

g=9.81; %Gravity

theta_0=[0;3]; %initial_displacement,velocity

t_span = linspace(0,20,500); %time

%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('velocity')

xlabel('time')

% displacement = results(:,1);

% velocity= results(:,2);

ct=1;

figure(2)

for i=1:length(results(:,1));

x0=0;

y0=0;

x1=l*sin(results(i,1));

y1= -l*cos(results(i,1));

figure(2)

plot([-2 2],[0 0])

hold on

plot([x0,x1],[y0,y1],'LineWidth',1,'color','r');

axis([-2 2 -2 2])

plot(x1,y1,'o','MarkerSize',20,'Color','b');

pause(0.01)

M(ct)= getframe(gcf);

hold off

ct=ct+1;

end

movie(M)

videofile=VideoWriter('ode_anime.avi','Uncompressed AVI');

open(videofile)

writeVideo(videofile,M)

close(videofile)

Result video link

Code Explanation:

Main program

the coding/programing starts with the main program with initial values of Damping coefficient, length,mass,gravity,initial velocity, displacement and time span.

ODE45

ode45 if the pre-feeded function in the mat lab which has the detailed sollution of diffirential eqn. sollution is written in the mattrix form so that the reults from this function should be a matrix. Ode45 is the deduced form of a Differential Equation,results of this function are used in main function for plotting and movie making.

in ode_func we are giving the values of b,g,l,m and formulas are written in this function.

Plotting

initially the result of Time and displacement of theta1 is plotted in figure1

then the plot is kept on hold and the results of Theta2 are plotted in same figure with lables of x and y axis

Movie

so for making a movie we have started by creating a loop of i=1 to length of result1 which is displacement values. so that the loop will operate according to the values of result1 (it means loop will repeat it self result1 times)

now giving the cordinates of length of pendulum and its motion equation.

plots are created in figure 2 for pendulum rod and the mass.

M(ct)=getframe(gcf) in this lines we are gathering the frames to creat movie/anime

then the loop is ended