How do I solve system of differentials using ode45 and eom functions

Hello and thank you for taking the time to read my question!

I am trying to solve a system of differential equations for analytically and then numerically for and and plot the two for comparison. I also have to analytically find tension T in terms of theta then numerically find it in terms of time. The given equations for motion in normal-tangiential form are , and . The starting position is at rad with m=2 kg, g=9.81 m/s^2, and r=0.5 m.

Below is my code so far. I keep getting an error trying to run the ode45 function saying that inputs must be floats, namely a single or double. I have to use the ode45 function as a requirment for the assignment, but I'm not sure if I have to use the eom bit.

Error:

Error using superiorfloat (line 13)

Inputs must be floats,

namely single or double.

Error in odearguments (line 114)

dataType = superiorfloat(t0,y0,f0);

Error in ode45 (line 104)

odearguments(odeIsFuncHandle,odeTreatAsMFile, solver_name, ode, tspan, y0, options, varargin);

Error in ESCI204_M1_Myers (line 74)

[T,S] = ode45(@eom,linspace(0,2,101),[0,0]);

Code:

%Givens:
theta0 = -pi/6;
m = 2;
g = 9.81;
r = 0.5;
syms Tens theta t thetaA
thetaA = linspace(-pi/6,pi/6,101);


% % DETERMINE ANALYTICALLY:------------------------------------------------
%Equation for angular velocity in terms of theta
thetadA = sqrt((-360*g/(pi*r))*(cos(theta0)-cos(thetaA)));

%Equation for Tension in terms of theta
TensA = 3*m*g*cos(thetaA)-2*m*g;

hold on
plot(thetaA,thetadA, 'black')     %analytical solution for thetad(theta)
xlabel('Angular Position (rad) or time (t)')
ylabel('Angular Velocity (rad/sec)')


% % DETERMINE NUMERICALLY:
% 1) The angular velocity and angular position of the mass about point C 
% as a function of time by integration of the equation of motion, thetadd
% 2) The tension in the string as a function of time.

function ds = eom(t,s) %S is current states - vector of current position and vel
%Givens:
theta0 = -pi/6;
m = 2;
g = 9.81;
r = 0.5;
syms Tens

ds(1,1) = sqrt((Tens-m*g*cos(s(1)))/(m*r));   % Derivative of position (s(1) = thetad
ds(2,1) = (-g*sin(s(1)))/r;                   % Derivative of velocity = thetadd
end

[T,S] = ode45(@eom,linspace(0,2,101),[0,0]);

Answer 1

[Aquatris]

You defined a symbolic variable within your ode function, which is a nono if you wanna solve it numerically :D

function ds = eom(t,s) %S is current states - vector of current position and vel
%Givens:
theta0 = -pi/6;
m = 2;
g = 9.81;
r = 0.5;
syms Tens %======= THIS IS A NONO, plug in the values: Tens = 3*m*g*... =======%
ds(1,1) = sqrt((Tens-m*g*cos(s(1)))/(m*r));   % Derivative of position (s(1) = thetad
ds(2,1) = (-g*sin(s(1)))/r;                   % Derivative of velocity = thetadd
end

So plug in your values there

[Steven Lord]

Technically you can use a symbolic variable inside your ODE function, so long as the value returned by your ODE function is numeric. You could do this while still using syms Tens inside eom so long as you used subs at the end of the function to substitute a numeric value for Tens.

More likely than not, since the values I think that @Clara wanted to use in place of Tens were defined in the earlier part of the script, they would need to parameterize the function.

One other consideration is the sizes of the various variables. All the values defined inside eom are scalar, as is the t input. The s input is a 2 element vector. But unless Tens is also scalar (which the TensA variable defined in the script is not) the output of that computation will not be scalar and so can't be assigned into one element of the ds vector.

theta0 = -pi/6;
m = 2;
g = 9.81;
r = 0.5;
thetaA = linspace(-pi/6,pi/6,101);
% % DETERMINE ANALYTICALLY:------------------------------------------------
%Equation for angular velocity in terms of theta
thetadA = sqrt((-360*g/(pi*r))*(cos(theta0)-cos(thetaA)));
%Equation for Tension in terms of theta
TensA = 3*m*g*cos(thetaA)-2*m*g;
whos theta0 m g r thetaA thetadA TensA

  Name         Size             Bytes  Class     Attributes

  TensA        1x101              808  double              
  g            1x1                  8  double              
  m            1x1                  8  double              
  r            1x1                  8  double              
  theta0       1x1                  8  double              
  thetaA       1x101              808  double              
  thetadA      1x101              808  double              

I think based on the equations that the computation of thetadA and TensA should probably move into the eom function, using the value of theta (s(1)) that was passed into it by ode45 to compute those values before they're used in the computation of ds.

[Clara]

@Steven Lord, thank you for your reply! TensA and thetadA were both found on paper and aren't supposed to be used for the numerical solution. They are in the code simply to be plotted on the same plot as thetad, which is what I'm trying to get out of the ode45 function. Tens is the tension on the string at a given point in time, which is one of the other things I need to find for this project. I tried running the script with the Tens function rearranged, but that led to errors with ds not being defined.

I found the analytical solution just using the equation , so I tried commenting out the equation with Tens and that ran, but it made S all 0s. I used a video that my school made, so I'm not sure what S is supposed to be, but I think it is the state, which should be a vector of the position in colomn 1 row 1, and velocity in colomn 1 row 2. The velocity should read similar to, if not exactly, as the black line in the plot created below.

Would it help anything by defining theta as a function of time and everything else as a function of time and theta? The initial state given was theta=-pi/6 at time t=0.

[Sam Chak]

Hi @Clara

The square in this differential equation suggests that there is another branch of solution.

The solution of is in the form of the Jacobi amplitude function.

[Aquatris]

The reason why you are getting 0s as solution is because you defined the initial condition as 0 in your ode with the [0,0] argument. It should be [theta0, 0] instead where theta0 is -pi/6. Also the second argument for ode45 function is the time vector so you should use a time vector with small enough step size, how small depends on the problem.

[t,x] = ode45('odeFunction here','time vec here','initial states here')

Second things is, I am not sure since I do not know the problem exactly, but from what I understand, the equation of motion is defined by the equation only. This completely defines your dynamics so you should solve this one analytically and numerically and compare. In analytical solution, I recommend you use small angle theorem to simplify the equation if it is a pendulum :)

The second equation, is there to calculate the T term after you solve your equation. So after you obtain your analytical and numerical solution, you should compare the resulting analytical T and numeric T by doing a simple .