I am facing problem with the code which i giving below.

function [theta_dot] = second_order_ode(t,theta,g,l,m);

theta_dot = zeros(2,1);

% Rate of change of angle = angular velocity
theta_dot(1) = theta(2);

% Rate of change of angular velocity = angular acceleration
theta_dot(2) = -(((g/(l*m))^2)*sin(theta));

end

% Solving ODE for Pendulum
 clear all
 close all
  clc
  t = [0 30];
  theta0 = [30 0];
  l = 5;
  m = 2;
  g = 9.8;

    [t, theta] = ode45(@(t,theta)  second_order_ode(t,theta, g, l, m), t, theta0);

    subplot(1,2,1)
    plot(t,theta(:,1),'color','b','linewidth', 2)
    xlabel('Time')
    ylabel('Angular Position')

    subplot(1,2,2)
    plot(t,theta(:,2),'color','r','linewidth', 2)
    xlabel('Time')
    ylabel('Angular velocity')

[James Tursa]

What is the problem? Does the code generate an error? Does it give the wrong answer? Or ...?

[Deepanshu Gautam]

yes it gives the following errors

In an assignment  A(I) = B, the number of elements in B and I must be the same.

Error in second_order_ode (line 9)
theta_dot(2) = -(((g/(l*m))^2)*sin(theta));

Error in @(t,theta)second_order_ode(t,theta,g,l,m)

Error in odearguments (line 87)
f0 = feval(ode,t0,y0,args{:});   % ODE15I sets args{1} to yp0.

Error in ode45 (line 113)
[neq, tspan, ntspan, next, t0, tfinal, tdir, y0, f0, odeArgs, odeFcn, ...

Error in Second_order_ODE_solver_for_Pendulum (line 11)
  [t, theta] = ode45(@(t,theta)  second_order_ode(t,theta, g, l, m), t, theta0);

Answer 1

[James Tursa]

This looks like your initial conditions are in degrees:

theta0 = [30 0];

But the code assumes the input is in radians because of the sin( ) function:

theta_dot(2) = -(((g/(l*m))^2)*sin(theta));

Also, you should be using theta(1) there instead of theta, since theta is a 2-element vector and only the 1st element is the position. So try doing this:

theta0 = [30 0] * (pi/180); % convert from degrees to radians

and this:

theta_dot(2) = -(((g/(l*m))^2)*sin(theta(1)));

And then convert the answer back to degrees for plotting:

    [t, theta] = ode45(@(t,theta)  second_order_ode(t,theta, g, l, m), t, theta0);
    theta = theta * (180/pi); % convert back to degrees

ADVICE: Every time you put constants in your code, add a comment on that same line that specifies what the constant is and what the units of the constant are. E.g.,

t = [0 30]; % sec
theta0 = [30 0]*(pi/180); % degrees (converted to radians)
l = 5; % meters
m = 2; % kg ???
g = 9.8; % meters/sec^2

Finally, this looks like a pendulum problem, but your derivative equation doesn't look correct. You should double check what you have in your code vs what the ODE should be for your problem. (Are you sure that m should be there? Are you sure that ^2 should be there?)

[Deepanshu Gautam]

should i use sind(x) for angles in degree.

[Walter Roberson]

Yes.

[James Tursa]

@Deepanshu: You can work in degrees, but I would strongly advise against this because that will introduce extra conversion constants in your derivative equations. E.g., if theta is in degrees, then your current code for theta_dot(1) and theta_dot(2) will not work because those calculation results on the right hand side were based on using radians for theta in the original ODE. To correct this you would need to introduce unit conversion constants into those lines. Why muck up the derivative equations with that? It is not worth it. Just keep your units as natural radians and save the unit conversion for the plotting stage. So my advice is NO, do not use sind( ) for this ... just use radians for your angle measure in the calling statement and in your derivative equations. Do any unit conversions that you want outside of that code.

To illustrate:

 d(sind(x))/dx
 = d(sin((pi/180)x))/dx
 = (pi/180)cos((pi/180)x)
 = (pi/180)cosd(x)

See how that constant conversion factor got into the derivative equation when using sind( )? Why would you want to muck up your code with that? And you would get another of these constants when you took the derivative of cosd( ). Don't do it ...

Answer 2

[VBBV]

% Solving ODE for Pendulum
 clear all
 close all
  clc
  t = [0 30];
  theta0 = [30 0 0];
  L = 5;
  m = 2;
  g = 9.8;

    [t, theta] = ode45(@(t,theta)  second_order_ode(t,theta, g, L, m), t, theta0);

    subplot(1,2,1)
    plot(t,theta(:,1),'color','b','linewidth', 2)
    xlabel('Time')
    ylabel('Angular Position')

    subplot(1,2,2)
    plot(t,theta(:,2),'color','r','linewidth', 2)
    xlabel('Time')
    ylabel('Angular velocity')

function [theta_dot] = second_order_ode(t,theta,g,L,m);

theta_dot = zeros(3,1);

% Rate of change of angle = angular velocity
theta_dot(1) = theta(1);

% Rate of change of angular velocity = angular acceleration
theta_dot(2) = -(((g/(L*m))^2)*sin(theta(1)*pi/180));   % given one input value here 
theta_dot(3) = -(((g/(L*m))^2)*sin(theta(2)*pi/180));   % given one input value here 

end

Give one input value of theta , it seems you are passing theta as vector of values