Week 3 - Solving second order ODEs

Solution of ODE:

theta=theta1

dtheta/dt=dtheta1/dt=theta2

d^2theta/dt^2=d^2theta1/dt^2=dtheta2/dt

by this ,

dtheta1/dt=theta2

dtheta2/dt=-(b/m)theta2-(g/l)sin(theta)

solution for theta

d/dt[theta1;theta2]=[theta2;-(b/m)theta2-(g/l)sin(theta)]

Input:

For this program the function ode_func is first created

function [dtheta_dt] = ode_func(t,theta,b,m,g,L) % creating function named ode_func
  theta1=theta(1); % Assign value theta(1) to theta1
  theta2=theta(2); % Assign value theta(2) to theta2
  
  % find derivatio of theta1 and theta2
  dtheta1_dt = theta(2); 
  dtheta2_dt = -(b/m)*theta2-(g/L)*sin(theta1);
  
  % creating array of dtheta_dt for function which is coloumn matrix
  [dtheta_dt] = [dtheta1_dt; dtheta2_dt];
  
end % end function

By this function i creat program for solve ode and animation for pendulum

% Solving 2nd order ode for pendulum 
clear all
close all
clc

% Input data
L = 1    % meter   
m = 1   % kg
b = 0.05
g = 9.81 % m/s^2

% time span for pendulum
t_span = linspace(0,20,500); % total time taken is 20s where 500 inbetween values are considered

% Intial condition
theta_0 = [0;3]  % 1st coloumn represents displacement while 2nd represents velocity

% Solving ode by ode45
[t,results] = ode45(@(t,theta) ode_func(t,theta,b,m,g,L),t_span,theta_0); 

% Plotting for analysing pandulum transient behavior and pendulum animation

figure(1) % Creating figure window for plotting
plot(t,results(:,1)) % Plotting a graph for time vs results(:,1) which is displacement
hold on              % Hold the figure window 1 for furture plotting
plot(t,results(:,2)) % Plotting a graph for time vs results(:,2) which is velocity
xlabel('Time')                  % Named x-axis as Time
ylabel('Displacement/Velocity') % Named x-axis as Displacement/Velocity
grid on  % make grid visible
hold off % Hold off the figure window 2

ct=1; % Assign value 1 to counter
for i=1:length(results(:,1)) % Assign velue to i 
  
  % Origin or starting position of pendulum which is fixed
  x0=0;
  y0=0;
  
  % Varied value for end point of pendulum
  x1=L*sin(results(i,1));
  y1=-L*cos(results(i,1));
  
  figure(2) % Creating figure window for plotting an animation of pendulum 
  hold on   % hold the figure window 2
  
  plot([-1 1],[0 0],'linewidth',8);     % Plotting fixed roof for pendulum's fixed point
  plot([x0 x1],[y0 y1],'linewidth',3,); % Plotting link or thred by which pendulum's ball is held
  plot(x1,y1,'marker','o','markersize',30,'markerfacecolor','b');     % Plotting the ball of pendulum
  plot(x0,y0,'marker','o','markersize',9,'markerfacecolor',[1 1 1]);  % Plotting the fixed point of pendulum
  line([x0 x0],[y0 -1.5],'linestyle','--','color','r','linewidth',3); % Plotting of centre line 
  
  axis('tight');  % To adjust axis as our requriments
  pause(0.00003); % update the figure at interval of 0.00003s
  
  grid on % Make the grid visible
  
  text = sprintf('test%d.png',ct); % Creat .png file for each image created by plotting
  saveas(gca,text);                % save all .png file 
  
  clf();   % Clearing the current graphic image 
  
  ct=ct+1; % Increment in figure counter
 

end % end for loop

Output:

Error and solution:

there are many errors occurs while generating loop for animation.

The main cause of error i faced is a typo which was "results(i:1)" which is not correct.

the correct form for this program is "results(i,1)".

I faced many more error because typo's in creating animation in octave.