Week 3 - Solving second order ODEs

Given equation

`(d^2 theta)/dt^2 +b/m(d theta)/dt+g/lsin theta =0`

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.

The given equation is second order but the function ode45 can solve only first order equation so convert this second order equation to two first order equation , for this

substitute

`theta = theta_1 and (d theta)/dt=theta_2`

form the solution vector

`(d)/dt[(theta_1) ,(theta_2)]=[[theta_2],[-b/mtheta_2-g/lsin theta_1]]`

ode45

MATLAB’s standard solver for ordinary differential equations (ODEs) is the function ode45. This function implements a Runge-Kutta method with a variable time step for efficient computation.The important thing to remember is that ode45 can only solve a first order ODE. Therefore to solve a higher order ODE, the ODE has to be first converted to a set of first order ODE’s

[t,y] = ode45(odefunction,tspan,y0)

where ,

odefunction - Function to solve. The function dydt=odefunction(t,y), for a scalar t and a column vector y, must return a column vector dydt that corresponds to f(t,y).

tspan = [t0 tf], integrates the system of differential equations y′=f(t,y) from to to tf with initial conditions y0

Each row in the solution array y corresponds to a value returned in column vector t

Write code

1)Define given values , inital condition and time span to variable and call ode45 function

2) write a function pendulum_ode to form solution vector

3) Make subplot for displacement vs time , velocity vs time and distance vs velocity with proper label and title

3) Add points of different graph for each value of time using for loop

4) plot the position of pendulum for different values of time

5) Pause the program and save the current frame to struct 'M' and then increment the counter 'ct'

6) Save the animation to a video file

Errors faced while programming.

1)

solution:- Removeing the comma ',' before calling function pendulum_ode

2)

pendulum is oscillating upwards solve by adding negative sign at line 80

correct line is 'y1=-l*cos(result(i,1));'

3)

solve by using correct built-in function name VideoWriter

4)

Error in profile name solve by going to 'help videoWriter' and find the correct profile name which is 'Uncompressed AVI'

Final code after solving errors

Script pendulum

close all
clear all
clc

% given values
b=0.05;  %damping coefficient
g=9.82;  %gravity in m/s2
l=1;     %length of the pendulum in m
m=0.1;   % mass of the ball in kg
 
%inital condition
theta_0=[0 3];

% time points
t_span=linspace(0,20,200);


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


% subplot for displacement vs time
subplot(2,2,3)
a=animatedline;   % adding the animated line to the subplot
xlabel('time')         % adding label
ylabel('angular position')    % adding label
title('Displacement vs time graph') %adding title
grid on


% subplot for velocity vs time
subplot(2,2,4)
b=animatedline;  % adding the animated line to the subplot
xlabel('time')     % adding label
ylabel('angular velocity')   % adding label
title('velocity vs time graph') %adding title
grid on 



% subplot for distance vs velocity
subplot(2,2,2);
h=animatedline;  % adding the animated line to the subplot
xlabel('position in radian')   % adding label
ylabel('angular velocity')    % adding label
title('Displacement vs velocity graph')  %adding title
grid on

  


ct=1;  %frame count 1
for i=1:length(t)
    
    % adding points in displacement vs velocity curve
    subplot(2,2,2);
    addpoints(h,result(i,1),result(i,2))
    axis([-1 1 -3 3 ])   %adjusting axis 
      
    
     % adding points in displacement curve
    subplot(2,2,3)
    addpoints(a,t(i),result(i,1))
    axis([0 20 -3 4 ])   %adjusting axis 
    
    % adding points in velocity curve
    subplot(2,2,4)
    addpoints(b,t(i),result(i,2))
    axis([0 20 -3 4 ])  %adjusting axis 
    
    
    % plot the pendulum 
    
    % fixed positionof the pendulum
    x0=0;  
    y0=0;
    
    % moving ball position
    x1=l*sin(result(i,1));
    y1=-l*cos(result(i,1));
    
    % ploting
    subplot(2,2,1)
    
    % plot a fixed suface 
    plot([-1 1],[0 0],'LineWidth',4,'color','k')
    axis([-1 1 -1.5 0.5])
    hold on 
    % plot string
    line([x0 x1],[y0 y1],'LineWidth',2,'color','b')
    % plot the ball
    plot(x1,y1,'-o','markers',10,'markerfacecolor','r')
    title('Pendulum')
    hold off
    
    %pause the program
    pause(0.001)
    
    % save the current frame to M
    M(ct)=getframe(gcf);
    
    %increase frame counter
    ct=ct+1;
end


% make movie from the frames
movie(M);

%create videofile
videofile=VideoWriter('pendulum.avi','Uncompressed AVI');

%open file
open(videofile)

%write movie of M in videofile
writeVideo(videofile,M)


%colse file
close(videofile)

Function pendulum_ode

function dtheta_dt=pendulum_ode(theta,b,g,l,m)
    theta1=theta(1);
    theta2=theta(2);
    
   % put theta=theta1 and dtheta_dt=theta2 to get below equation
    dtheta1_dt=theta2;
    dtheta2_dt=-(b/m)*theta2-(g/l)*sin(theta1);
    % make the solution column vector
    dtheta_dt=[dtheta1_dt;dtheta2_dt];
    
end

Result

https://youtu.be/zJjiTHmAMs0