Solving Ordinary Differential Equations
Given Second Order ODE:
Code 1:
%Matlab program to animate a second order ODE
close all;
clear all;
clc;
%Given conditions related to ODE
L=1;
m=1;
b=0.05;
g=9.81;
%Time period is from 0-20 seconds;creating a span of 500 entities for finer
%animation
t_span=linspace(0,20,500);
%1st term indicates angular displacement;2nd term indicates angular
%velocity
theta_initial=[0;3];
%ODE solver
[t,result]=ode45(@(t,theta)ode_func(t,theta,b,m,g,L),t_span,theta_initial);
figure(1);
%result(1) contains all the angular displacements which are plotted against
%time
plot(t,result(:,1));
hold on;
%result(2) contains all the angular velocity which are plotted against time
plot(t,result(:,2));
grid on;
legend('Angular Displacement','Angular Velocity');
xlabel('Time');
ylabel('Angular Displacement/Angular Velocity');
title('Time vs Angular Displacement/Angular Velocity');
ct=1;
for i=1:length(t)
%Fixed point
x0=0;
y0=0;
%Moving end point
x1=L*sin(result(i,1));
y1=-L*cos(result(i,1));
figure(2);
plot([1,-1],[0,0],'linewidth',5);
hold on;
plot([x0,x1],[y0,y1],'linewidth',1.5);
hold on;
plot([0,0],[-5,5],'--');
hold on;
plot(x1,y1,'o','markersize',15,'markerfacecolor','r')
hold off;
axis([-2 2 -1.5 1.5]);
pause(0.02);
M(ct)=getframe(gcf);
ct=ct+1;
end
%Functions to create animation
movie(M)
videofile=VideoWriter('second_ode.avi','uncompressed AVI');
open(videofile);
writeVideo(videofile,M);
close(videofile);
Code 2:
%Function to solve second order ODE by dividing second order ODE into 2
%first order ODE.
function[dtheta_dt]=ode_func(t,theta,b,m,g,L)
dtheta1_dt=theta(2);
dtheta2_dt=-(b/m)*theta(2)-(g/L)*sin(theta(1));
dtheta_dt=[dtheta1_dt;dtheta2_dt];
end
Plot:
Inference:
Hence we can observe that the angular displacement and velocity peak gets lowered as time progresses ,i.e. it gets dampens out ,which is in aggrement with the observations from the physical experiment.
Animation video: