Week 3 - Solving second order ODEs

Objective:

To write a program that solves the following ODE. This ODE represents the equation of motion of a simple pendulum with damping

In the above equation,

g = gravity in m/s2.

L = length of the pendulum in m.

m = mass of the ball in kg.

b=damping coefficient.

ODE's:

An odinary differential equation is a differential equation containing one or more function of one independent variable and the derivaties of the function. These equations are used to calculate the movement of electricity ,motion of an object i.e, a simple pendulam.

Input Values:

L=1 metre.

m=1 kg.

b=0.05.

g=9.81 m/s2.

Code for 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

Code Explanation:

First of all every second or higher degree differential equation should be break down into several first order equations.

θ = θ1

dθ/dt = dθ1/dt = θ2.........(1)

dθ2/dt = -b/m*θ2 - g/l sinθ1 ........(2)

This how we break a first degree equation form a second degree equation.

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

We are creating a function ode_func(t,theta,b,g,l,m) and assigning it to a variable called [dtheta_dt]. The inputs are t,theta,b,g,l,m.

Code for simple pendulam:

clear all
close all
clc

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

%Initial condition
theta_0 = [0;3];

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

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

%Plotting
plot(t,result(:,1))
hold on
plot(t,result(:,2))
xlabel('Time')
ylabel('Plot')

%Animation
x0=0
y0=0
x1=-2
y1=0
x2=2
y2=0
ct=1
for i = 1:length(t)
    x=l*sin(result(i))
    y=-l*cos(result(i))

    plot([x1,x0,x2],[y1,y0,y2],'linewidth',2)
    hold on
    plot([x0,x],[y0,y],'linewidth',2)
    hold on
    plot(x,y,'o','markersize',15,'markerfacecolor','b')
    axis([-3 3 -3 0])
    hold off
    pause(0.1)
    d(ct) = getframe(gcf);
    ct=ct+1;
end
movie(d)
videofile=VideoWriter("Pendulam.avi","Uncompressed AVI");
open(videofile)
writeVideo(videofile,d)
close(videofile)

Code Explanation:

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

First the inputs are given for the pendulam.

%Initial condition
theta_0 = [0;3];

Initial conditions are given for the pendulam, where 0 is the starting time and 3 is the initial velocity.

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

The time range is given by the linspace command for 20 sec with 200 values in between.

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

Now, to solve the ode the solver ode45 and the function ode_func is used along with the input values.

%Plotting
plot(t,result(:,1))
hold on
plot(t,result(:,2))
xlabel('Time')
ylabel('Plot')

To display the output plot command is used for the velocity and time. The title for the x and y axis is given using the xlabel and y label.

%Animation
x0=0
y0=0
x1=-2
y1=0
x2=2
y2=0
ct=1
for i = 1:length(t)
x=l*sin(result(i))
y=-l*cos(result(i))

plot([x1,x0,x2],[y1,y0,y2],'linewidth',2)
hold on
plot([x0,x],[y0,y],'linewidth',2)
hold on
plot(x,y,'o','markersize',15,'markerfacecolor','b')
axis([-3 3 -3 0])
hold off
pause(0.1)
d(ct) = getframe(gcf);
ct=ct+1;

Now for makind the animation the values of x0,x1,x2 and y0,y1,y2 are assigned.

For loop is used for looping process.

To plot the values plot command is used for the previously assigned values.

The size of the ball is created by the mareker size command and the colour of the ball is given as blue.

The axis values are given using axis command.

To get the frames for the animation getframe command is used.

The value of counter is assigned as 1 previously in the inputs section here the increment is given as +1 for the loop process.

movie(d)
videofile=VideoWriter("Pendulam.avi","Uncompressed AVI");
open(videofile)
writeVideo(videofile,d)
close(videofile)

To save the animation as a video and save it in the directory videofile=VideoWriter("Pendulam.avi","Uncompressed AVI"); this command is used.open(videofile) is used to open the video file.

Output:

The output obtained from the program for velocity and time.

Errors:

Forgot to add the value for counter in for loop. Got this error message and then rectified it.

In solving the ode function i forgot to give the initial values goe time and velocity which is theta_0.

Link for the video:

https://youtu.be/nYDyL51UVKI

Conclusion:

Thus the following ODE is solved.