Solving Second order differential equations using MATLAB

Introduction

The presented project report attempts to provide a comprehensive study and develop a script for solving second order differential equations for the transient motion using MATLAB. MATLAB is software featuring a high-level programming language, primarily intended for numerical computations. Octave helps in solving linear and nonlinear problems numerically, and for performing other numerical experiments to describe the transient behavior of a system.

A simple example is a pendulum. The way the pendulum moves depends on the Newton's second law. When this law is written down, we get a second order Ordinary Differential Equation that describes the position of the "ball" w.r.t time.

The simple pendulum

The simple pendulum consists of a point mass m hanging on a string of length L from some support. One pulls the mass back to some starting angle, θ0, and releases it. The goal is to find the angular position as a function of time.

We could either use Newton’s Second Law of Motion, F = ma, or its rotational analogue in terms of torque, τ = Iα. We will use the former only to limit the amount of physics background needed. There are two forces acting on the point mass. The first is gravity. This points downward and has a magnitude of mg, where g is the standard symbol for the acceleration due to gravity. The other force is the tension in the string. In these forces and their sum are shown. The magnitude of the sum is easily found as F = mg sin θ using the addition of these two vectors.

Now, Newton’s Second Law of Motion tells us that the net force is the mass times the acceleration. So, we can write mx¨ = −mg sin θ. Next, we need to relate x and θ. x is the distance traveled, which is the length of the arc traced out by the point mass.

This ODE represents the equation of motion of a simple pendulum with damping.

Given Parameters:

g (m/s2)	gravity	9.81 m/s2
m (kg)	mass of the ball	1 kg
b	damping coefficient	0.05
L (meter)	The length of the pendulum	1 meter
Time Interval	Simulation motion	0-20 sec

angular displacement=0,

angular velocity=3 rad/sec at time t=0.

Written Script:

%Damping coefficent = b
b = 0.05;
%gravity = g m/s^2
g = 9.81;
%length of pendulum string = l metre
l = 1;
%mass in kg = m
m = 1;

%Intial condition 
%for solving first ode it is required to provide value for displacement and velocity
% initial displacement of the pendulum is 0 
%intial velocity is 3m/s

theta_0 = [0 3];

%time intervals
t_span = linspace(0,20,100);
% solving ODE
[t, results] = ode45(@(t,theta) ode_func(t, theta, b,g,l,m),t_span,theta_0);
outcome = results(:,1); 
%contor_command
ct = 1;
for i= 1:length(outcome)
    pos = outcome(i);
    %origin_points
    x0 = 0; 
    y0 = 0;
    %max positions
    x1 = l*sin(pos);
    y1 = -l*cos(pos);
    %plotting
    plot([x0 x1], [y0 y1],x1,y1,'g.','markersize',15,'linewidth',2)
    axis([-1.5 1.5 -1.5 1.5])
    hold off
    grid on
    M(ct) = getframe(gcf);
    ct = ct+1;
end
%Animation
movie(M)
videofile = VideoWriter('second_order_ pendulum');
open(videofile)
writeVideo(videofile,M)
close(videofile)

Script to produce damping graph.

%Plotting
figure(1)
plot(t,results(:,1))
hold on;
xlabel('Time')
plot(t,results(:,2))
hold off;
axis([0 22 -3.5 3.5])
ylabel('plot')
legend('Angular displacement/time','Angular velocity/time')

Script explanation:

Defining Parameters

%Damping coefficent = b

%gravity = g m/s^2

%length of pendulum string = l metre

%mass in kg = m

b = 0.05;

g = 9.81;

l = 1;

m = 1;

%Intial condition

%for solving first ode it is required to provide value for displacement and velocity

%initial displacement of pendulum is 0

%initial velocity is 3m/s

theta_0 = [0 3];

In above code we just defined the given parameters.

%time intervals

t_span = linspace(0,20,100);

% solving ODE

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

outcome = results(:,1);

Time interval is the time span for which our pendulum shows transient motion and generating results by using command ode45 which is an inbuilt function for solving differential equations. outcome is nothing but the results following 1st column.

Ode_func command here is hyperlinked code written to solve differential equations.

Ode_func script:

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

here we just solved the 2^nd order differential equations and encoded them into code.

%contor_command

ct = 1;

for i= 1:length(outcome)

    pos = outcome(i);

    %origin_points

    x0 = 0;

    y0 = 0;

    %max positions

    x1 = l*sin(pos);

    y1 = -l*cos(pos);

    %plotting

    plot([x0 x1], [y0 y1],x1,y1,'g.','markersize',15,'linewidth',2)

    axis([-1.5 1.5 -1.5 1.5])

    hold off

    grid on

    M(ct) = getframe(gcf);

    ct = ct+1;

end

contor command is used to change loop irretation after every motion. X = 0 and Y = 0 depicts initial position of pendulum bob and using trigonometry function we can generate transient motion of the pendulum and maximum positions of the pendulum motion. Using the marker or viscircles command we can attach a bob to pendulum string.

%Animation

movie(M)

videofile = VideoWriter('second_order_ pendulum');

open(videofile)

writeVideo(videofile,M)

close(videofile)

Command for generating animation files to show transient behavior of pendulum.

Thank You