Week 3 - Solving second order ODEs

Function code:-

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

1) theta 1 & theta2 are where theta is our solution array.

2) dtheta1/dt=theta2 our first ode.

3) dtheta2/dt=-(b/m)*theta2-(g/l)*sin(theta1); Second ode.

Main course:-

clear all
close all
clc
b=0.05
g=9.81
l=1
m=1
%Initial condition
theta_0=[0;3];
%time
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)
plot(t,results(:,1))
hold on
plot(t,results(:,2))
xlabel('Time')
ylabel('plot')

1)Defined :

Damping coeffitient

acceleration due to gravity.

length of pendulum

and mass of the bob.

2) Initial condition

theta_0=[0;3];

initial displacement is 0 and initial velocity is 3.

3)Declared time array taking 500 values between 0-20 sec.

4)Call ODE function using ode45

5)ode_func : function.A separate matlab programme written for defining function

the line -

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

shows ode as function of t and theta and ode_func shows how it is defined.