Solving second order ODEs

Steps for the program:

The steps below show to various stages in which the program was executed.

Step 1:

The first step in the process is to use the clear all, close all and clc commands to clear the old data and values and to clear the screen

Step 2:

In this step we need to define the initial values for the Damping Coefficient(b) , Acceleration due to gravity(g), Length of the pendulum(l), mass of the pendulum(m) and radius(r) which is equal to the length of the pendulum. Also, the initial values for displacement and velocity are setup to the variable ‘thetainital’. Also, the timespan for the linkages are set using linspace command, which gives a range of values between 0 and 40 in 500 equal parts.

Also, the origin for the first linkage is given as x0 and y0. The ‘count’ counter is also initiated to a value of 1.

Step 3:

Now we setup a function called ‘odefunction’ having the inputs time(t), angle of pendulum(theta), Damping Coefficient(b) , Acceleration due to gravity(g), Length of the pendulum(l) and mass of the pendulum(m).

Within the function, we map the theta1 to the displacement values of the input theta(in this case thetainitial). Next, we map the theta2 to the velocity values of the input theta(in this case thetainitial). Now, from the analysis we know that

Also, we know that

Now, we code this in words in our program and make an array of the results, and , ie dthetabydt.

Now we link the output of the function to the above-mentioned array.

Step 4:

Now in the main program we input a function call using ode45. Here we require two parameters t and theta. Now within the ode45 function, we link our function odefunction with the parameters for t being timespan and theta being thetainitial. Now the ode45 function computes the differential equation based on the inputs given and gives out the output.

Step 5:

Next we use the plot command to plot a graph between the displacement and the timespan. Also, we plot a graph between the velocity and timespan on the same plot using the hold on command. The legend, axes and the title were setup too using the appropriate commands.

Step 6:

Now in this step, we will be doing the animation for the pendulum.

Firstly, we use a term delta that converts the displacement to the angle using the formula as shown below.

Now we set the limits of the ‘for’ loop from 1 to the size of the ‘delta’ array. Now, we use another variable ‘psi’ to access each values from the array ‘delta’. Now, we can get the values of the end point of the string x1 and y1. For this, we use trigonometry to calculate the instantaneous point for the final coordinates. Now we use the plot command to plot the string of the pendulum. Now similarly, we can plot the pendulum ball using the final points i.e. x1 and y1. We also changed the thickness of the ball marker to make it solid.

Now we setup the axes, legend and the title for the plot. We also, set up the axes limit to trim the plot area more focussed.

Now we use the ‘getframe’ command to collect the frames to make the movie of the animation.

Next we increase the counter and repeat the for loop again. The rest of the program shows the syntax that extracts the frames from the previous command and uses them to make a movie out of it.

clear all
close all
clc

% Inputs

%Damping Coefficient
b=0.05;

%Acceleration due to gravity
g=9.81;

%Length of the pendulum
l=1;

%Mass of the pendulum
m=1;

%Radius=Length of Pendulum String
r=l;

%Starting Point for the String
x0=0;
y0=0;

%Counter Initial value
count=1;

% Initial Values for Displacement and Velocity
thetainital=[0;3];

% Time Span Points from 0 to 20 with 500 values
timespan=linspace(0,20,500);

% Solve ODE using ode45 which takes the inputs time and theta angle
[t, results]= ode45(@(t,theta) odefunction(t,theta,b,g,l,m),timespan, thetainital);


%Plot showing the Displacement and Velocity results against the time span
figure(1)
hold on
plot(t,results(:,1),'linewidth',1)
plot(t,results(:,2),'linewidth',1)
xlabel('Time')
ylabel('Plot')
legend('Displacement Plot','Velocity Plot')
title('Plot versus time span')
hold off

% Animation for the pendulum movement
figure(2)

%Formula for converting the displacement into angle(rads)
delta=results(:,1)/r;

%For loop for reading all values in the array of the results
for i=1:length(delta)
    
    psi=delta(i);
    
    %New points for the end of the string
    x1=r*sin(psi);
    
    %Negative value to denote that the string is downwards
    y1=-r*cos(psi);
    
    %Plotting the string in blue color
    plot([x0 x1],[y0 y1],'linewidth',1,'color','b')
    hold on
    
    %Plotting the pendulum in red
    plot(x1,y1,'O','color','r','linewidth',8)
    
    xlabel('X axis')
    ylabel('Y axis')
    legend('Pendulum String','Pendulum Ball')
    title('Pendulum Animation')
    %Setting the axis limits for the animation
    axis([-1.5 1.5 -1.2 0])

    %Using the given syntax to capture the frames for the movie animation
    M(count)=getframe(gcf);
    
    %Increasing the counter by 1
    count=count+1;
    hold off
    
end

%Syntax for making the movie
movie(M)
videofile=VideoWriter('ODE.avi','Uncompressed AVI');
open(videofile)
writeVideo(videofile,M)
close(videofile)

Error faced:

Error 1:

The plot below was the error that was faced during the programming.

Error 2:

The errors were overcome by observing the code and trying to debug it and also by watching the tutorial video again. The error turned to to be caused due to the omission of square brackets in the plot command.

Animation Youtube Link:

https://youtu.be/4mC4HO8esYg

The plots obtained are as shown below:

Plot 1:

Plot 2: