Week 3 - Solving second order ODEs

MATLAB PROGRAM TO SOLVE THE SECOND ORDER ODE

AIM:

To write the MATLAB program to solve the second order ODE.

MATLAB CODE:

Main Program:

%Program to solve second order ODE
clear all
close all
clc

%Inputs
b=0.05;
g=9.81;
l=1;
m=0.1;

%Initial conditions
theta_0=[0,3];

%Time Period
t_span=linspace(0,20,200);

%ODE 
[t,results]=ode45(@(t,theta)ode_func(t,theta,b,g,l,m),t_span,theta_0);

%Plotting results
 j=1;
for loop=1:length(t_span)
    x0=0;
    y0=0;
    x1=l*sin(results(j,1));
    y1=-l*cos(results(j,1));
    figure()
    line([x0 x1],[y0 y1],'linewidth',2);
    axis([-1.5 1.5 -1.5 0.1]);
    hold on;
    plot(x1,y1,'-o','markers',20,'markerfacecolor','r');
    anim(j)=getframe(gcf);
    j=j+1;
end
movie(anim);
videofile=VideoWriter('Pendulam.avi','Uncompressed AVI');
open(videofile);
writeVideo(videofile,anim);
close(videofile);
 

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

CODE EXPLANATION:

MAIN PROGRAM:

Line 1-4:

Resetting the MATLAB program window.

Line 7-10:

Inputs for the function.

Line 12:

Initial conditions defined as a matrix ‘0’ refers to the initial displacement and ‘5’ is initial velocity in m/s.

Line 16:

t_span=linspace(0,10,100)

Indicates the period for which the function is going to evaluate. Here the functions will be evaluated for 100 values of time between 0 to 10s.

Line 19:

[t,results]=ode45(@(t,theta)ode_func(theta,b,g,l,m),t_span,theta_0);

ode45 is predefined function to solve the differential equations in MATLAB.

@(t,theta) refers that the function depend on t,theta.

t, theta values were assigned by the MATLAB from the values which we sent in t_span and theta_0 array respectively.

ode_func(t,theta,b,g,l,m) will call the function stored in another window.

t_span is the period for which the function is going to evaluate.

theta_0 is the initial conditions for the ODE.

SUB PROGRAM:

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

Line 1:

function [dtheta_dt]=ode_func(t,theta,b,g,l,m)

This is the sub program which has the 2^nd order derivatives broke down into two first order derivatives.

Line 2&3:

theta1=theta(1);

theta2=theta(2);

Here the initial conditions were assigned to the theta1 and the theta2 values.

Line 4&5:

dtheta1_dt=theta2;

dtheta2_dt=-(b/m)*theta2-(g/l)*sin(theta1);

These are the 1^st order derivatives which are derived from the one 2^nd order derivatives.

Line 6:

dtheta_dt=[dtheta1_dt;dtheta2_dt];

In this the dtheta_dt is an array which contains 2x1 array which contains the displacement and the velocity at the time t.

This function is based on the Runge kutta method. Here the ode45 function calculates the values for k1,k2,k3,k4 at the background. For the initial conditions we have provided that displacement is 0 and the initial velocity is 3m/s. Then using the initial conditions and the k1 to k4 values matlab evaluates the displacement and the velocity values for the next step. Now these values were stored to the results array in the main program. Now theta carries the recent values of displacement and velocity. Considering the recent array of theta this function calculates the velocity and the displacement for the further step. This iteration continues until the time limit which we have provided. Now all the values of displacement and the velocity for the time t will be stored in the results array in the main program.

MAIN PROGRAM CONTUATION:

Line 22-35:

j=1;

for loop=1:length(t_span)

    x0=0;

    y0=0;

    x1=1*sin(results(j,1));

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

    figure()

    line([x0 x1],[y0 y1],'linewidth',2);

    axis([-1.5 1.5 -1.5 0.1]);

    hold on;

    plot(x1,y1,'-o','markers',20,'markerfacecolor','r');

    anim(j)=getframe(gcf);

    j=j+1;

end

j=1 is defined to refer the values in the results array and to capture the figure plotted for each values of j. Here the figures were captured using anim variable. After the loop is executed completely the anim variable will have number of frames same as the number of data stored in the t_span array. For example if we have 200 values in t_span array then anim will also will have 200 frames. loop variable here will act as a condition for the for loop to execute. loop variable will allow the for loop to execute from 1 to 200 as t_span contains 200 number of values.

x0,y0 values were assigned to be zero as the one end of the string is fixed. x1, y1 were the end of the string where the ball of mass m will be mounted. In programming to represent the ball we have used the line  “ plot(x1,y1,'-o','markers',20,'markerfacecolor','r');”.x1, y1 values were plotted using the displacement values which are stored in the first column and the row will be provided as variable j as we need to refer all the values in the displacement we will increment the j value at the end of the loop. The coordinates of the lines were shown below.

x0=0;

y0=0;

 x1=1*sin(results(j,1));

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

The line will be plotted using the line “line([x0 x1],[y0 y1],'linewidth',2);” in the program.

The frames which we have captured in the variable anim will be converted in to a video by using the following commands.

movie(anim);

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

open(videofile);

writeVideo(videofile,anim);

close(videofile);

RESULTS:

The animation generated is shown in the links below

Youtube:  https://youtu.be/jMZ-42e_oII

Gdrive: https://drive.google.com/file/d/1Jq1C9lmrITvPBdAzzCRm9WFE76PtSLOx/view?usp=sharing